From 620f0594aa77337f678292466e8d9b073dea8624 Mon Sep 17 00:00:00 2001 From: Al Date: Wed, 27 Apr 2016 17:27:30 -0400 Subject: [PATCH] [points] haversine distance in a different method --- scripts/geodata/distance/haversine.py | 32 +++++++++++++++++++++++++++ scripts/geodata/points/index.py | 32 ++------------------------- 2 files changed, 34 insertions(+), 30 deletions(-) create mode 100644 scripts/geodata/distance/haversine.py diff --git a/scripts/geodata/distance/haversine.py b/scripts/geodata/distance/haversine.py new file mode 100644 index 00000000..576aaecf --- /dev/null +++ b/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 diff --git a/scripts/geodata/points/index.py b/scripts/geodata/points/index.py index 5f606eb7..0100c95e 100644 --- a/scripts/geodata/points/index.py +++ b/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)