From 9fae7d9a89d2e885044fe21b5f196a1bfbc3553a Mon Sep 17 00:00:00 2001 From: VibhorGautam Date: Sun, 8 Mar 2026 23:46:06 +0530 Subject: [PATCH 1/2] 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/2] 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