TheAlgorithms-Python

1
from math import atan, cos, radians, sin, tan
2

3
from .haversine_distance import haversine_distance
4

5
AXIS_A = 6378137.0
6
AXIS_B = 6356752.314245
7
EQUATORIAL_RADIUS = 6378137
8

9

10
def lamberts_ellipsoidal_distance(
11
    lat1: float, lon1: float, lat2: float, lon2: float
12
) -> float:
13
    """
14
    Calculate the shortest distance along the surface of an ellipsoid between
15
    two points on the surface of earth given longitudes and latitudes
16
    https://en.wikipedia.org/wiki/Geographical_distance#Lambert's_formula_for_long_lines
17

18
    NOTE: This algorithm uses geodesy/haversine_distance.py to compute central angle,
19
        sigma
20

21
    Representing the earth as an ellipsoid allows us to approximate distances between
22
    points on the surface much better than a sphere. Ellipsoidal formulas treat the
23
    Earth as an oblate ellipsoid which means accounting for the flattening that happens
24
    at the North and South poles. Lambert's formulae provide accuracy on the order of
25
    10 meteres over thousands of kilometeres. Other methods can provide
26
    millimeter-level accuracy but this is a simpler method to calculate long range
27
    distances without increasing computational intensity.
28

29
    Args:
30
        lat1, lon1: latitude and longitude of coordinate 1
31
        lat2, lon2: latitude and longitude of coordinate 2
32
    Returns:
33
        geographical distance between two points in metres
34

35
    >>> from collections import namedtuple
36
    >>> point_2d = namedtuple("point_2d", "lat lon")
37
    >>> SAN_FRANCISCO = point_2d(37.774856, -122.424227)
38
    >>> YOSEMITE = point_2d(37.864742, -119.537521)
39
    >>> NEW_YORK = point_2d(40.713019, -74.012647)
40
    >>> VENICE = point_2d(45.443012, 12.313071)
41
    >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *YOSEMITE):0,.0f} meters"
42
    '254,351 meters'
43
    >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *NEW_YORK):0,.0f} meters"
44
    '4,138,992 meters'
45
    >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *VENICE):0,.0f} meters"
46
    '9,737,326 meters'
47
    """
48

49
    # CONSTANTS per WGS84 https://en.wikipedia.org/wiki/World_Geodetic_System
50
    # Distance in metres(m)
51
    # Equation Parameters
52
    # https://en.wikipedia.org/wiki/Geographical_distance#Lambert's_formula_for_long_lines
53
    flattening = (AXIS_A - AXIS_B) / AXIS_A
54
    # Parametric latitudes
55
    # https://en.wikipedia.org/wiki/Latitude#Parametric_(or_reduced)_latitude
56
    b_lat1 = atan((1 - flattening) * tan(radians(lat1)))
57
    b_lat2 = atan((1 - flattening) * tan(radians(lat2)))
58

59
    # Compute central angle between two points
60
    # using haversine theta. sigma =  haversine_distance / equatorial radius
61
    sigma = haversine_distance(lat1, lon1, lat2, lon2) / EQUATORIAL_RADIUS
62

63
    # Intermediate P and Q values
64
    p_value = (b_lat1 + b_lat2) / 2
65
    q_value = (b_lat2 - b_lat1) / 2
66

67
    # Intermediate X value
68
    # X = (sigma - sin(sigma)) * sin^2Pcos^2Q / cos^2(sigma/2)
69
    x_numerator = (sin(p_value) ** 2) * (cos(q_value) ** 2)
70
    x_demonimator = cos(sigma / 2) ** 2
71
    x_value = (sigma - sin(sigma)) * (x_numerator / x_demonimator)
72

73
    # Intermediate Y value
74
    # Y = (sigma + sin(sigma)) * cos^2Psin^2Q / sin^2(sigma/2)
75
    y_numerator = (cos(p_value) ** 2) * (sin(q_value) ** 2)
76
    y_denominator = sin(sigma / 2) ** 2
77
    y_value = (sigma + sin(sigma)) * (y_numerator / y_denominator)
78

79
    return EQUATORIAL_RADIUS * (sigma - ((flattening / 2) * (x_value + y_value)))
80

81

82
if __name__ == "__main__":
83
    import doctest
84

85
    doctest.testmod()
86