[points] haversine distance in a different method

scripts/geodata/distance/haversine.py

@@ -0,0 +1,32 @@
+import math
+
+EARTH_RADIUS_KM = 6373
+
+
+def haversine_distance(lat1, lon1, lat2, lon2, radius=EARTH_RADIUS_KM):
+    """Calculate the Haversine distance between two lat/lon pairs, given by:
+    a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
+    c = 2 ⋅ atan2( √a, √(1−a) )
+    d = R ⋅ c
+
+    where R is the radius of the Earth (in kilometers). By default we use 6373 km,
+    a radius optimized for calculating distances at approximately 39 degrees from
+    the equator i.e. Washington, DC
+
+    :param lat1: first latitude
+    :param lon1: first longitude (use negative range for longitudes West of the Prime Meridian)
+    :param lat2: second latitude
+    :param lon2: second longitude (use negative range for longitudes West of the Prime Meridian)
+    :param radius: radius of the Earth in (miles|kilometers) depending on the desired units
+    """
+    lat1 = math.radians(lat1)
+    lat2 = math.radians(lat2)
+    lon1 = math.radians(lon1)
+    lon2 = math.radians(lon2)
+
+    dlon = lon2 - lon1
+    dlat = lat2 - lat1
+    a = (math.sin(dlat / 2.0)) ** 2 + math.cos(lat1) * math.cos(lat2) * (math.sin(dlon/2.0)) ** 2
+    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
+    d = radius * c
+    return d

scripts/geodata/points/index.py

@@ -9,7 +9,7 @@ from collections import defaultdict, OrderedDict
 
 from leveldb import LevelDB
 
-EARTH_RADIUS_KM = 6373  # in km
+from geodata.distance.haversine import haversine_distance
 
 
 class PointIndex(object):
@@ -103,34 +103,6 @@ class PointIndex(object):
         self.save_index()
         self.save_properties(os.path.join(self.save_dir, self.DEFAULT_PROPS_FILENAME))
 
-    def haversine_distance(self, lat1, lon1, lat2, lon2, radius=EARTH_RADIUS_KM):
-        """Calculate the Haversine distance between two lat/lon pairs, given by:
-        a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
-        c = 2 ⋅ atan2( √a, √(1−a) )
-        d = R ⋅ c
-
-        where R is the radius of the Earth (in kilometers). By default we use 6373 km,
-        a radius optimized for calculating distances at approximately 39 degrees from
-        the equator i.e. Washington, DC
-
-        :param lat1: first latitude
-        :param lon1: first longitude (use negative range for longitudes West of the Prime Meridian)
-        :param lat2: second latitude
-        :param lon2: second longitude (use negative range for longitudes West of the Prime Meridian)
-        :param radius: radius of the Earth in (miles|kilometers) depending on the desired units
-        """
-        lat1 = math.radians(lat1)
-        lat2 = math.radians(lat2)
-        lon1 = math.radians(lon1)
-        lon2 = math.radians(lon2)
-
-        dlon = lon2 - lon1
-        dlat = lat2 - lat1
-        a = (math.sin(dlat / 2.0)) ** 2 + math.cos(lat1) * math.cos(lat2) * (math.sin(dlon/2.0)) ** 2
-        c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
-        d = radius * c
-        return d
-
     def get_candidate_points(self, latitude, longitude):
         code = geohash.encode(latitude, longitude)[:self.precision]
         candidates = OrderedDict()
@@ -144,7 +116,7 @@ class PointIndex(object):
 
     def point_distances(self, latitude, longitude):
         candidates = self.get_candidate_points(latitude, longitude)
-        return [(i, lat, lon, self.haversine_distance(latitude, longitude, lat, lon)) for i, lat, lon in candidates]
+        return [(i, lat, lon, haversine_distance(latitude, longitude, lat, lon)) for i, lat, lon in candidates]
 
     def nearest_n_points(self, latitude, longitude, n=2):
         distances = self.point_distances(latitude, longitude)