fix: use geodetic latitudes in haversine distance formula (#14351)

* 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

* 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

* Fix typos

Updated the docstring for the haversine_distance function to improve clarity and fix minor grammatical issues.

* Fix typos in docstring and variable names

* Clarify note on using haversine_distance.py

Updated the note to clarify the use of haversine_distance.py.

---------

Co-authored-by: Christian Clauss <cclauss@me.com>

Author: VibhorGautam

geodesy/haversine_distance.py

@@ -1,28 +1,28 @@
-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
     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`, `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__":

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
@@ -12,18 +12,17 @@ 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,
-        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
+    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:
@@ -59,11 +58,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'
     """
 
     # Validate latitude values
@@ -86,7 +85,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
@@ -95,8 +94,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)