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 the geodesic distances between points on Earth (i.e. WSG 84 ellipsoid) using Vincenty's inverse method.

Supports single or sequences of coordiates.

Parameters:

Name Type Description Default
a ArrayLike

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

 required
b ArrayLike

Second coordinates, same format/shape as a.

 required
units str

Output units, "km" (default) or "m".

 'km'
tolerance float

Convergence threshold in radians. Default is 1e-12.

 1e-12
max_iterations int

Maximum iterations before failure. Default is 10.

 10

Returns:

Type Description
float64 | NDArray[float64]

float or np.ndarray: The geodesic distance or distances between the point in a and b.

Raises:

Type Description
ValueError

If input shapes are incompatible or units are invalid.
Note

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

Examples:

>>> geodesic([51.352757, 12.43392], [38.559, 68.856])
4548.675334434374
>>> geodesic([0, 0], [[0, 0], [10, 0], [20, 0]])
array([   0.        , 1105.85483324, 2212.36625417])
>>> 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 the geodesic distances between points on Earth (i.e. WSG 84 ellipsoid) using Vincenty's inverse method.

    Supports single or sequences of coordiates.

    Args:
        a (ArrayLike): Coordinates [lat, lon] or array of shape (N, 2), in decimal degrees.
        b (ArrayLike): Second coordinates, same format/shape as `a`.
        units (str, optional): Output units, "km" (default) or "m".
        tolerance (float, optional): Convergence threshold in radians. Default is 1e-12.
        max_iterations (int, optional): Maximum iterations before failure. Default is 10.

    Returns:
        float or np.ndarray: The geodesic distance or distances between the point in `a` and `b`.

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

    Note:
        Uses WGS84 (a=6378137.0 m, f=1/298.257223563). May fail for nearly antipodal points.

    Examples:
        >>> geodesic([51.352757, 12.43392], [38.559, 68.856])
        4548.675334434374
        >>> geodesic([0, 0], [[0, 0], [10, 0], [20, 0]])
        array([   0.        , 1105.85483324, 2212.36625417])
        >>> 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 the great-circle (spherical) distance between pairs of latitude/longitude coordinates using the haversine formula.

Parameters:

Name Type Description Default
a ArrayLike

An array-like object of shape (..., 2) containing latitude and longitude coordinates in degrees. The last dimension must be 2: (lat, lon).

 required
b ArrayLike

An array-like object of the same shape as a, containing corresponding latitude and longitude coordinates.

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

Unit of the output distance. Must be either "km" for kilometers or "m" for meters. Defaults to "km".

 'km'
radius float

Radius of the sphere to use for distance calculation. Defaults to MEAN_EARTH_RADIUS_METERS (based on WSG 84 ellipsoid: ~6371008.77 meters). Note: If units="km", this value is automatically converted to kilometers.

 MEAN_EARTH_RADIUS_METERS

Returns:

Type Description

np.ndarray: Array of great-circle distances between a and b, in the specified units. The shape matches the input shape excluding the last dimension.

Raises:

Type Description
ValueError

If the shapes of a and b are incompatible or units is not one of "m" or "km".

Examples:

>>> haversine([51.352757, 12.43392], [38.559, 68.856])
4537.564747442274
>>> haversine([0, 0], [[0, 0], [10, 0], [20, 0]])
array([   0.        , 1111.95079735, 2223.90159469])
>>> 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
86
87
88
89
90
def haversine(
    a: ArrayLike,
    b: ArrayLike,
    units: Literal["m", "km"] = "km",
    radius: float = MEAN_EARTH_RADIUS_METERS,
):
    """
    Calculates the great-circle (spherical) distance between pairs of latitude/longitude coordinates
    using the haversine formula.

    Args:
        a (ArrayLike): An array-like object of shape (..., 2) containing latitude and longitude
            coordinates in degrees. The last dimension must be 2: (lat, lon).
        b (ArrayLike): An array-like object of the same shape as `a`, containing corresponding
            latitude and longitude coordinates.
        units (Literal["m", "km"], optional): Unit of the output distance. Must be either
            "km" for kilometers or "m" for meters. Defaults to "km".
        radius (float, optional): Radius of the sphere to use for distance calculation.
            Defaults to MEAN_EARTH_RADIUS_METERS (based on WSG 84 ellipsoid: ~6371008.77 meters).
            Note: If `units="km"`, this value is automatically converted to kilometers.

    Returns:
        np.ndarray: Array of great-circle distances between `a` and `b`, in the specified units.
            The shape matches the input shape excluding the last dimension.

    Raises:
        ValueError: If the shapes of `a` and `b` are incompatible or `units` is not one of "m" or "km".

    Examples:
        >>> haversine([51.352757, 12.43392], [38.559, 68.856])
        4537.564747442274
        >>> haversine([0, 0], [[0, 0], [10, 0], [20, 0]])
        array([   0.        , 1111.95079735, 2223.90159469])
        >>> 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 the geodesic distances between points on Earth (i.e. WSG 84 ellipsoid) using Vincenty's inverse method.

Supports single or sequences of coordiates.

Parameters:

Name Type Description Default
a ArrayLike

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

 required
b ArrayLike

Second coordinates, same format/shape as a.

 required
units str

Output units, "km" (default) or "m".

 'km'
tolerance float

Convergence threshold in radians. Default is 1e-12.

 1e-12
max_iterations int

Maximum iterations before failure. Default is 10.

 10

Returns:

Type Description
float64 | NDArray[float64]

float or np.ndarray: The geodesic distance or distances between the point in a and b.

Raises:

Type Description
ValueError

If input shapes are incompatible or units are invalid.
Note

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

Examples:

>>> geodesic([51.352757, 12.43392], [38.559, 68.856])
4548.675334434374
>>> geodesic([0, 0], [[0, 0], [10, 0], [20, 0]])
array([   0.        , 1105.85483324, 2212.36625417])
>>> 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 the geodesic distances between points on Earth (i.e. WSG 84 ellipsoid) using Vincenty's inverse method.

    Supports single or sequences of coordiates.

    Args:
        a (ArrayLike): Coordinates [lat, lon] or array of shape (N, 2), in decimal degrees.
        b (ArrayLike): Second coordinates, same format/shape as `a`.
        units (str, optional): Output units, "km" (default) or "m".
        tolerance (float, optional): Convergence threshold in radians. Default is 1e-12.
        max_iterations (int, optional): Maximum iterations before failure. Default is 10.

    Returns:
        float or np.ndarray: The geodesic distance or distances between the point in `a` and `b`.

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

    Note:
        Uses WGS84 (a=6378137.0 m, f=1/298.257223563). May fail for nearly antipodal points.

    Examples:
        >>> geodesic([51.352757, 12.43392], [38.559, 68.856])
        4548.675334434374
        >>> geodesic([0, 0], [[0, 0], [10, 0], [20, 0]])
        array([   0.        , 1105.85483324, 2212.36625417])
        >>> 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