From 36c1792d19ae48ea9f2b9d38ec3c8b05e94fdbb5 Mon Sep 17 00:00:00 2001 From: eightysixth <25541207+eightysixth@users.noreply.github.com> Date: Sat, 15 Feb 2020 23:32:56 -0700 Subject: [PATCH 1/6] implemented haversine --- geodesy/haversine_distance.py | 49 +++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 geodesy/haversine_distance.py diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py new file mode 100644 index 000000000000..20eb57488744 --- /dev/null +++ b/geodesy/haversine_distance.py @@ -0,0 +1,49 @@ +from math import asin, atan, cos, sin, sqrt, tan, pow, radians + + +def haversine_distance(lat1, lon1, lat2, lon2): + """ + Calculate great circle distance between two points in a sphere, + given longitudes and latitudes. + (https://en.wikipedia.org/wiki/Haversine_formula) + + Args: + lat1, lon1: latitude and longitude of coordinate 1 + lat2, lon2: latitude and longitude of coordinate 2 + returnAngle: Toggle to return distance or angle + + Returns: + geographical distance between two points in metres + + >>> int(haversine_distance(37.774856, -122.424227, 37.864742, -119.537521)) # From SF to Yosemite + 254352 + + """ + + # CONSTANTS per WGS84 https://en.wikipedia.org/wiki/World_Geodetic_System + AXIS_A = 6378137.0 + AXIS_B = 6356752.314245 + RADIUS = 6378137 + + # Equation parameters + # Equation https://en.wikipedia.org/wiki/Haversine_formula#Formulation + flattening = (AXIS_A - AXIS_B) / AXIS_A + phi_1 = atan((1 - flattening) * tan(radians(lat1))) + phi_2 = atan((1 - flattening) * tan(radians(lat2))) + lambda_1 = radians(lon1) + lambda_2 = radians(lon2) + + # Equation + sin_sq_phi = pow(sin((phi_2 - phi_1) / 2), 2) + sin_sq_lambda = pow(sin((lambda_2 - lambda_1) / 2), 2) + h_value = sqrt(sin_sq_phi + (cos(phi_1) * cos(phi_2) * sin_sq_lambda)) + + distance = 2 * RADIUS * asin(h_value) + + return distance + + +if __name__ == "__main__": + import doctest + + doctest.testmod() From 395a806db6f84401d3a6dfded940d33cf7b09889 Mon Sep 17 00:00:00 2001 From: eightysixth <25541207+eightysixth@users.noreply.github.com> Date: Sun, 16 Feb 2020 00:19:20 -0700 Subject: [PATCH 2/6] updated docstring Only calculate distance --- geodesy/haversine_distance.py | 1 - 1 file changed, 1 deletion(-) diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py index 20eb57488744..edfc8ea111e4 100644 --- a/geodesy/haversine_distance.py +++ b/geodesy/haversine_distance.py @@ -10,7 +10,6 @@ def haversine_distance(lat1, lon1, lat2, lon2): Args: lat1, lon1: latitude and longitude of coordinate 1 lat2, lon2: latitude and longitude of coordinate 2 - returnAngle: Toggle to return distance or angle Returns: geographical distance between two points in metres From e440ed2c2fe0d7b004e9a16ea7bd7a4b3ea2766a Mon Sep 17 00:00:00 2001 From: eightysixth <25541207+eightysixth@users.noreply.github.com> Date: Sun, 16 Feb 2020 00:25:47 -0700 Subject: [PATCH 3/6] added type hints --- geodesy/haversine_distance.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py index 20eb57488744..478166e5e965 100644 --- a/geodesy/haversine_distance.py +++ b/geodesy/haversine_distance.py @@ -1,7 +1,7 @@ from math import asin, atan, cos, sin, sqrt, tan, pow, radians -def haversine_distance(lat1, lon1, lat2, lon2): +def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float: """ Calculate great circle distance between two points in a sphere, given longitudes and latitudes. From 3dfd133a741ad5b011cff91cc6092571eebb4e4a Mon Sep 17 00:00:00 2001 From: eightysixth <25541207+eightysixth@users.noreply.github.com> Date: Sun, 16 Feb 2020 00:27:12 -0700 Subject: [PATCH 4/6] added type hints --- geodesy/haversine_distance.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py index edfc8ea111e4..bbf35004e5ac 100644 --- a/geodesy/haversine_distance.py +++ b/geodesy/haversine_distance.py @@ -1,7 +1,7 @@ from math import asin, atan, cos, sin, sqrt, tan, pow, radians -def haversine_distance(lat1, lon1, lat2, lon2): +def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float: """ Calculate great circle distance between two points in a sphere, given longitudes and latitudes. From 73e8401f834366ca95cb9de03ea0ca610cc5de08 Mon Sep 17 00:00:00 2001 From: eightysixth <25541207+eightysixth@users.noreply.github.com> Date: Sun, 16 Feb 2020 01:48:33 -0700 Subject: [PATCH 5/6] improved docstring and math usage --- geodesy/haversine_distance.py | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py index bbf35004e5ac..18330ee4f153 100644 --- a/geodesy/haversine_distance.py +++ b/geodesy/haversine_distance.py @@ -1,11 +1,20 @@ -from math import asin, atan, cos, sin, sqrt, tan, pow, radians +from math import asin, atan, cos, sin, sqrt, tan, radians def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float: """ Calculate great circle distance between two points in a sphere, - given longitudes and latitudes. - (https://en.wikipedia.org/wiki/Haversine_formula) + given longitudes and latitudes https://en.wikipedia.org/wiki/Haversine_formula + + We know that the globe is "sort of" spherical, so a path between two points + isn't exactly a straight line. We need to account for the Earth's curvature + when calculating distance from point A to B. This effect is negligible for + small distances but adds up as distance increases. The Haversine method treats + the earth as a sphere which allows us to "project" the two points A and B + onto the surface of that sphere and approximate the spherical distance between + them. Since the Earth is not a perfect sphere, other methods which model the + Earth's ellipsoidal nature are more accurate but a quick and modifiable + computation like Haversine can be handy for shorter range distances. Args: lat1, lon1: latitude and longitude of coordinate 1 @@ -14,12 +23,13 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl Returns: geographical distance between two points in metres - >>> int(haversine_distance(37.774856, -122.424227, 37.864742, -119.537521)) # From SF to Yosemite + >>> int(haversine_distance(37.774856, -122.424227, 37.864742, -119.537521)) # From SF to Yosemite in metres. 254352 """ # CONSTANTS per WGS84 https://en.wikipedia.org/wiki/World_Geodetic_System + # Distance in metres(m) AXIS_A = 6378137.0 AXIS_B = 6356752.314245 RADIUS = 6378137 @@ -33,8 +43,12 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl lambda_2 = radians(lon2) # Equation - sin_sq_phi = pow(sin((phi_2 - phi_1) / 2), 2) - sin_sq_lambda = pow(sin((lambda_2 - lambda_1) / 2), 2) + sin_sq_phi = sin((phi_2 - phi_1) / 2) + sin_sq_lambda = sin((lambda_2 - lambda_1) / 2) + # Square both values + sin_sq_phi *= sin_sq_phi + sin_sq_lambda *= sin_sq_lambda + h_value = sqrt(sin_sq_phi + (cos(phi_1) * cos(phi_2) * sin_sq_lambda)) distance = 2 * RADIUS * asin(h_value) From b5d09f3d44c404423ccc70ab504575dafb8efe9e Mon Sep 17 00:00:00 2001 From: Christian Clauss Date: Sun, 16 Feb 2020 10:53:46 +0100 Subject: [PATCH 6/6] f"{haversine_distance(*SAN_FRANCISCO, *YOSEMITE):0,.0f} meters" --- geodesy/haversine_distance.py | 34 ++++++++++++++-------------------- 1 file changed, 14 insertions(+), 20 deletions(-) diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py index 18330ee4f153..de8ac7f88302 100644 --- a/geodesy/haversine_distance.py +++ b/geodesy/haversine_distance.py @@ -1,4 +1,4 @@ -from math import asin, atan, cos, sin, sqrt, tan, radians +from math import asin, atan, cos, radians, sin, sqrt, tan def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float: @@ -6,34 +6,33 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl Calculate great circle distance between two points in a sphere, given longitudes and latitudes https://en.wikipedia.org/wiki/Haversine_formula - We know that the globe is "sort of" spherical, so a path between two points - isn't exactly a straight line. We need to account for the Earth's curvature - when calculating distance from point A to B. This effect is negligible for + We know that the globe is "sort of" spherical, so a path between two points + isn't exactly a straight line. We need to account for the Earth's curvature + when calculating distance from point A to B. This effect is negligible for small distances but adds up as distance increases. The Haversine method treats - the earth as a sphere which allows us to "project" the two points A and B + the earth as a sphere which allows us to "project" the two points A and B onto the surface of that sphere and approximate the spherical distance between them. Since the Earth is not a perfect sphere, other methods which model the Earth's ellipsoidal nature are more accurate but a quick and modifiable - computation like Haversine can be handy for shorter range distances. + computation like Haversine can be handy for shorter range distances. Args: lat1, lon1: latitude and longitude of coordinate 1 lat2, lon2: latitude and longitude of coordinate 2 - - Returns: + Returns: geographical distance between two points in metres - - >>> int(haversine_distance(37.774856, -122.424227, 37.864742, -119.537521)) # From SF to Yosemite in metres. - 254352 - + >>> from collections import namedtuple + >>> point_2d = namedtuple("point_2d", "lat lon") + >>> SAN_FRANCISCO = point_2d(37.774856, -122.424227) + >>> YOSEMITE = point_2d(37.864742, -119.537521) + >>> f"{haversine_distance(*SAN_FRANCISCO, *YOSEMITE):0,.0f} meters" + '254,352 meters' """ - # CONSTANTS per WGS84 https://en.wikipedia.org/wiki/World_Geodetic_System # Distance in metres(m) AXIS_A = 6378137.0 AXIS_B = 6356752.314245 RADIUS = 6378137 - # Equation parameters # Equation https://en.wikipedia.org/wiki/Haversine_formula#Formulation flattening = (AXIS_A - AXIS_B) / AXIS_A @@ -41,19 +40,14 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl phi_2 = atan((1 - flattening) * tan(radians(lat2))) lambda_1 = radians(lon1) lambda_2 = radians(lon2) - # Equation sin_sq_phi = sin((phi_2 - phi_1) / 2) sin_sq_lambda = sin((lambda_2 - lambda_1) / 2) # Square both values sin_sq_phi *= sin_sq_phi sin_sq_lambda *= sin_sq_lambda - h_value = sqrt(sin_sq_phi + (cos(phi_1) * cos(phi_2) * sin_sq_lambda)) - - distance = 2 * RADIUS * asin(h_value) - - return distance + return 2 * RADIUS * asin(h_value) if __name__ == "__main__":