Get geographic distances

distance.py

# A Pythonofication of http://www.bdcc.co.uk/Gmaps/BdccGeo.js
# Originally in Javascript by Bill Chadwick

import sys
import math
import re
from urllib import urlopen
import json
from pprint import pprint
import plistlib

class Coordinates(object):
    F = (1.0 - 1.0 / 298.257223563)**2
    EQ_RADIUS = 6378137.0; # WGS84 Equatorial Radius in Meters
    def __init__(self, lon, lat):
        theta = lon * math.pi / 180.0
        rlat = self.geocentric_latitude(lat * math.pi / 180.0)
        c = math.cos(rlat)
        self.x = c * math.cos(theta)
        self.y = c * math.sin(theta)
        self.z = math.sin(rlat)

    # Class methods

    @classmethod
    def geocentric_latitude(cls, geographic_latitude):
        return math.atan(math.tan(geographic_latitude) * cls.F)

    @classmethod
    def geographic_latitude(cls, geocentric_latitude):
        return math.atan(math.tan(geocentric_latitude) / cls.F)

    @classmethod
    def get_intersection(cls, geo1, geo2, geo3, geo4):
        '''Returns the two antipodal points of intersection of two great circles
        defined by the arcs geo1 to geo2 and geo3 to geo4. Returns a point as a Geo,
        use .antipode to get the other point'''
        geo_cross1 = geo1.cross_normalize(geo2)
        geo_cross2 = geo3.cross_normalize(geo4)
        return geo_cross1.cross_normalize(geo_cross2)

    @classmethod
    def radians_to_meters(cls, rad):
        return rad * cls.EQ_RADIUS

    @classmethod
    def meters_to_radians(cls, m):
        return m / cls.EQ_RADIUS

    # Properties

    @property
    def latitude_radians(self):
        return self.geographic_latitude(math.atan2(self.z,
                                                   math.sqrt(self.x**2 + self.y**2)))
    @property
    def longitude_radians(self):
        return math.atan2(self.y, self.x)

    @property
    def latitude(self):
        return self.latitude_radians * 180.0 / math.pi

    @property
    def longitude(self):
        return self.longitude_radians * 180.0 / math.pi

    @property
    def antipode(self):
        return self.scale(-1.0)

    # Methods

    def __str__(self):
        return '%f, %f' % (self.longitude, self.latitude)

    def dot(self, b):
        return (self.x * b.x) + (self.y * b.y) + (self.z * b.z)

    def cross_length(self, b):
        x = (self.y * b.z) - (self.z * b.y)
        y = (self.z * b.x) - (self.x * b.z)
        z = (self.x * b.y) - (self.y * b.x)
        return math.sqrt(x**2 + y**2 + z**2)

    def scale(self, s):
        r = Coordinates(0,0)
        r.x = self.x * s
        r.y = self.y * s
        r.z = self.z * s
        return r

    def cross_normalize(self, b):
        x = (self.y * b.z) - (self.z * b.y)
        y = (self.z * b.x) - (self.x * b.z)
        z = (self.x * b.y) - (self.y * b.x)
        L = math.sqrt(x**2 + y**2 + z**2)
        r = Coordinates(0,0)
        r.x = x / L
        r.y = y / L
        r.z = z / L
        return r

    def distance_radians(self, b):
        return math.atan2(b.cross_length(self), b.dot(self))

    def distance(self, b):
        return self.radians_to_meters(self.distance_radians(b))

    def distance_to_line_segment(self, geo1, geo2):
        p2 = geo1.cross_normalize(geo2)
        ip = self.get_intersection(geo1, geo2, self, p2)

        d = geo1.distance(geo2)
        d1p = geo1.distance(ip)
        d2p = geo2.distance(ip)
        if d >= d1p and d >= d2p:
            return self.distance(ip)

        ip = ip.antipode
        d1p = geo1.distance(ip)
        d2p = geo2.distance(ip)
        if d >= d1p and d >= d2p:
            return self.distance(ip)

        return min(geo1.distance(self), geo2.distance(self))

    def distance_to_linestring(self, polyline):
        d = sys.maxint
        c = polyline
        for l1, l2 in [(c[i], c[i+1]) for i in xrange(len(c) - 1)]:
            dp = self.distance_to_line_segment(l1, l2)
            if dp < d:
                d = dp;
        return d

    def hit_test(self, polygon):
        counter = 0
        c = polygon
        n = len(p)
        for p1, p2 in [(c[i], c[i+1]) for i in xrange(len(c) - 1)]:
            if self.y > min(p1.y, p2.y):
                if self.y <= max(p1.y, p2.y):
                    if self.x <= max(p1.x, p2.x):
                        if p1.y != p2.y:
                            xinters = (self.y - p1.y) * (p2.x - p1.x) / \
                                (p2.y-p1.y) + p1.x
                            if p1.x == p2.x or self.x <= xinters:
                                counter += 1

        return counter % 2 != 0

    def distance_to_polygon(self, polygon):
        if self.hit_test(polygon):
            return 0.0
        return self.distance_to_linestring(polygon)

