fix: use geodetic latitudes in haversine distance formula

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)