Skip to content

API reference#

earthcarekit.geo.distance

Algorithms to calculate geospatial distances.

Notes#

This module depends on other internal modules:

geodesic #

geodesic(
    a: ArrayLike,
    b: ArrayLike,
    units: str = "km",
    tolerance: float = 1e-12,
    max_iterations: int = 10,
) -> float64 | NDArray[float64]

Calculates geodesic distances on the WGS 84 ellipsoid using Vincenty's inverse method.

Parameters:

Name Type Description Default
a ArrayLike

Coordinates [lat, lon] or (N, 2) array in decimal degrees.

 required
b ArrayLike

Second coordinates, same format/shape as a.

 required
units str

Output units ("km" or "m"); defaults to "km".

 'km'
tolerance float

Convergence threshold in radians; defaults to 1e-12.

 1e-12
max_iterations int

Maximum iterations before failure; defaults to 10.

 10

Returns:

Type Description
float64 | NDArray[float64]

Geodesic distance(s) in the specified units.

Raises:

Type Description
ValueError

If input shapes are incompatible, units are invalid, or convergence fails.
Notes

Uses WGS 84 parameters (a=6378137.0 m, f=1/298.257223563). May fail for nearly antipodal points.

Examples:

>>> import earthcarekit as eck
>>> eck.geodesic([51.352757, 12.43392], [38.559, 68.856])
4548.675334434374

>>> eck.geodesic([0, 0], [[0, 0], [10, 0], [20, 0]])
array([   0.        , 1105.85483324, 2212.36625417])

>>> eck.geodesic([[0, 0], [10, 0], [20, 0]], [[0, 0], [10, 0], [20, 0]])
array([0., 0., 0.])
References
  • Vincenty, T. (1975). Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations. Survey Review, 23(176), 88-93. https://doi.org/10.1179/sre.1975.23.176.88
