Logo Search packages:      
Sourcecode: d-rats version File versions  Download package

distance.py

from math import *
import util

# Average great-circle radius in kilometers, from Wikipedia.
# Using a sphere with this radius results in an error of up to about 0.5%.
EARTH_RADIUS = 6372.795

# From http://www.movable-type.co.uk/scripts/LatLongVincenty.html:
#   The most accurate and widely used globally-applicable model for the earth
#   ellipsoid is WGS-84, used in this script. Other ellipsoids offering a
#   better fit to the local geoid include Airy (1830) in the UK, International
#   1924 in much of Europe, Clarke (1880) in Africa, and GRS-67 in South
#   America. America (NAD83) and Australia (GDA) use GRS-80, functionally
#   equivalent to the WGS-84 ellipsoid.
#
#             model             major (km)   minor (km)     flattening
ELLIPSOIDS = {'WGS-84':        (6378.137,    6356.7523142,  1 / 298.257223563),
              'GRS-80':        (6378.137,    6356.7523141,  1 / 298.257222101),
              'Airy (1830)':   (6377.563396, 6356.256909,   1 / 299.3249646),
              'Intl 1924':     (6378.388,    6356.911946,   1 / 297.0),
              'Clarke (1880)': (6378.249145, 6356.51486955, 1 / 293.465),
              'GRS-67':        (6378.1600,   6356.774719,   1 / 298.25),
              }

def arc_degrees(arcminutes=0, arcseconds=0):
    """Calculate the decimal equivalent of the sum of ``arcminutes`` and
    ``arcseconds`` in degrees."""
    if arcminutes is None:
        arcminutes = 0
    if arcseconds is None:
        arcseconds = 0
    arcmin = float(arcminutes)
    arcsec = float(arcseconds)
    return arcmin * 1 / 60. + arcsec * 1 / 3600.

def kilometers(miles=0, feet=0, nautical=0):
    d = 0
    if feet:
        miles += feet / ft(1.0)
    if nautical:
        d += nautical / nm(1.0)
    d += miles * 1.609344
    return d

def miles(kilometers=0, feet=0, nautical=0):
    d = 0
    if nautical:
        kilometers += nautical / nm(1.0)
    if feet:
        d += feet / ft(1.0)
    d += kilometers * 0.621371192
    return d

def feet(miles=0, kilometers=0, nautical=0):
    d = 0
    if nautical:
        kilometers += nautical / nm(1.0)
    if kilometers:
        miles += mi(kilometers)
    d += miles * 5280
    return d

def nautical(kilometers=0, miles=0, feet=0):
    d = 0
    if feet:
        miles += feet / ft(1.0)
    if miles:
        kilometers += km(miles)
    d += kilometers / 1.852
    return d

km = kilometers
mi = miles
ft = feet
nm = nautical


00078 class Distance(object):
00079     def __init__(self, kilometers=0, miles=0, feet=0, nautical=0):
        """Initialize a Distance whose length is the sum of all the units
        measured in the constructor (kilometers, miles, feet, nautical miles).
        """
        kilometers += km(miles=miles, feet=feet, nautical=nautical)
        self._kilometers = kilometers

    @property
    def kilometers(self):
        return self._kilometers
    
    @property
    def miles(self):
        return miles(self.kilometers)

    @property
    def feet(self):
        return feet(self.miles)

    @property
    def nautical(self):
        return nautical(self.kilometers)

    # Sadly, just aliasing the above properties with their abbreviations does
    # not work when they are subclassed. The easiest way I could find to
    # make this work without using a metaclass was to write more full-fledged
    # definitions...

    @property
    def mi(self):
        return self.miles
    
    @property
    def km(self):
        return self.kilometers
    
    @property
    def ft(self):
        return self.feet

    @property
    def nm(self):
        return self.nautical

00123     def __add__(self, other):
        """Return a new Distance of length ``self`` + ``other``."""
        if isinstance(other, Distance):
            return Distance(self.kilometers + other.kilometers)
        else:
            raise TypeError("Distance must be added with Distance.")

