""" geodetic_dist.py V.0.1 by leonardo maffi. To compute distances on Earth given two point latitudes and longitudes. sperical_dist, approx_ellipsoid_dist and ellipsoid_dist give incresignly precise results, but note that this code isn't much tested. For most purposes approx_ellipsoid_dist is enough. Original sources: http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393241 http://www.codeguru.com/Cpp/Cpp/algorithms/article.php/c5115/ """ from math import sin, cos, sqrt, atan2, asin DE2RA = 0.01745329252 # Degrees to radians ERAD = 6378.137 FLATTENING = 1.000000 / 298.257223563 # Earth flattening (WGS84) EPS = 0.000000000005 def sperical_dist((lat1, lon1), (lat2, lon2)): """sperical_dist((lat1, lon1), (lat2, lon2)): returns the distance in kilometers between two points on the exactly sperical Earth. Input: lat1, lon1, lat2, lon2: latitude and longitude (in degrees) of two points on Earth. Output: distance in kilometers "crow fly" between the two points. If you want more precision use approx_geodesic_len. If you want even more precision use ellipsoid_distance.""" # Modified from Kevin Ryan code: # http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393241 lon1 *= DE2RA lon2 *= DE2RA lat1 *= DE2RA lat2 *= DE2RA dlon = lon2 - lon1 dlat = lat2 - lat1 a = (sin(dlat / 2.0))**2 + cos(lat1) * cos(lat2) * (sin(dlon / 2.0))**2 return 6366.56 * 2 * asin( min(1, sqrt(a)) ) def approx_ellipsoid_dist((lat1, lon1), (lat2, lon2)): """approx_ellipsoid_dist((lat1, lon1), (lat2, lon2)): Input: lat1, lon1, lat2, lon2: latitude and longitude (in degrees) of two points on Earth. Output: distance in kilometers "crow fly" between the two points. If you want less precision use sperical_distance. If you want more precision use ellipsoid_distance.""" # http://www.codeguru.com/Cpp/Cpp/algorithms/article.php/c5115/ # By Andy McGovern, translated to Python if lon1 == lon2 and lat1 == lat2: return 0.0 lat1 = DE2RA * lat1 lon1 = -DE2RA * lon1 lat2 = DE2RA * lat2 lon2 = -DE2RA * lon2 F = (lat1 + lat2) / 2.0 G = (lat1 - lat2) / 2.0 L = (lon1 - lon2) / 2.0 sing = sin(G) cosl = cos(L) cosf = cos(F) sinl = sin(L) sinf = sin(F) cosg = cos(G) S = sing*sing*cosl*cosl + cosf*cosf*sinl*sinl C = cosg*cosg*cosl*cosl + sinf*sinf*sinl*sinl W = atan2(sqrt(S), sqrt(C)) R = sqrt(S*C) / W H1 = (3 * R - 1.0) / (2.0 * C) H2 = (3 * R + 1.0) / (2.0 * S) D = 2 * W * ERAD return D * (1 + FLATTENING * H1 * sinf*sinf*cosg*cosg - FLATTENING*H2*cosf*cosf*sing*sing) def ellipsoid_dist((lat1, lon1), (lat2, lon2)): """ellipsoid_dist((lat1, lon1), (lat2, lon2)): Input: lat1, lon1, lat2, lon2: latitude and longitude (in degrees) of two points on Earth. Output: distance in kilometers "crow fly" between the two points. If you want less precision use approx_geodesic_len. If you want even less precision use sperical_distance.""" # http://www.codeguru.com/Cpp/Cpp/algorithms/article.php/c5115/ # By Andy McGovern, translated to Python distance = 0.0 r = 1.0 - FLATTENING distance = 0.0 if lon1 == lon2 and lat1 == lat2: return 0.0 lon1 *= DE2RA lon2 *= DE2RA lat1 *= DE2RA lat2 *= DE2RA cosy1 = cos(lat1) cosy2 = cos(lat2) if cosy1 == 0.0: cosy1 = 0.0000000001 if cosy2 == 0.0: cosy2 = 0.0000000001 tu1 = r * sin(lat1) / cosy1 tu2 = r * sin(lat2) / cosy2 cu1 = 1.0 / sqrt(tu1 * tu1 + 1.0) su1 = cu1 * tu1 cu2 = 1.0 / sqrt(tu2 * tu2 + 1.0) x = lon2 - lon1 distance = cu1 * cu2 baz = distance * tu2 faz = baz * tu1 while True: sx = sin(x) cx = cos(x) tu1 = cu2 * sx tu2 = baz - su1 * cu2 * cx sy = sqrt(tu1 * tu1 + tu2 * tu2) cy = distance * cx + faz y = atan2(sy, cy) sa = distance * sx / sy c2a = -sa * sa + 1.0 cz = faz + faz if c2a > 0.0: cz = -cz / c2a + cy e = cz * cz * 2. - 1.0 c = ((-3.0 * c2a + 4.0) * FLATTENING + 4.0) * c2a * FLATTENING / 16.0 d = x x = ((e * cy * c + cz) * sy * c + y) * sa x = (1.0 - c) * x * FLATTENING + lon2 - lon1 if abs(d - x) <= EPS: break x = sqrt((1.0 / r / r - 1.0) * c2a + 1.0) + 1.0 x = (x - 2.0) / x c = 1.0 - x c = (x * x / 4.0 + 1.0) / c d = (0.375 * x * x - 1.0) * x x = e * cy distance = 1.0 - e - e distance = ((((sy * sy * 4.0 - 3.0) * distance * cz * d / 6.0 - x) * d / 4.0 + cz) * sy * d + y) * c * ERAD * r return distance if __name__ == '__main__': testdata = ( [(34.086159000000002, -118.375984), (34.086159000000002, -118.375984), 0.000000, 0.000000, 0.000000], [(34.086159000000002, -118.375984), (32.779541000000002, -117.146344), 184.607705, 184.597416, 184.597132], [(34.086159000000002, -118.375984), (37.304051000000001, -121.87273399999999), 476.811534, 477.035894, 477.034999], [(34.086159000000002, -118.375984), (37.759881, -122.437392), 547.780078, 548.069548, 548.068515], [(34.086159000000002, -118.375984), (33.804133, -118.158028), 37.225313, 37.210025, 37.209926], [(34.086159000000002, -118.375984), (36.781548999999998, -119.792113), 325.781374, 325.513752, 325.512094], [(34.086159000000002, -118.375984), (38.555605, -121.468926), 568.518904, 568.331239, 568.328792], [(34.086159000000002, -118.375984), (37.795226999999997, -122.228111), 538.363056, 538.575883, 538.574713], [(34.086159000000002, -118.375984), (33.740716999999997, -117.88140799999999), 59.610044, 59.666414, 59.666396], [(34.086159000000002, -118.375984), (33.836165000000001, -117.889769), 52.722451, 52.805946, 52.805954], [(34.086159000000002, -118.375984), (33.948065, -117.39612700000001), 91.541596, 91.792418, 91.792417], [(34.086159000000002, -118.375984), (35.357275999999999, -119.031661), 153.412194, 153.268196, 153.267458], [(34.086159000000002, -118.375984), (37.975622999999999, -121.30086799999999), 505.759483, 505.644757, 505.642844], [(34.086159000000002, -118.375984), (37.542943000000001, -121.982786), 503.059631, 503.257601, 503.256542], [(34.086159000000002, -118.375984), (34.170938999999997, -118.25008099999999), 14.928440, 14.944005, 14.944001] ) for po1, po2, d1, d2, d3 in testdata: assert abs(sperical_dist(po1, po2) - d1) < 1e-6 assert abs(approx_ellipsoid_dist(po1, po2) - d2) < 1e-6 assert abs(ellipsoid_dist(po1, po2) - d3) < 1e-6