Source code in earthcarekit/geo/distance/_vincenty.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def vincenty(
    a: ArrayLike,
    b: ArrayLike,
    units: str = "km",
    tolerance: float = 1e-12,
    max_iterations: int = 10,
) -> np.float64 | NDArray[np.float64]:
    """Calculates geodesic distances on the WGS 84 ellipsoid using Vincenty's inverse method.

    Args:
        a: Coordinates [lat, lon] or (N, 2) array in decimal degrees.
        b: Second coordinates, same format/shape as `a`.
        units: Output units ("km" or "m"); defaults to "km".
        tolerance: Convergence threshold in radians; defaults to 1e-12.
        max_iterations: Maximum iterations before failure; defaults to 10.

    Returns:
        Geodesic distance(s) in the specified units.

    Raises:
        ValueError: If input shapes are incompatible, units are invalid, or convergence fails.

    Notes:
        Uses WGS 84 parameters (a=6378137.0 m, f=1/298.257223563). May fail for nearly antipodal points.

    Examples:
        >>> import earthcarekit as eck
        >>> eck.geodesic([51.352757, 12.43392], [38.559, 68.856])
        4548.675334434374

        >>> eck.geodesic([0, 0], [[0, 0], [10, 0], [20, 0]])
        array([   0.        , 1105.85483324, 2212.36625417])

        >>> eck.geodesic([[0, 0], [10, 0], [20, 0]], [[0, 0], [10, 0], [20, 0]])
        array([0., 0., 0.])

    References:
        - Vincenty, T. (1975). Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations.
        Survey Review, 23(176), 88-93. https://doi.org/10.1179/sre.1975.23.176.88
    """
    _a, _b = map(np.asarray, [a, b])
    coord_a, coord_b = map(np.atleast_2d, [_a, _b])
    coord_a, coord_b = map(np.radians, [coord_a, coord_b])

    if (coord_a.shape[1] != 2) or (coord_b.shape[1] != 2):
        raise ValueError(
            f"At least one passed array has a wrong shape (a={_a.shape}, b={_b.shape}). 1d arrays should be of length 2 (i.e. [lat, lon]) and 2d array should have the shape (n, 2)."
        )
    if (coord_a.shape[0] < 1) or (coord_b.shape[0] < 1):
        raise ValueError(
            f"At least one passed array contains no values (a={_a.shape}, b={_b.shape})."
        )
    if coord_a.shape[0] != coord_b.shape[0]:
        if (coord_a.shape[0] != 1) and (coord_b.shape[0] != 1):
            raise ValueError(
                f"The shapes of passed arrays dont match (a={_a.shape}, b={_b.shape}). Either both should contain the same number of coordinates or at least one of them should contain a single coordinate."
            )

    lat_1, lon_1 = coord_a[:, 0], coord_a[:, 1]
    lat_2, lon_2 = coord_b[:, 0], coord_b[:, 1]

    # WGS84 ellipsoid constants
    a = 6378137.0  # semi-major axis (equatorial radius) in meters
    f = 1 / 298.257223563  # flattening
    b = (1 - f) * a  # semi-minor axis (polar radius) in meters

    # Reduced latitudes
    beta_1 = np.arctan((1 - f) * np.tan(lat_1))
    beta_2 = np.arctan((1 - f) * np.tan(lat_2))

    initial_lon_diff = lon_2 - lon_1

    # Initialize variables for iterative solution
    lon_diff = initial_lon_diff
    sin_beta_1, cos_beta_1 = np.sin(beta_1), np.cos(beta_1)
    sin_beta_2, cos_beta_2 = np.sin(beta_2), np.cos(beta_2)
    # Track convergence for each point pair
    converged = np.full_like(lat_1, False, dtype=bool)

    for _ in range(max_iterations):
        sin_lon_diff, cos_lon_diff = np.sin(lon_diff), np.cos(lon_diff)

        sin_sigma = np.sqrt(
            (cos_beta_2 * sin_lon_diff) ** 2
            + (cos_beta_1 * sin_beta_2 - sin_beta_1 * cos_beta_2 * cos_lon_diff) ** 2
        )
        cos_sigma = (sin_beta_1 * sin_beta_2) + (cos_beta_1 * cos_beta_2 * cos_lon_diff)
        sigma = np.arctan2(sin_sigma, cos_sigma)

        with warnings.catch_warnings():
            warnings.filterwarnings("ignore")
            sin_alpha = cos_beta_1 * cos_beta_2 * sin_lon_diff / sin_sigma
        sin_alpha = np.nan_to_num(sin_alpha, nan=0.0)
        cos2_alpha = 1 - sin_alpha**2

        with warnings.catch_warnings():
            warnings.filterwarnings("ignore")
            cos2_sigma_m = np.where(
                cos2_alpha != 0.0,
                cos_sigma - ((2 * sin_beta_1 * sin_beta_2) / cos2_alpha),
                0.0,
            )
        cos2_sigma_m = np.nan_to_num(cos2_sigma_m, nan=0.0)

        C = f / 16 * cos2_alpha * (4 + f * (4 - 3 * cos2_alpha))

        previous_lon_diff = lon_diff
        lon_diff = initial_lon_diff + (1 - C) * f * sin_alpha * (
            sigma + C * sin_sigma * (cos2_sigma_m + C * cos_sigma * (-1 + 2 * cos2_sigma_m**2))
        )
        converged = converged | (np.abs(lon_diff - previous_lon_diff) < tolerance)
        if np.all(converged):
            break

    u2 = cos2_alpha * (a**2 - b**2) / b**2
    A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
    B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)))

    delta_sigma = (
        B
        * sin_sigma
        * (
            cos2_sigma_m
            + B
            / 4
            * (
                cos_sigma * (-1 + 2 * cos2_sigma_m**2)
                - B / 6 * cos2_sigma_m * (-3 + 4 * sin_sigma**2) * (-3 + 4 * cos2_sigma_m**2)
            )
        )
    )

    distance = b * A * (sigma - delta_sigma)

    if units == "km":
        distance = distance / 1000.0
    elif units != "m":
        raise ValueError(f"{vincenty.__name__}() Invalid units : {units}. Use 'm' or 'km' instead.")

    if len(_a.shape) == 1 and len(_b.shape) == 1:
        return distance[0]

    return distance

haversine #

haversine(
    a: ArrayLike,
    b: ArrayLike,
    units: Literal["m", "km"] = "km",
    radius: float = MEAN_EARTH_RADIUS_METERS,
)

