Skip to content

fix: use geodetic latitudes in haversine distance formula #14351

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
VibhorGautam wants to merge 2 commits into
base: master
Choose a base branch
from
+37 −23
Open

fix: use geodetic latitudes in haversine distance formula #14351

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9fae7d9
fix: use geodetic latitudes in haversine distance formula
VibhorGautam Mar 8, 2026
635e9ca
fix: update Lambert's to use corrected haversine radius for central a…
VibhorGautam Mar 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
50 changes: 32 additions & 18 deletions geodesy/haversine_distance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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

Expand All @@ -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__":
Expand Down
10 changes: 5 additions & 5 deletions geodesy/lamberts_ellipsoidal_distance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down