geodesy.py

from math import fmod, modf, sqrt, degrees, radians, sin, cos, atan2


# Mean Earth radius defined by IUGG. (Unit: meters)
# ref. https://en.wikipedia.org/wiki/Earth_radius#Mean_radius
EARTH_RADIUS = 6371009.


def distance((lat1, lng1), (lat2, lng2)):
    """Get approximate geographical distance between 2 coordinates in meters.
    """
    # Based on Haversine Formula: https://en.wikipedia.org/wiki/Haversine_formula
    lat1, lng1, lat2, lng2 = map(float, [lat1, lng1, lat2, lng2])
    a = sin(radians(lat2 - lat1) / 2) ** 2 + \
        cos(radians(lat1)) * cos(radians(lat2)) * \
        sin(radians(lng2 - lng1) / 2) ** 2
    c = 2 * atan2(sqrt(a), sqrt(1 - a))
    d = EARTH_RADIUS * c

    return d


def search_area(center, radius=1000):
    """Get a square area, each edge length is (2 * radius) meters.
    """
    assert radius > 0
    lat, lng = map(float, center)

    # Get north and south edge
    dlat = degrees(float(radius) / EARTH_RADIUS)
    north = Degree(decimal=fmod(lat + dlat + 90, 90 * 2) - 90)
    south = Degree(decimal=fmod(lat - dlat + 90, 90 * 2) - 90)

    # Get east and west edge
    # Based on Aviation Formulary: http://williams.best.vwh.net/avform.htm#LL
    lat_rad = radians(lat)
    dist = float(radius) / EARTH_RADIUS
    dlng = degrees(atan2(sin(dist) * cos(lat_rad), cos(dist) - sin(lat_rad) ** 2))
    east = Degree(decimal=fmod(lng + dlng + 180, 180 * 2) - 180)
    west = Degree(decimal=fmod(lng - dlng + 180, 180 * 2) - 180)

    return north, east, south, west


class Degree(object):
    def __init__(self, deg=(0, 1), min=(0, 1), sec=(0, 1), degree_ref=None,
                 decimal=None, precision=1000):
        """Create a Degree object.

        If decimal is used, deg/min/sec and degree_ref will be ignored.
        """
        if decimal is not None:
            assert precision > 0
            remainder, deg = modf(abs(decimal))
            remainder, min = modf(remainder * 60)
            self.sign = -1 if decimal < 0 else 1
            self.degrees = (int(deg), 1)
            self.minutes = (int(min), 1)
            self.seconds = (int(round(remainder * 60 * precision)), precision)
        else:
            self.sign = -1 if degree_ref is not None and degree_ref[0].upper() in 'SW' else 1
            self.degrees = deg
            self.minutes = min
            self.seconds = sec

    @staticmethod
    def fraction_to_decimal(fraction):
        numerator, denominator = fraction
        if numerator == 0 or denominator == 0:
            return 0.
        else:
            return float(numerator) / denominator

    def __float__(self):
        d = self.fraction_to_decimal(self.degrees)
        m = self.fraction_to_decimal(self.minutes)
        s = self.fraction_to_decimal(self.seconds)
        return self.sign * (d + m / 60 + s / 3600)

    def __str__(self):
        d = self.fraction_to_decimal(self.degrees)
        m = self.fraction_to_decimal(self.minutes)
        s = self.fraction_to_decimal(self.seconds)
        return '%s\xa2X %s\' %s"' % (d, m, s)

    def __unicode__(self):
        return str(self).decode('mbcs')

    def __repr__(self):
        return str(self)


if __name__ == '__main__':
    from random import random, randrange

    for _ in xrange(10):
        lat1 = Degree(decimal=90 * random() * randrange(-1, 2, 2))
        lng1 = Degree(decimal=180 * random() * randrange(-1, 2, 2))
        lat2 = Degree(decimal=90 * random() * randrange(-1, 2, 2))
        lng2 = Degree(decimal=180 * random() * randrange(-1, 2, 2))

        print 'Distance between:'
        print '  %s' % str((lat1, lng1))
        print '  %s' % str((lat2, lng2))
        print 'Result:'
        print '  %s (meters)' % distance((lat1, lng1), (lat2, lng2))
        print ''

    for _ in xrange(10):
        lat1 = Degree(decimal=90 * random() * randrange(-1, 2, 2))
        lng1 = Degree(decimal=180 * random() * randrange(-1, 2, 2))

        print 'Search area:'
        print '  %s' % str((lat1, lng1))
        print 'Result:'
        north, east, south, west = search_area((lat1, lng1))
        print '  upper_right:', (north, east)
        print '  lower_left:', (south, west)
        print ''