From 9fae7d9a89d2e885044fe21b5f196a1bfbc3553a Mon Sep 17 00:00:00 2001 From: VibhorGautam Date: Sun, 8 Mar 2026 23:46:06 +0530 Subject: [PATCH 1/5] fix: use geodetic latitudes in haversine distance formula The implementation was incorrectly using reduced latitudes (via a flattening factor from WGS84 ellipsoid constants) instead of raw geodetic latitudes. Reduced latitudes are appropriate for ellipsoidal models like Lambert's formula, but the Haversine formula operates on a sphere and should use geodetic latitudes directly. Changes: - Use radians(lat) directly instead of computing reduced latitudes with atan((1 - flattening) * tan(radians(lat))) - Replace equatorial radius (6378137m) with mean Earth radius (6371000m) for better spherical approximation - Remove unused WGS84 ellipsoid constants (AXIS_A, AXIS_B) - Remove unused imports (atan, tan) - Add edge case and cross-continental doctests Fixes #11308 --- geodesy/haversine_distance.py | 50 ++++++++++++++++++++++------------- 1 file changed, 32 insertions(+), 18 deletions(-) diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py index 39cd250af965..014a65451508 100644 --- a/geodesy/haversine_distance.py +++ b/geodesy/haversine_distance.py @@ -1,13 +1,11 @@ -from math import asin, atan, cos, radians, sin, sqrt, tan +from math import asin, cos, radians, sin, sqrt -AXIS_A = 6378137.0 -AXIS_B = 6356752.314245 -RADIUS = 6378137 +EARTH_RADIUS = 6371000 def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float: """ - Calculate great circle distance between two points in a sphere, + Calculate great circle distance between two points on 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 @@ -21,8 +19,10 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl 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 + lat1: latitude of coordinate 1 in degrees + lon1: longitude of coordinate 1 in degrees + lat2: latitude of coordinate 2 in degrees + lon2: longitude of coordinate 2 in degrees Returns: geographical distance between two points in metres @@ -31,25 +31,39 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> fl >>> 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' + '253,748 meters' + >>> NEW_YORK = point_2d(40.712776, -74.005974) + >>> LOS_ANGELES = point_2d(34.052235, -118.243683) + >>> f"{haversine_distance(*NEW_YORK, *LOS_ANGELES):0,.0f} meters" + '3,935,746 meters' + >>> LONDON = point_2d(51.507351, -0.127758) + >>> PARIS = point_2d(48.856614, 2.352222) + >>> f"{haversine_distance(*LONDON, *PARIS):0,.0f} meters" + '343,549 meters' + >>> haversine_distance(0, 0, 0, 0) + 0.0 + >>> from math import isclose + >>> quarter_equator = haversine_distance(0, 0, 0, 90) + >>> isclose(quarter_equator, 10_007_543, rel_tol=1e-3) + True """ - # 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))) + # Convert geodetic coordinates from degrees to radians. + # The Haversine formula operates on a sphere, so we use the raw geodetic + # latitudes directly rather than reduced latitudes (which apply to + # ellipsoidal models like Lambert's formula). + # Reference: https://en.wikipedia.org/wiki/Haversine_formula#Formulation + phi_1 = radians(lat1) + phi_2 = radians(lat2) lambda_1 = radians(lon1) lambda_2 = radians(lon2) - # Equation + + # Haversine 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) + return 2 * EARTH_RADIUS * asin(h_value) if __name__ == "__main__": From 635e9ca7f574459fe630d5b7a1c91532e41335b5 Mon Sep 17 00:00:00 2001 From: VibhorGautam Date: Sun, 8 Mar 2026 23:57:26 +0530 Subject: [PATCH 2/5] fix: update Lambert's to use corrected haversine radius for central angle Lambert's ellipsoidal distance computes the central angle sigma by dividing the haversine distance by a radius. Previously both functions used the same equatorial radius (6378137m), so the values cancelled out. After correcting haversine to use the mean Earth radius (6371000m), Lambert's must divide by the same radius to recover the correct angle. Also update the expected doctest values to match the corrected haversine output. Fixes #11308 --- geodesy/lamberts_ellipsoidal_distance.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/geodesy/lamberts_ellipsoidal_distance.py b/geodesy/lamberts_ellipsoidal_distance.py index 4805674e51ab..c1b1bbdf388b 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 EARTH_RADIUS, haversine_distance AXIS_A = 6378137.0 AXIS_B = 6356752.314245 @@ -39,11 +39,11 @@ 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 @@ -58,7 +58,7 @@ def lamberts_ellipsoidal_distance( # Compute central angle between two points # using haversine theta. sigma = haversine_distance / equatorial radius - sigma = haversine_distance(lat1, lon1, lat2, lon2) / EQUATORIAL_RADIUS + sigma = haversine_distance(lat1, lon1, lat2, lon2) / EARTH_RADIUS # Intermediate P and Q values p_value = (b_lat1 + b_lat2) / 2 From a0b8c6b1abcaa0b122d06c01a08bd021eb585efb Mon Sep 17 00:00:00 2001 From: Christian Clauss Date: Sat, 12 Sep 2026 09:56:54 +0200 Subject: [PATCH 3/5] Fix typos Updated the docstring for the haversine_distance function to improve clarity and fix minor grammatical issues. --- geodesy/haversine_distance.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py index 014a65451508..48ff8508fba5 100644 --- a/geodesy/haversine_distance.py +++ b/geodesy/haversine_distance.py @@ -5,18 +5,18 @@ def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float: """ - Calculate great circle distance between two points on a sphere, + Calculate great-circle distance between two points on 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 + 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. + Earth's ellipsoidal nature are more accurate, but a quick and modifiable + computation like Haversine can be handy for shorter-range distances. Args: lat1: latitude of coordinate 1 in degrees From fa54ec81042477ef87d0e4180ec59f6e19a24e72 Mon Sep 17 00:00:00 2001 From: Christian Clauss Date: Sat, 12 Sep 2026 09:58:18 +0200 Subject: [PATCH 4/5] Fix typos in docstring and variable names --- geodesy/lamberts_ellipsoidal_distance.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/geodesy/lamberts_ellipsoidal_distance.py b/geodesy/lamberts_ellipsoidal_distance.py index bc237709ca57..89bce5924a93 100644 --- a/geodesy/lamberts_ellipsoidal_distance.py +++ b/geodesy/lamberts_ellipsoidal_distance.py @@ -12,18 +12,18 @@ def lamberts_ellipsoidal_distance( ) -> float: """ Calculate the shortest distance along the surface of an ellipsoid between - two points on the surface of earth given longitudes and latitudes + 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, + 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 + 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 + 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 + 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: @@ -95,8 +95,8 @@ def lamberts_ellipsoidal_distance( # 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) From 9d503fcff52e5ee9f00d02892f6a2a05895ac378 Mon Sep 17 00:00:00 2001 From: Christian Clauss Date: Sat, 12 Sep 2026 10:00:41 +0200 Subject: [PATCH 5/5] Clarify note on using haversine_distance.py Updated the note to clarify the use of haversine_distance.py. --- geodesy/lamberts_ellipsoidal_distance.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/geodesy/lamberts_ellipsoidal_distance.py b/geodesy/lamberts_ellipsoidal_distance.py index 89bce5924a93..18be863f4a79 100644 --- a/geodesy/lamberts_ellipsoidal_distance.py +++ b/geodesy/lamberts_ellipsoidal_distance.py @@ -15,8 +15,7 @@ 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 the central angle, - sigma + NOTE: 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