class Shape(object):
    def __init__(self, coordinates, name=None):
        self.name = name
        self._coordinates = [Coordinates(lon, lat) for lon, lat in coordinates]

    @property
    def coordinates(self):
        return self._coordinates

class Linestring(Shape):
    def __init__(self, coordinates, name=None):
        if len(coordinates) < 2:
            raise TypeError, "Need at least two points for a line."
        Shape.__init__(self, coordinates, name)

class Polygon(Shape):
    def __init__(self, coordinates, name=None):
        if len(coordinates) < 3:
            raise TypeError, "Need at least three points for a polygon."
        if coordinates[0] != coordinates[-1]:
            raise TypeError, "Polygon coordinates must start and end with the same one"
        Shape.__init__(self, coordinates, name)

class Point(Shape):
    def __init__(self, coordinates, name=None):
        if len(coordinates) != 1:
            raise TypeError, "A point is one set of coordinates"
        Shape.__init__(self, coordinates, name)

class ShapeFactory(object):
    def __init__(self):
        self.wkt_pattern = re.compile('^\s*(\w+)\s*\(+(.*?)\)+$')
        self.shape_classes = {"POINT" : Point,
                              "POLYGON" : Polygon,
                              "LINESTRING" : Linestring}

    def __call__(self, wkt, name=None):
        shape, coords = self.wkt_pattern.match(wkt).groups()
        cls = self.shape_classes[shape]
        return cls([[float(c) for c in b.strip().split(' ')] \
                    for b in coords.split(',')], name)

if __name__ == "__main__":
    jdata = urlopen(
        'http://staging.monotonous.org/device/distances/sound-areas').read()

    features = json.loads(jdata)['features']

    pprint(features)
    
    points = {}
    cpoints = {}

    for f in features:
        if f['geometry']['type'] == 'Point':
            coords = f['geometry']['coordinates']
            points[f['properties']['name']] = coords
            cpoints[f['properties']['name']] = Coordinates(*coords)

    plistlib.writePlist(points, 'points.plist')

    lines = {}
    clines = {}

    for f in features:
        if f['geometry']['type'] == 'LineString':
            coords = f['geometry']['coordinates']
            lines[f['properties']['name']] = coords
            clines[f['properties']['name']] = \
                [Coordinates(*c) for c in coords]

    plistlib.writePlist(lines, 'lines.plist')

    polys = {}
    cpolys = {}

    for f in features:
        if f['geometry']['type'] == 'Polygon':
            coords = f['geometry']['coordinates'][0]
            polys[f['properties']['name']] = coords
            cpolys[f['properties']['name']] = \
                [Coordinates(*c) for c in coords]

    plistlib.writePlist(polys, 'polygons.plist')

    for name, point in cpoints.items():
        for n, p in cpoints.items():
            print '%s -> %s: %f' % (name, n, point.distance(p))
        for n, l in clines.items():
           print '%s -> %s: %f' % (name, n, point.distance_to_linestring(l))
        for n, p in cpolys.items():
            print '%s -> %s: %f' % (name, n, point.distance_to_polygon(p))