00130     def __sub__(self, other):
        """Return a new Distance of length ``self`` - ``other``."""
        if isinstance(other, Distance):
            return Distance(self.kilometers - other.kilometers)
        else:
            raise TypeError("Distance must be subtracted from Distance.")
        
00137     def __mul__(self, other):
        """Return a new Distance ``other`` times the length of ``self``,
        ``other`` is a number (int, float, long, or Decimal).
        """
        if isinstance(other, (int, float, long, Decimal)):
            other = float(other)
            return Distance(self.kilometers * other)
        else:
            raise TypeError("Distance must be multiplied by number.")

00147     def __div__(self, other):
        """If ``other`` is a number (int, float, long, or Decimal), return
        a new Distance of length ``self`` / ``other``.
        
        If ``other`` is a Distance, return the fraction given by
        ``self`` / ``other``.
        """
        if isinstance(other, Distance):
            return float(self.kilometers) / other.kilometers
        elif isinstance(other, (int, float, long, Decimal)):
            other = float(other)
            return Distance(self.kilometers / other)
        else:
            raise TypeError("Distance must be divided by Distance or number.")

00162     def __nonzero__(self):
        """Return whether or not this Distance is 0 units in length."""
        return bool(self.kilometers)


00167 class GeodesicDistance(Distance):
00168     def __init__(self, a, b):
        """Initialize a Distance whose length is the distance between the two
        geodesic points ``a`` and ``b``, using the ``calculate`` method to
        determine this distance.
        """
        if isinstance(a, basestring):
            a = util.parse_geo(a)
        if isinstance(b, basestring):
            b = util.parse_geo(b)
        
        self.a = a
        self.b = b
        
        if a and b:
            self.calculate()

00184     def calculate(self):
        """Calculate and set the distance between ``self.a`` and ``self.b``,
        which should be two geodesic points. Since there are multiple formulas
        to calculate this, implementation is left up to the subclass.
        """
        raise NotImplementedError

    @property
    def kilometers(self):
        raise NotImplementedError

00195 class GreatCircleDistance(GeodesicDistance):
    """Use spherical geometry to calculate the surface distance between two
    geodesic points. This formula can be written many different ways, including
    just the use of the spherical law of cosines or the haversine formula.
    
    The class member ``RADIUS`` indicates which radius of the earth to use,
    in kilometers. The default is to use the module constant ``EARTH_RADIUS``,
    which uses the average great-circle radius.
    """
    
    RADIUS = EARTH_RADIUS
    
00207     def calculate(self):
        lat1, lng1 = map(radians, self.a)
        lat2, lng2 = map(radians, self.b)
        
        sin_lat1, cos_lat1 = sin(lat1), cos(lat1)
        sin_lat2, cos_lat2 = sin(lat2), cos(lat2)
        
        delta_lng = lng2 - lng1
        cos_delta_lng, sin_delta_lng = cos(delta_lng), sin(delta_lng)
        
        central_angle = acos(sin_lat1 * sin_lat2 +
                             cos_lat1 * cos_lat2 * cos_delta_lng)
        
        # From http://en.wikipedia.org/wiki/Great_circle_distance:
        #   Historically, the use of this formula was simplified by the
        #   availability of tables for the haversine function. Although this
        #   formula is accurate for most distances, it too suffers from
        #   rounding errors for the special (and somewhat unusual) case of
        #   antipodal points (on opposite ends of the sphere). A more
        #   complicated formula that is accurate for all distances is: (below)
        
        d = atan2(sqrt((cos_lat2 * sin_delta_lng) ** 2 +
                       (cos_lat1 * sin_lat2 -
                        sin_lat1 * cos_lat2 * cos_delta_lng) ** 2),
                  sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_delta_lng)
        
        self.radians = d
    
    @property
    def kilometers(self):
        return self.RADIUS * self.radians
    

00240 class VincentyDistance(GeodesicDistance):
    """Calculate the geodesic distance between two points using the formula
    devised by Thaddeus Vincenty, with an accurate ellipsoidal model of the
    earth.
    
    The class attribute ``ELLIPSOID`` indicates which ellipsoidal model of the
    earth to use. If it is a string, it is looked up in the ELLIPSOIDS
    dictionary to obtain the major and minor semiaxes and the flattening.
    Otherwise, it should be a tuple with those values. The most globally
    accurate model is WGS-84. See the comments above the ELLIPSOIDS dictionary
    for more information.
    """

    ELLIPSOID = ELLIPSOIDS['WGS-84']
    
