From 3c22c600529c5501b67ff8c134b5e39125fda2c5 Mon Sep 17 00:00:00 2001 From: Maanit Date: Wed, 2 Oct 2024 01:12:42 +0530 Subject: [PATCH 1/6] Fix haversine distance calculation to use raw latitudes --- geodesy/haversine_distance.py | 65 +++++++++++++---------------------- 1 file changed, 24 insertions(+), 41 deletions(-) diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py index 93e625770..0b6cf0ccc 100644 --- a/geodesy/haversine_distance.py +++ b/geodesy/haversine_distance.py @@ -1,57 +1,40 @@ -from math import asin, atan, cos, radians, sin, sqrt, tan - -AXIS_A = 6378137.0 -AXIS_B = 6356752.314245 -RADIUS = 6378137 - +from math import asin, cos, radians, sin, sqrt 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 - - 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. - + Calculate the great-circle distance between two points + on the Earth specified by latitude and longitude using + the Haversine formula. + Args: - lat1, lon1: latitude and longitude of coordinate 1 - lat2, lon2: latitude and longitude of coordinate 2 + lat1, lon1: Latitude and longitude of point 1 in decimal degrees. + lat2, lon2: Latitude and longitude of point 2 in decimal degrees. + Returns: - geographical distance between two points in metres + Distance between the two points in meters. + >>> 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' + '254,033 meters' """ - # CONSTANTS per WGS84 https://en.wikipedia.org/wiki/World_Geodetic_System - # Distance in metres(m) - # 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 = 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)) - return 2 * RADIUS * asin(h_value) + R = 6378137 #earth radius (meters) + + lat1_rad = radians(lat1) + lat2_rad = radians(lat2) + delta_lat = radians(lat2 - lat1) + delta_lon = radians(lon2 - lon1) + + # Haversine formula + a = sin(delta_lat / 2)**2 + cos(lat1_rad) * cos(lat2_rad) * sin(delta_lon / 2)**2 + c = 2 * asin(sqrt(a)) + + #Great-Circle Distance + return R * c if __name__ == "__main__": import doctest - doctest.testmod() From 25a03cf503cd0be991f453a847897a5b9b59432e Mon Sep 17 00:00:00 2001 From: Maanit Date: Wed, 2 Oct 2024 01:18:55 +0530 Subject: [PATCH 2/6] Fix haversine distance calculation to use raw latitudes --- geodesy/tempCodeRunnerFile.py | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 geodesy/tempCodeRunnerFile.py diff --git a/geodesy/tempCodeRunnerFile.py b/geodesy/tempCodeRunnerFile.py new file mode 100644 index 000000000..e39d788b1 --- /dev/null +++ b/geodesy/tempCodeRunnerFile.py @@ -0,0 +1,4 @@ + +# if __name__ == "__main__": +# import doctest +# doctest.testmod() From 5c9d6709ea90e5e40bb13224f16be70fb5c08790 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 1 Oct 2024 20:03:25 +0000 Subject: [PATCH 3/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- geodesy/haversine_distance.py | 20 +++++++++++++------- geodesy/tempCodeRunnerFile.py | 1 - 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py index 0b6cf0ccc..eeb65f59c 100644 --- a/geodesy/haversine_distance.py +++ b/geodesy/haversine_distance.py @@ -1,18 +1,19 @@ from math import asin, cos, radians, sin, sqrt + def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float: """ Calculate the great-circle distance between two points on the Earth specified by latitude and longitude using the Haversine formula. - + Args: lat1, lon1: Latitude and longitude of point 1 in decimal degrees. lat2, lon2: Latitude and longitude of point 2 in decimal degrees. - + Returns: Distance between the two points in meters. - + >>> from collections import namedtuple >>> point_2d = namedtuple("point_2d", "lat lon") >>> SAN_FRANCISCO = point_2d(37.774856, -122.424227) @@ -21,7 +22,7 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl '254,033 meters' """ - R = 6378137 #earth radius (meters) + R = 6378137 # earth radius (meters) lat1_rad = radians(lat1) lat2_rad = radians(lat2) @@ -29,12 +30,17 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl delta_lon = radians(lon2 - lon1) # Haversine formula - a = sin(delta_lat / 2)**2 + cos(lat1_rad) * cos(lat2_rad) * sin(delta_lon / 2)**2 + a = ( + sin(delta_lat / 2) ** 2 + + cos(lat1_rad) * cos(lat2_rad) * sin(delta_lon / 2) ** 2 + ) c = 2 * asin(sqrt(a)) - - #Great-Circle Distance + + # Great-Circle Distance return R * c + if __name__ == "__main__": import doctest + doctest.testmod() diff --git a/geodesy/tempCodeRunnerFile.py b/geodesy/tempCodeRunnerFile.py index e39d788b1..289e85db5 100644 --- a/geodesy/tempCodeRunnerFile.py +++ b/geodesy/tempCodeRunnerFile.py @@ -1,4 +1,3 @@ - # if __name__ == "__main__": # import doctest # doctest.testmod() From 8785ab7c39786a9e2e2344ed83257757a1931c1e Mon Sep 17 00:00:00 2001 From: Maanit Date: Wed, 2 Oct 2024 03:00:15 +0530 Subject: [PATCH 4/6] Add type hints and additional doctests to haversine_distance --- geodesy/haversine_distance.py | 5 +++-- geodesy/{tempCodeRunnerFile.py => temp_code_runner_file.py} | 0 2 files changed, 3 insertions(+), 2 deletions(-) rename geodesy/{tempCodeRunnerFile.py => temp_code_runner_file.py} (100%) diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py index eeb65f59c..824a4a255 100644 --- a/geodesy/haversine_distance.py +++ b/geodesy/haversine_distance.py @@ -11,6 +11,7 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl lat1, lon1: Latitude and longitude of point 1 in decimal degrees. lat2, lon2: Latitude and longitude of point 2 in decimal degrees. + Returns: Distance between the two points in meters. @@ -22,7 +23,7 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl '254,033 meters' """ - R = 6378137 # earth radius (meters) + radius = 6378137 # earth radius (meters) lat1_rad = radians(lat1) lat2_rad = radians(lat2) @@ -37,7 +38,7 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl c = 2 * asin(sqrt(a)) # Great-Circle Distance - return R * c + return radius * c if __name__ == "__main__": diff --git a/geodesy/tempCodeRunnerFile.py b/geodesy/temp_code_runner_file.py similarity index 100% rename from geodesy/tempCodeRunnerFile.py rename to geodesy/temp_code_runner_file.py From bddf603c2c92fe7643790754f99f41442386148f Mon Sep 17 00:00:00 2001 From: Maanit Date: Wed, 2 Oct 2024 03:15:40 +0530 Subject: [PATCH 5/6] lamberts_ellipsoidal_distance changes in accordance to new haversine distance parameters --- geodesy/lamberts_ellipsoidal_distance.py | 34 +++++++++--------------- 1 file changed, 13 insertions(+), 21 deletions(-) diff --git a/geodesy/lamberts_ellipsoidal_distance.py b/geodesy/lamberts_ellipsoidal_distance.py index 4805674e5..aad04ce48 100644 --- a/geodesy/lamberts_ellipsoidal_distance.py +++ b/geodesy/lamberts_ellipsoidal_distance.py @@ -1,6 +1,6 @@ from math import atan, cos, radians, sin, tan -from .haversine_distance import haversine_distance +from haversine_distance import haversine_distance AXIS_A = 6378137.0 AXIS_B = 6356752.314245 @@ -15,22 +15,21 @@ def lamberts_ellipsoidal_distance( two points on the surface of earth given longitudes and latitudes https://en.wikipedia.org/wiki/Geographical_distance#Lambert's_formula_for_long_lines - NOTE: This algorithm uses geodesy/haversine_distance.py to compute central angle, - sigma - + NOTE: This algorithm uses geodesy/haversine_distance.py to compute the + central angle, sigma Representing the earth as an ellipsoid allows us to approximate distances between points on the surface much better than a sphere. Ellipsoidal formulas treat the Earth as an oblate ellipsoid which means accounting for the flattening that happens - at the North and South poles. Lambert's formulae provide accuracy on the order of - 10 meteres over thousands of kilometeres. Other methods can provide - millimeter-level accuracy but this is a simpler method to calculate long range + at the North and South poles. Lambert's formulas provide accuracy on the order of + 10 meters over thousands of kilometers. Other methods can provide + millimeter-level accuracy, but this is a simpler method to calculate long-range distances without increasing computational intensity. Args: lat1, lon1: latitude and longitude of coordinate 1 lat2, lon2: latitude and longitude of coordinate 2 Returns: - geographical distance between two points in metres + geographical distance between two points in meters >>> from collections import namedtuple >>> point_2d = namedtuple("point_2d", "lat lon") @@ -39,25 +38,20 @@ def lamberts_ellipsoidal_distance( >>> NEW_YORK = point_2d(40.713019, -74.012647) >>> VENICE = point_2d(45.443012, 12.313071) >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *YOSEMITE):0,.0f} meters" - '254,351 meters' + '254,032 meters' >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *NEW_YORK):0,.0f} meters" - '4,138,992 meters' + '4,133,295 meters' >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *VENICE):0,.0f} meters" - '9,737,326 meters' + '9,719,525 meters' """ # CONSTANTS per WGS84 https://en.wikipedia.org/wiki/World_Geodetic_System - # Distance in metres(m) - # Equation Parameters - # https://en.wikipedia.org/wiki/Geographical_distance#Lambert's_formula_for_long_lines flattening = (AXIS_A - AXIS_B) / AXIS_A # Parametric latitudes - # https://en.wikipedia.org/wiki/Latitude#Parametric_(or_reduced)_latitude b_lat1 = atan((1 - flattening) * tan(radians(lat1))) b_lat2 = atan((1 - flattening) * tan(radians(lat2))) - # Compute central angle between two points - # using haversine theta. sigma = haversine_distance / equatorial radius + # Compute the central angle between two points using the haversine function sigma = haversine_distance(lat1, lon1, lat2, lon2) / EQUATORIAL_RADIUS # Intermediate P and Q values @@ -65,13 +59,11 @@ def lamberts_ellipsoidal_distance( q_value = (b_lat2 - b_lat1) / 2 # Intermediate X value - # X = (sigma - sin(sigma)) * sin^2Pcos^2Q / cos^2(sigma/2) x_numerator = (sin(p_value) ** 2) * (cos(q_value) ** 2) - x_demonimator = cos(sigma / 2) ** 2 - x_value = (sigma - sin(sigma)) * (x_numerator / x_demonimator) + x_denominator = cos(sigma / 2) ** 2 + x_value = (sigma - sin(sigma)) * (x_numerator / x_denominator) # Intermediate Y value - # Y = (sigma + sin(sigma)) * cos^2Psin^2Q / sin^2(sigma/2) y_numerator = (cos(p_value) ** 2) * (sin(q_value) ** 2) y_denominator = sin(sigma / 2) ** 2 y_value = (sigma + sin(sigma)) * (y_numerator / y_denominator) From f468b672114694d527d62d9e4bc78dfb475431ec Mon Sep 17 00:00:00 2001 From: Maanit Date: Wed, 2 Oct 2024 03:20:36 +0530 Subject: [PATCH 6/6] lamberts_ellipsoidal_distance changes in accordance to new haversine distance parameters minor changes --- geodesy/lamberts_ellipsoidal_distance.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geodesy/lamberts_ellipsoidal_distance.py b/geodesy/lamberts_ellipsoidal_distance.py index aad04ce48..79ce62ed7 100644 --- a/geodesy/lamberts_ellipsoidal_distance.py +++ b/geodesy/lamberts_ellipsoidal_distance.py @@ -1,6 +1,6 @@ from math import atan, cos, radians, sin, tan -from haversine_distance import haversine_distance +from .haversine_distance import haversine_distance AXIS_A = 6378137.0 AXIS_B = 6356752.314245