Calculates great-circle distances between latitude/longitude pairs using the haversine formula.

Parameters:

Name Type Description Default
a ArrayLike

Coordinates as (..., 2) array (lat, lon) in degrees.

 required
b ArrayLike

Corresponding coordinates as (..., 2) array (lat, lon) in degrees.

 required
units Literal['m', 'km']

Output unit ("m" or "km"); defaults to "km".

 'km'
radius float

Sphere radius in meters; defaults to mean earth radius (~6371008.77 m).

 MEAN_EARTH_RADIUS_METERS

Returns:

Type Description

Great-circle distances in the specified units; shape matches input excluding last dimension.

Raises:

Type Description
ValueError

If shapes of a and b are incompatible or units is invalid.

Examples:

>>> import earthcarekit as eck
>>> eck.haversine([51.352757, 12.43392], [38.559, 68.856])
4537.564747442274

>>> eck.haversine([0, 0], [[0, 0], [10, 0], [20, 0]])
array([   0.        , 1111.95079735, 2223.90159469])

>>> eck.haversine([[0, 0], [10, 0], [20, 0]], [[0, 0], [10, 0], [20, 0]])
array([0., 0., 0.])
Source code in earthcarekit/geo/distance/_haversine.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
def haversine(
    a: ArrayLike,
    b: ArrayLike,
    units: Literal["m", "km"] = "km",
    radius: float = MEAN_EARTH_RADIUS_METERS,
):
    """Calculates great-circle distances between latitude/longitude pairs using the haversine formula.

    Args:
        a: Coordinates as (..., 2) array (lat, lon) in degrees.
        b: Corresponding coordinates as (..., 2) array (lat, lon) in degrees.
        units: Output unit ("m" or "km"); defaults to "km".
        radius: Sphere radius in meters; defaults to mean earth radius (~6371008.77 m).

    Returns:
        Great-circle distances in the specified units; shape matches input excluding last dimension.

    Raises:
        ValueError: If shapes of `a` and `b` are incompatible or `units` is invalid.

    Examples:
        >>> import earthcarekit as eck
        >>> eck.haversine([51.352757, 12.43392], [38.559, 68.856])
        4537.564747442274

        >>> eck.haversine([0, 0], [[0, 0], [10, 0], [20, 0]])
        array([   0.        , 1111.95079735, 2223.90159469])

        >>> eck.haversine([[0, 0], [10, 0], [20, 0]], [[0, 0], [10, 0], [20, 0]])
        array([0., 0., 0.])
    """

    if units not in ["m", "km"]:
        raise ValueError(
            f"{haversine.__name__}() Invalid units : {units}. Use 'm' or 'km' instead."
        )

    if units == "km":
        radius = radius / 1000.0

    a = np.array(a)
    b = np.array(b)

    coord_a = np.atleast_2d(a)
    coord_b = np.atleast_2d(b)

    if (coord_a.shape[1] != 2) or (coord_b.shape[1] != 2):
        raise ValueError(
            f"At least one passed array has a wrong shape (a={a.shape}, b={b.shape}). 1d arrays should be of length 2 (i.e. [lat, lon]) and 2d array should have the shape (n, 2)."
        )
    if (coord_a.shape[0] < 1) or (coord_b.shape[0] < 1):
        raise ValueError(
            f"At least one passed array contains no values (a={a.shape}, b={b.shape})."
        )
    if coord_a.shape[0] != coord_b.shape[0]:
        if (coord_a.shape[0] != 1) and (coord_b.shape[0] != 1):
            raise ValueError(
                f"The shapes of passed arrays dont match (a={a.shape}, b={b.shape}). Either both should contain the same number of coordinates or at least one of them should contain a single coordinate."
            )

    coord_a = np.radians(coord_a)
    coord_b = np.radians(coord_b)

    phi_1, lambda_1 = coord_a[:, 0], coord_a[:, 1]
    phi_2, lambda_2 = coord_b[:, 0], coord_b[:, 1]

    def hav(theta):
        return (1 - np.cos(theta)) / 2

    h = hav(phi_2 - phi_1) + np.cos(phi_1) * np.cos(phi_2) * hav(lambda_2 - lambda_1)

    d = 2 * radius * np.arcsin(np.sqrt(h))

    if len(a.shape) == 1 and len(b.shape) == 1:
        return d[0]

    return d