00255     def calculate(self):
        lat1, lng1 = map(radians, self.a)
        lat2, lng2 = map(radians, self.b)
        
        if isinstance(self.ELLIPSOID, basestring):
            major, minor, f = ELLIPSOIDS[self.ELLIPSOID]
        else:
            major, minor, f = self.ELLIPSOID
        
        delta_lng = lng2 - lng1
        
        reduced_lat1 = atan((1 - f) * tan(lat1))
        reduced_lat2 = atan((1 - f) * tan(lat2))
        
        sin_reduced1, cos_reduced1 = sin(reduced_lat1), cos(reduced_lat1)
        sin_reduced2, cos_reduced2 = sin(reduced_lat2), cos(reduced_lat2)
        
        lambda_lng = delta_lng
        lambda_prime = 2 * pi
        
        iter_limit = 20
        
        while abs(lambda_lng - lambda_prime) > 10e-12 and iter_limit > 0:
            sin_lambda_lng, cos_lambda_lng = sin(lambda_lng), cos(lambda_lng)
            
            sin_sigma = sqrt((cos_reduced2 * sin_lambda_lng) ** 2 +
                             (cos_reduced1 * sin_reduced2 - sin_reduced1 *
                              cos_reduced2 * cos_lambda_lng) ** 2)
            
            if sin_sigma == 0:
                # Coincident points
                self._kilometers = self.initial_bearing = self.final_bearing = 0
                return
            
            cos_sigma = (sin_reduced1 * sin_reduced2 +
                         cos_reduced1 * cos_reduced2 * cos_lambda_lng)
            
            sigma = atan2(sin_sigma, cos_sigma)
            
            sin_alpha = cos_reduced1 * cos_reduced2 * sin_lambda_lng / sin_sigma
            cos_sq_alpha = 1 - sin_alpha ** 2
            
            if cos_sq_alpha != 0:
                cos2_sigma_m = cos_sigma - 2 * (sin_reduced1 * sin_reduced2 /
                                                cos_sq_alpha)
            else:
                cos2_sigma_m = 0.0 # Equatorial line
            
            C = f / 16. * cos_sq_alpha * (4 + f * (4 - 3 * cos_sq_alpha))
            
            lambda_prime = lambda_lng
            lambda_lng = (delta_lng + (1 - C) * f * sin_alpha *
                          (sigma + C * sin_sigma *
                           (cos2_sigma_m + C * cos_sigma * 
                            (-1 + 2 * cos2_sigma_m ** 2))))
            iter_limit -= 1
            
        if iter_limit == 0:
            raise ValueError("Vincenty formula failed to converge!")
        
        u_sq = cos_sq_alpha * (major ** 2 - minor ** 2) / minor ** 2
        
        A = 1 + u_sq / 16384. * (4096 + u_sq * (-768 + u_sq *
                                                (320 - 175 * u_sq)))
        
        B = u_sq / 1024. * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq)))
        
        delta_sigma = (B * sin_sigma *
                       (cos2_sigma_m + B / 4. *
                        (cos_sigma * (-1 + 2 * cos2_sigma_m ** 2) -
                         B / 6. * cos2_sigma_m * (-3 + 4 * sin_sigma ** 2) *
                         (-3 + 4 * cos2_sigma_m ** 2))))
        
        s = minor * A * (sigma - delta_sigma)
        
        sin_lambda, cos_lambda = sin(lambda_lng), cos(lambda_lng)
        
        alpha_1 = atan2(cos_reduced2 * sin_lambda,
                        cos_reduced1 * sin_reduced2 -
                        sin_reduced1 * cos_reduced2 * cos_lambda)
        
        alpha_2 = atan2(cos_reduced1 * sin_lambda,
                        cos_reduced1 * sin_reduced2 * cos_lambda -
                        sin_reduced1 * cos_reduced2)
        
        self._kilometers = s
        self.initial_bearing = (360 + degrees(alpha_1)) % 360
        self.final_bearing = (360 + degrees(alpha_2)) % 360

    @property
    def kilometers(self):
        return self._kilometers

    @property
    def forward_azimuth(self):
        return self.initial_bearing