vincenty #

vincenty(
    a: ArrayLike,
    b: ArrayLike,
    units: str = "km",
    tolerance: float = 1e-12,
    max_iterations: int = 10,
) -> float64 | NDArray[float64]

Calculates geodesic distances on the WGS 84 ellipsoid using Vincenty's inverse method.

Parameters:

Name Type Description Default
a ArrayLike

Coordinates [lat, lon] or (N, 2) array in decimal degrees.

 required
b ArrayLike

Second coordinates, same format/shape as a.

 required
units str

Output units ("km" or "m"); defaults to "km".

 'km'
tolerance float

Convergence threshold in radians; defaults to 1e-12.

 1e-12
max_iterations int

Maximum iterations before failure; defaults to 10.

 10

Returns:

Type Description
float64 | NDArray[float64]

Geodesic distance(s) in the specified units.

Raises:

Type Description
ValueError

If input shapes are incompatible, units are invalid, or convergence fails.
Notes

Uses WGS 84 parameters (a=6378137.0 m, f=1/298.257223563). May fail for nearly antipodal points.

Examples:

>>> import earthcarekit as eck
>>> eck.geodesic([51.352757, 12.43392], [38.559, 68.856])
4548.675334434374

>>> eck.geodesic([0, 0], [[0, 0], [10, 0], [20, 0]])
array([   0.        , 1105.85483324, 2212.36625417])

>>> eck.geodesic([[0, 0], [10, 0], [20, 0]], [[0, 0], [10, 0], [20, 0]])
array([0., 0., 0.])
References
  • Vincenty, T. (1975). Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations. Survey Review, 23(176), 88-93. https://doi.org/10.1179/sre.1975.23.176.88
Source code in earthcarekit/geo/distance/_vincenty.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def vincenty(
    a: ArrayLike,
    b: ArrayLike,
    units: str = "km",
    tolerance: float = 1e-12,
    max_iterations: int = 10,
) -> np.float64 | NDArray[np.float64]:
    """Calculates geodesic distances on the WGS 84 ellipsoid using Vincenty's inverse method.

    Args:
        a: Coordinates [lat, lon] or (N, 2) array in decimal degrees.
        b: Second coordinates, same format/shape as `a`.
        units: Output units ("km" or "m"); defaults to "km".
        tolerance: Convergence threshold in radians; defaults to 1e-12.
        max_iterations: Maximum iterations before failure; defaults to 10.

    Returns:
        Geodesic distance(s) in the specified units.

    Raises:
        ValueError: If input shapes are incompatible, units are invalid, or convergence fails.

    Notes:
        Uses WGS 84 parameters (a=6378137.0 m, f=1/298.257223563). May fail for nearly antipodal points.

    Examples:
        >>> import earthcarekit as eck
        >>> eck.geodesic([51.352757, 12.43392], [38.559, 68.856])
        4548.675334434374

        >>> eck.geodesic([0, 0], [[0, 0], [10, 0], [20, 0]])
        array([   0.        , 1105.85483324, 2212.36625417])

        >>> eck.geodesic([[0, 0], [10, 0], [20, 0]], [[0, 0], [10, 0], [20, 0]])
        array([0., 0., 0.])

    References:
        - Vincenty, T. (1975). Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations.
        Survey Review, 23(176), 88-93. https://doi.org/10.1179/sre.1975.23.176.88
    """
    _a, _b = map(np.asarray, [a, b])
    coord_a, coord_b = map(np.atleast_2d, [_a, _b])
    coord_a, coord_b = map(np.radians, [coord_a, coord_b])

    if (coord_a.shape[1] != 2) or (coord_b.shape[1] != 2):
        raise ValueError(
            f"At least one passed array has a wrong shape (a={_a.shape}, b={_b.shape}). 1d arrays should be of length 2 (i.e. [lat, lon]) and 2d array should have the shape (n, 2)."
        )
    if (coord_a.shape[0] < 1) or (coord_b.shape[0] < 1):
        raise ValueError(
            f"At least one passed array contains no values (a={_a.shape}, b={_b.shape})."
        )
    if coord_a.shape[0] != coord_b.shape[0]:
        if (coord_a.shape[0] != 1) and (coord_b.shape[0] != 1):
            raise ValueError(
                f"The shapes of passed arrays dont match (a={_a.shape}, b={_b.shape}). Either both should contain the same number of coordinates or at least one of them should contain a single coordinate."
            )

    lat_1, lon_1 = coord_a[:, 0], coord_a[:, 1]
    lat_2, lon_2 = coord_b[:, 0], coord_b[:, 1]

    # WGS84 ellipsoid constants
    a = 6378137.0  # semi-major axis (equatorial radius) in meters
    f = 1 / 298.257223563  # flattening
    b = (1 - f) * a  # semi-minor axis (polar radius) in meters

    # Reduced latitudes
    beta_1 = np.arctan((1 - f) * np.tan(lat_1))
    beta_2 = np.arctan((1 - f) * np.tan(lat_2))

    initial_lon_diff = lon_2 - lon_1

    # Initialize variables for iterative solution
    lon_diff = initial_lon_diff
    sin_beta_1, cos_beta_1 = np.sin(beta_1), np.cos(beta_1)
    sin_beta_2, cos_beta_2 = np.sin(beta_2), np.cos(beta_2)
    # Track convergence for each point pair
    converged = np.full_like(lat_1, False, dtype=bool)

    for _ in range(max_iterations):
        sin_lon_diff, cos_lon_diff = np.sin(lon_diff), np.cos(lon_diff)

        sin_sigma = np.sqrt(
            (cos_beta_2 * sin_lon_diff) ** 2
            + (cos_beta_1 * sin_beta_2 - sin_beta_1 * cos_beta_2 * cos_lon_diff) ** 2
        )
        cos_sigma = (sin_beta_1 * sin_beta_2) + (cos_beta_1 * cos_beta_2 * cos_lon_diff)
        sigma = np.arctan2(sin_sigma, cos_sigma)

        with warnings.catch_warnings():
            warnings.filterwarnings("ignore")
            sin_alpha = cos_beta_1 * cos_beta_2 * sin_lon_diff / sin_sigma
        sin_alpha = np.nan_to_num(sin_alpha, nan=0.0)
        cos2_alpha = 1 - sin_alpha**2

        with warnings.catch_warnings():
            warnings.filterwarnings("ignore")
            cos2_sigma_m = np.where(
                cos2_alpha != 0.0,
                cos_sigma - ((2 * sin_beta_1 * sin_beta_2) / cos2_alpha),
                0.0,
            )
        cos2_sigma_m = np.nan_to_num(cos2_sigma_m, nan=0.0)

        C = f / 16 * cos2_alpha * (4 + f * (4 - 3 * cos2_alpha))

        previous_lon_diff = lon_diff
        lon_diff = initial_lon_diff + (1 - C) * f * sin_alpha * (
            sigma + C * sin_sigma * (cos2_sigma_m + C * cos_sigma * (-1 + 2 * cos2_sigma_m**2))
        )
        converged = converged | (np.abs(lon_diff - previous_lon_diff) < tolerance)
        if np.all(converged):
            break

    u2 = cos2_alpha * (a**2 - b**2) / b**2
    A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
    B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)))

    delta_sigma = (
        B
        * sin_sigma
        * (
            cos2_sigma_m
            + B
            / 4
            * (
                cos_sigma * (-1 + 2 * cos2_sigma_m**2)
                - B / 6 * cos2_sigma_m * (-3 + 4 * sin_sigma**2) * (-3 + 4 * cos2_sigma_m**2)
            )
        )
    )

    distance = b * A * (sigma - delta_sigma)

    if units == "km":
        distance = distance / 1000.0
    elif units != "m":
        raise ValueError(f"{vincenty.__name__}() Invalid units : {units}. Use 'm' or 'km' instead.")

    if len(_a.shape) == 1 and len(_b.shape) == 1:
        return distance[0]

    return distance