# Set the default distance formula to the most generally accurate.
distance = VincentyDistance


def destination(start, bearing, distance, radius=EARTH_RADIUS):
    lat1, lng1 = map(radians, start)
    bearing = radians(bearing)
    
    if isinstance(distance, Distance):
        distance = distance.kilometers
        
    d_div_r = float(distance) / radius
    
    lat2 = asin(sin(lat1) * cos(d_div_r) +
                cos(lat1) * sin(d_div_r) * cos(bearing))
    
    lng2 = lng1 + atan2(sin(bearing) * sin(d_div_r) * cos(lat1),
                        cos(d_div_r) - sin(lat1) * sin(lat2))
    
    return (degrees(lat2), degrees(lng2))


def vincenty_destination(start, bearing, distance,
                         ellipsoid=ELLIPSOIDS['WGS-84']):
    lat1, lng1 = map(radians, start)
    bearing = radians(bearing)
    
    if isinstance(distance, Distance):
        distance = distance.kilometers
    
    if isinstance(ellipsoid, basestring):
        ellipsoid = ELLIPSOIDS[ellipsoid]
    
    major, minor, f = ellipsoid
    
    tan_reduced1 = (1 - f) * tan(lat1)
    cos_reduced1 = 1 / sqrt(1 + tan_reduced1 ** 2)
    sin_reduced1 = tan_reduced1 * cos_reduced1
    sin_bearing, cos_bearing = sin(bearing), cos(bearing)
    sigma1 = atan2(tan_reduced1, cos_bearing)
    sin_alpha = cos_reduced1 * sin_bearing
    cos_sq_alpha = 1 - sin_alpha ** 2
    u_sq = cos_sq_alpha * (major ** 2 - minor ** 2) / minor ** 2
    
    A = 1 + u_sq / 16384. * (4096 + u_sq * (-768 + u_sq * (320 - 175 * u_sq)))
    B = u_sq / 1024. * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq)))

    sigma = distance / (minor * A)
    sigma_prime = 2 * pi
    
    while abs(sigma - sigma_prime) > 10e-12:
        cos2_sigma_m = cos(2 * sigma1 + sigma)
        sin_sigma, cos_sigma = sin(sigma), cos(sigma)
        delta_sigma = B * sin_sigma * (cos2_sigma_m + B / 4. *
                                       (cos_sigma * (-1 + 2 * cos2_sigma_m) -
                                        B / 6. * cos2_sigma_m *
                                        (-3 + 4 * sin_sigma ** 2) *
                                        (-3 + 4 * cos2_sigma_m ** 2)))
        sigma_prime = sigma
        sigma = distance / (minor * A) + delta_sigma
    
    sin_sigma, cos_sigma = sin(sigma), cos(sigma)
    
    lat2 = atan2(sin_reduced1 * cos_sigma +
                 cos_reduced1 * sin_sigma * cos_bearing,
                 (1 - f) * sqrt(sin_alpha ** 2 +
                                (sin_reduced1 * sin_sigma -
                                 cos_reduced1 * cos_sigma * cos_bearing) ** 2))
    
    lambda_lng = atan2(sin_sigma * sin_bearing,
                       cos_reduced1 * cos_sigma -
                       sin_reduced1 * sin_sigma * cos_bearing)
    
    C = f / 16. * cos_sq_alpha * (4 + f * (4 - 3 * cos_sq_alpha))
    
    delta_lng = (lambda_lng - (1 - C) * f * sin_alpha *
                 (sigma + C * sin_sigma *
                  (cos2_sigma_m + C * cos_sigma *
                   (-1 + 2 * cos2_sigma_m ** 2))))
    
    final_bearing = atan2(sin_alpha,
                          cos_reduced1 * cos_sigma * cos_bearing -
                          sin_reduced1 * sin_sigma)
    
    lng2 = lng1 + delta_lng
    
    return (degrees(lat2), degrees(lng2))

Generated by  Doxygen 1.6.0   Back to index