Skip to content

API reference#

earthcarekit.geo

Geospatial utilities for handling coordinates, coordinate transformations, distance calculation, and interpolation.

Notes#

This module depends on other internal modules:

create_spherical_grid #

create_spherical_grid(
    nlat: int,
    nlon: int | None = None,
    reduced: bool = False,
    lat_spacing: Literal["regular", "sinusoidal", "gaussian"] = "regular",
    min_lat: float = -90.0,
    max_lat: float = 90.0,
    min_lon: float = -180.0,
    max_lon: float = 180.0,
) -> SphericalGrid

Generate a spherical gird with regular, sinusoidal or gaussian latitudes and uniform or reduced longitudes.

Parameters:

Name Type Description Default
nlat int

Number of latitude.

 required
nlon int | None

Nuber of longitudes at the equator. If None, set to nlat * 2. Defaults to None.

 None
reduced bool

If True, reduces longitudes near poles using ~cos(latitude) scaling. Defaults to False.

 False
lat_spacing regular or sinusoidal or gaussian

Method used to place latitudes. Defaults to "regular".

 'regular'

Returns:

Name Type Description
SphericalGrid SphericalGrid

A container storing
SphericalGrid
  • lat_centers: A numpy.array of latitude bin centers.
SphericalGrid
  • lon_centers: A numpy.array of longitude bin centers or if is_reduced=True a list of numpy.array per latitude center.
SphericalGrid
  • lat_bounds: A numpy.array of latitude bin bounds.
SphericalGrid
  • lon_bounds: A numpy.array of longitude bin bounds or if is_reduced=True a list of numpy.array per latitude bound.
SphericalGrid
  • is_reduced (bool): Whether the grid has reduced number of longitudes near the poles.
Source code in earthcarekit/geo/grid/_create_global_grid.py 
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
def create_spherical_grid(
    nlat: int,
    nlon: int | None = None,
    reduced: bool = False,
    lat_spacing: Literal["regular", "sinusoidal", "gaussian"] = "regular",
    min_lat: float = -90.0,
    max_lat: float = 90.0,
    min_lon: float = -180.0,
    max_lon: float = 180.0,
) -> SphericalGrid:
    """
    Generate a spherical gird with regular, sinusoidal or gaussian latitudes and uniform or reduced longitudes.

    Args:
        nlat (int): Number of latitude.
        nlon (int | None): Nuber of longitudes at the equator. If None, set to `nlat * 2`. Defaults to None.
        reduced (bool): If True, reduces longitudes near poles using `~cos(latitude)` scaling. Defaults to False.
        lat_spacing ("regular" or "sinusoidal" or "gaussian", optional): Method used to place latitudes. Defaults to "regular".

    Returns:
        SphericalGrid: A container storing

        - lat_centers: A `numpy.array` of latitude bin centers.
        - lon_centers: A `numpy.array` of longitude bin centers or if `is_reduced=True` a list of `numpy.array` per latitude center.
        - lat_bounds: A `numpy.array` of latitude bin bounds.
        - lon_bounds: A `numpy.array` of longitude bin bounds or if `is_reduced=True` a list of `numpy.array` per latitude bound.
        - is_reduced (bool): Whether the grid has reduced number of longitudes near the poles.
    """

    if nlon is None:
        nlon = nlat + nlat

    if lat_spacing == "regular":
        lat_c, lat_b = _get_regular_latitudes(nlat, min_lat=min_lat, max_lat=max_lat)
    elif lat_spacing == "sinusoidal":
        lat_c, lat_b = _get_sinusoidal_latitudes(nlat, min_lat=min_lat, max_lat=max_lat)
    elif lat_spacing == "gaussian":
        lat_c, lat_b = _get_gaussian_latitudes(nlat, min_lat=min_lat, max_lat=max_lat)
    else:
        raise ValueError("grid_type must be 'regular', 'gaussian', or 'sinusoidal'")

    lon_b_list = _get_lon_bounds_per_lat(
        nlon,
        lat_c,
        reduced,
        min_lon=min_lon,
        max_lon=max_lon,
    )

    lon_b: NDArray | list[NDArray]
    if reduced:
        lon_b = lon_b_list
        lon_c = [0.5 * (lb[:-1] + lb[1:]) for lb in lon_b_list]
    else:
        lon_b = lon_b_list[0]
        lon_c = 0.5 * (lon_b[:-1] + lon_b[1:])

    return SphericalGrid(
        lat_centers=lat_c,
        lon_centers=lon_c,
        lat_bounds=lat_b,
        lon_bounds=lon_b,
        is_reduced=reduced,
    )

ecef_to_geo #

ecef_to_geo(
    x: SupportsFloat,
    y: SupportsFloat,
    z: SupportsFloat,
    target_radius: float = 1.0,
    perfect_sphere: bool = True,
    semi_major: float = SEMI_MAJOR_AXIS_METERS,
    semi_minor: float = SEMI_MINOR_AXIS_METERS,
) -> tuple[float, float, float]

Converts Earth-centered, Earth-fixed (ECEF) coordinates (x, y, z) to geodetic coordinates.

Parameters:

Name Type Description Default
x SupportsFloat

Cartesian ECEF x-coordinate in meters.

 required
y SupportsFloat

Cartesian ECEF y-coordinate in meters.

 required
z SupportsFloat

Cartesian ECEF z-coordinate in meters.

 required
target_radius float

Target mean radius of the ECEF coordinate system; defaults to 1.0.

 1.0
perfect_sphere bool

If True, assume a spherical Earth; otherwise, use ellipsoidal (WGS 84).

 True
semi_major float

Semi-major axis of the ellipsoid in meters; defaults to WGS 84 (6378137).

 SEMI_MAJOR_AXIS_METERS
semi_minor float

Semi-minor axis of the ellipsoid in meters; defaults to WGS 84 (6356752.314245).

 SEMI_MINOR_AXIS_METERS

Returns:

Type Description
tuple[float, float, float]

Geodetic coordinates (latitude, longitude, altitude)

  • lat: Latitude in degrees (positive north, negative south).
  • lon: Longitude in degrees (positive east, negative west).
  • alt: Height above the ellipsoid in meters.
Source code in earthcarekit/geo/convertsions.py 
 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
150
151
152
def ecef_to_geo(
    x: SupportsFloat,
    y: SupportsFloat,
    z: SupportsFloat,
    target_radius: float = 1.0,
    perfect_sphere: bool = True,
    semi_major: float = SEMI_MAJOR_AXIS_METERS,
    semi_minor: float = SEMI_MINOR_AXIS_METERS,
) -> tuple[float, float, float]:
    """Converts Earth-centered, Earth-fixed (ECEF) coordinates (x, y, z) to geodetic coordinates.

    Args:
        x: Cartesian ECEF x-coordinate in meters.
        y: Cartesian ECEF y-coordinate in meters.
        z: Cartesian ECEF z-coordinate in meters.
        target_radius: Target mean radius of the ECEF coordinate system; defaults to 1.0.
        perfect_sphere: If True, assume a spherical Earth; otherwise, use ellipsoidal (WGS 84).
        semi_major: Semi-major axis of the ellipsoid in meters; defaults to WGS 84 (6378137).
        semi_minor: Semi-minor axis of the ellipsoid in meters; defaults to WGS 84 (6356752.314245).

    Returns:
        Geodetic coordinates (latitude, longitude, altitude)

            - lat: Latitude in degrees (positive north, negative south).
            - lon: Longitude in degrees (positive east, negative west).
            - alt: Height above the ellipsoid in meters.
    """
    x = float(x)
    y = float(y)
    z = float(z)

    # Undo scaling
    R = ((semi_major + (semi_major if perfect_sphere else semi_minor)) / 2) / target_radius
    x = -x * R
    y = -y * R
    z = z * R

    lon = math.atan2(y, x)

    if perfect_sphere:
        r = math.sqrt(x**2 + y**2 + z**2)
        lat = math.asin(z / r)
        alt = r - semi_major
    else:
        e2 = 1 - (semi_minor**2 / semi_major**2)
        p = math.sqrt(x**2 + y**2)
        # Initial guess
        lat = math.atan2(z, p * (1 - e2))
        # Iterative refienment
        for _ in range(5):
            N = semi_major / math.sqrt(1 - e2 * math.sin(lat) ** 2)
            alt = p / math.cos(lat) - N
            lat = math.atan2(z, p * (1 - e2 * (N / (N + alt))))
        N = semi_major / math.sqrt(1 - e2 * math.sin(lat) ** 2)
        alt = p / math.cos(lat) - N

    return math.degrees(lat), math.degrees(lon), alt

geo_to_ecef #

geo_to_ecef(
    lat: SupportsFloat,
    lon: SupportsFloat,
    alt: SupportsFloat | None = None,
    target_radius: float = 1.0,
    perfect_sphere: bool = True,
    semi_major: float = SEMI_MAJOR_AXIS_METERS,
    semi_minor: float = SEMI_MINOR_AXIS_METERS,
) -> tuple[float, float, float]

Converts geodetic coordinates (latitude, longitude, altitude) to Earth-centered, Earth-fixed (ECEF) coordinates.

Parameters:

Name Type Description Default
lat SupportsFloat

Latitude in degrees (positive north, negative south).

 required
lon SupportsFloat

Longitude in degrees (positive east, negative west).

 required
alt SupportsFloat | None

Height above the Earth ellipsoid in meters; None assumes ellipsoid surface.

 None
target_radius float

Target mean radius of the ECEF coordinate system; defaults to 1.0.

 1.0
perfect_sphere bool

If True, use a sphere with radius target_radius; otherwise, use ellipsoid.

 True
semi_major float

Semi-major axis of the ellipsoid in meters; defaults to WGS 84 (6378137).

 SEMI_MAJOR_AXIS_METERS
semi_minor float

Semi-minor axis of the ellipsoid in meters; defaults to WGS 84 (6356752.314245).

 SEMI_MINOR_AXIS_METERS

Returns:

Type Description
tuple[float, float, float]

ECEF coordinates (x, y, z) in meters

  • x: Axis through equator at prime meridian (lat=0°, lon=0°).
  • y: Axis through equator at 90°E (lat=0°, lon=90°).
  • z: Axis through north pole (lat=90°).
Source code in earthcarekit/geo/convertsions.py 
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
def geo_to_ecef(
    lat: SupportsFloat,
    lon: SupportsFloat,
    alt: SupportsFloat | None = None,
    target_radius: float = 1.0,
    perfect_sphere: bool = True,
    semi_major: float = SEMI_MAJOR_AXIS_METERS,
    semi_minor: float = SEMI_MINOR_AXIS_METERS,
) -> tuple[float, float, float]:
    """Converts geodetic coordinates (latitude, longitude, altitude) to Earth-centered, Earth-fixed (ECEF) coordinates.

    Args:
        lat: Latitude in degrees (positive north, negative south).
        lon: Longitude in degrees (positive east, negative west).
        alt: Height above the Earth ellipsoid in meters; None assumes ellipsoid surface.
        target_radius: Target mean radius of the ECEF coordinate system; defaults to 1.0.
        perfect_sphere: If True, use a sphere with radius `target_radius`; otherwise, use ellipsoid.
        semi_major: Semi-major axis of the ellipsoid in meters; defaults to WGS 84 (6378137).
        semi_minor: Semi-minor axis of the ellipsoid in meters; defaults to WGS 84 (6356752.314245).

    Returns:
        ECEF coordinates (x, y, z) in meters

            - x: Axis through equator at prime meridian (lat=0°, lon=0°).
            - y: Axis through equator at 90°E (lat=0°, lon=90°).
            - z: Axis through north pole (lat=90°).
    """
    lat = float(lat)
    lon = float(lon)
    if alt is None:
        alt = 0.0
    else:
        alt = float(alt)

    def sin(x):
        return math.sin(x)

    def cos(x):
        return math.cos(x)

    def sqrt(x):
        return math.sqrt(x)

    lat = math.radians(lat)
    lon = math.radians(lon)

    if perfect_sphere:
        # Calculate ECEF coordinates
        f = 1 - (semi_major / semi_major)  # Flattening of the ellipsoid
        N = semi_major / sqrt(1 - (f * sin(lat)) ** 2)  # Prime vertical radius of curvature

        x = (N + alt) * cos(lat) * cos(lon)
        y = (N + alt) * cos(lat) * sin(lon)
        z = ((semi_major**2 / semi_major**2) * N + alt) * sin(lat)

        # # Alternative
        # e2 = 1 - (semi_minor**2 / semi_major**2) # Square of the first numerical eccentricity of the ellipsoid
        # N = semi_major / sqrt(1 - e2 * (sin(lat) ** 2)) # Prime vertical radius of curvature

        # Scale ECEF coordinates to target radius
        R = ((semi_major + semi_major) / 2) / target_radius
        x = -x / R
        y = -y / R
        z = z / R
    else:
        # Calculate ECEF coordinates
        f = 1 - (semi_minor / semi_major)  # Flattening of the ellipsoid
        N = semi_major / sqrt(1 - (f * sin(lat)) ** 2)  # Prime vertical radius of curvature

        x = (N + alt) * cos(lat) * cos(lon)
        y = (N + alt) * cos(lat) * sin(lon)
        z = ((semi_minor**2 / semi_major**2) * N + alt) * sin(lat)

        # # Alternative
        # e2 = 1 - (semi_minor**2 / semi_major**2) # Square of the first numerical eccentricity of the ellipsoid
        # N = semi_major / sqrt(1 - e2 * (sin(lat) ** 2)) # Prime vertical radius of curvature

        # Scale ECEF coordinates to target radius
        R = ((semi_major + semi_minor) / 2) / target_radius
        x = -x / R
        y = -y / R
        z = z / R

    return x, y, z

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

get_central_coords #

get_central_coords(latitude: ArrayLike, longitude: ArrayLike) -> tuple[float, float]

Calculates the central lat/lon coordinates.

Source code in earthcarekit/geo/coordinates.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
def get_central_coords(
    latitude: ArrayLike,
    longitude: ArrayLike,
) -> tuple[float, float]:
    """Calculates the central lat/lon coordinates."""
    from .convertsions import ecef_to_geo, geo_to_ecef

    lats: NDArray = flatten_array(latitude)
    lons: NDArray = flatten_array(longitude)

    coords_ecef = np.array([geo_to_ecef(lat=lt, lon=ln) for lt, ln in zip(lats, lons)])
    coords_ecef_min = np.nanmin(coords_ecef, axis=0)
    coords_ecef_max = np.nanmax(coords_ecef, axis=0)
    coords_ecef_central = (coords_ecef_min + coords_ecef_max) * 0.5
    coords_geo_central = ecef_to_geo(
        coords_ecef_central[0], coords_ecef_central[1], coords_ecef_central[2]
    )
    return (coords_geo_central[0], coords_geo_central[1])

get_central_latitude #

get_central_latitude(latitude: ArrayLike) -> float

Calculates the central latitude coordinate.

Source code in earthcarekit/geo/coordinates.py 
59
60
61
62
63
64
65
66
67
68
def get_central_latitude(
    latitude: ArrayLike,
) -> float:
    """Calculates the central latitude coordinate."""
    lats: NDArray = flatten_array(latitude)
    lons: NDArray = np.zeros(lats.shape)

    central_coords = get_central_coords(lats, lons)

    return central_coords[0]

get_central_longitude #

get_central_longitude(longitude: ArrayLike) -> float

Calculates the central longitude coordinate.

Source code in earthcarekit/geo/coordinates.py 
71
72
73
74
75
76
77
78
79
80
def get_central_longitude(
    longitude: ArrayLike,
) -> float:
    """Calculates the central longitude coordinate."""
    lons: NDArray = flatten_array(longitude)
    lats: NDArray = np.zeros(lons.shape)

    central_coords = get_central_coords(lats, lons)

    return central_coords[1]

get_coord_between #

get_coord_between(coord1: ArrayLike, coord2: ArrayLike, f: float = 0.5) -> NDArray

Interpolates between two coordinates by fraction f (0 to 1).

Parameters:

Name Type Description Default
coord1 ArrayLike

The first lat/lon point.

 required
coord2 ArrayLike

The second lat/lon point.

 required
f float

A fractional value between 0 and 1. Defaults to 0.5, i.e., the mid point between coord1 and coord2.

 0.5

Returns:

Type Description
NDArray

A 2-element numpy.ndarray representing the interpolated lat/lon point.
Source code in earthcarekit/geo/interpolate.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
def get_coord_between(
    coord1: ArrayLike,
    coord2: ArrayLike,
    f: float = 0.5,
) -> NDArray:
    """Interpolates between two coordinates by fraction f (0 to 1).

    Args:
        coord1: The first lat/lon point.
        coord2: The second lat/lon point.
        f: A fractional value between 0 and 1. Defaults to 0.5, i.e., the mid point between coord1 and coord2.

    Returns:
        A 2-element `numpy.ndarray` representing the interpolated lat/lon point.
    """

    coord1 = np.array(coord1)
    coord2 = np.array(coord2)

    if coord1.shape != (2,):
        raise ValueError("coord1 must be a 2-element sequence (lat, lon)")

    if coord2.shape != (2,):
        raise ValueError("coord2 must be a 2-element sequence (lat, lon)")

    lon, lat = interpgeo(
        lat1=float(coord1[0]),
        lon1=float(coord1[1]),
        lat2=float(coord2[0]),
        lon2=float(coord2[1]),
        f=f,
    )
    return np.array([lat, lon])

get_coords #

get_coords(
    ds: Dataset,
    *,
    lat_var: str = TRACK_LAT_VAR,
    lon_var: str = TRACK_LON_VAR,
    flatten: bool = False
) -> NDArray

Takes a xarray.Dataset and returns the lat/lon coordinates as a numpy array.

Parameters:

Name Type Description Default
lat_var str

Name of the latitude variable.

 TRACK_LAT_VAR
lon_var str

Name of the longitude variable.

 TRACK_LON_VAR
flatten bool

If True, the coordinates will be flattened to a 2D array

  • 1st dimension: time
  • 2nd dimension: lat/lon
 False

Returns:

Type Description
NDArray

The extracted lat/lon coordinates.
Source code in earthcarekit/geo/coordinates.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
def get_coords(
    ds: xr.Dataset,
    *,
    lat_var: str = TRACK_LAT_VAR,
    lon_var: str = TRACK_LON_VAR,
    flatten: bool = False,
) -> NDArray:
    """Takes a `xarray.Dataset` and returns the lat/lon coordinates as a numpy array.

    Args:
        lat_var: Name of the latitude variable.
        lon_var: Name of the longitude variable.
        flatten:
            If True, the coordinates will be flattened to a 2D array

            - 1st dimension: time
            - 2nd dimension: lat/lon

    Returns:
        The extracted lat/lon coordinates.
    """
    lat = ds[lat_var].values
    lon = ds[lon_var].values
    coords = np.stack((lat, lon)).transpose()

    if len(coords.shape) > 2 and flatten:
        coords = coords.reshape(-1, 2)
    return coords

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

interpgeo #

interpgeo(lat1: float, lon1: float, lat2: float, lon2: float, f: float) -> tuple[float, float]

Interpolates along the geodesic from (lon1, lat1) to (lon2, lat2) by fraction f (0 to 1) and returns interpolated (lon, lat).

Source code in earthcarekit/geo/interpolate.py 
42
43
44
45
46
47
48
49
50
51
52
53
54
def interpgeo(
    lat1: float,
    lon1: float,
    lat2: float,
    lon2: float,
    f: float,
) -> tuple[float, float]:
    """
    Interpolates along the geodesic from (lon1, lat1) to (lon2, lat2) by fraction f (0 to 1) and returns interpolated (lon, lat).
    """
    azi1, azi2, dist = GEOD.inv(lon1, lat1, lon2, lat2)
    lon, lat, _ = GEOD.fwd(lon1, lat1, azi1, f * dist)
    return lon, lat

sequence_ecef_to_geo #

sequence_ecef_to_geo(
    x: NDArray | list[SupportsFloat],
    y: NDArray | list[SupportsFloat],
    z: NDArray | list[SupportsFloat],
    target_radius: float = 1.0,
    perfect_sphere: bool = True,
    semi_major: float = SEMI_MAJOR_AXIS_METERS,
    semi_minor: float = SEMI_MINOR_AXIS_METERS,
) -> NDArray

Converts sequences of Earth-centered, Earth-fixed (ECEF) coordinates to geodetic coordinates.

Parameters:

Name Type Description Default
x NDArray | list[SupportsFloat]

ECEF x-coordinates in meters.

 required
y NDArray | list[SupportsFloat]

ECEF y-coordinates in meters.

 required
z NDArray | list[SupportsFloat]

ECEF z-coordinates in meters.

 required
target_radius float

Target mean radius of the ECEF coordinate system; defaults to 1.0.

 1.0
perfect_sphere bool

If True, assume a spherical Earth; otherwise, use ellipsoidal (WGS 84).

 True
semi_major float

Semi-major axis of the ellipsoid in meters; defaults to WGS 84 (6378137).

 SEMI_MAJOR_AXIS_METERS
semi_minor float

Semi-minor axis of the ellipsoid in meters; defaults to WGS 84 (6356752.314245).

 SEMI_MINOR_AXIS_METERS

Returns:

Type Description
NDArray

Geodetic coordinates as a 2D array (N, 3); (latitude, longitude, altitude) for each input point.
Source code in earthcarekit/geo/convertsions.py 
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
def sequence_ecef_to_geo(
    x: NDArray | list[SupportsFloat],
    y: NDArray | list[SupportsFloat],
    z: NDArray | list[SupportsFloat],
    target_radius: float = 1.0,
    perfect_sphere: bool = True,
    semi_major: float = SEMI_MAJOR_AXIS_METERS,
    semi_minor: float = SEMI_MINOR_AXIS_METERS,
) -> NDArray:
    """Converts sequences of Earth-centered, Earth-fixed (ECEF) coordinates to geodetic coordinates.

    Args:
        x: ECEF x-coordinates in meters.
        y: ECEF y-coordinates in meters.
        z: ECEF z-coordinates in meters.
        target_radius: Target mean radius of the ECEF coordinate system; defaults to 1.0.
        perfect_sphere: If True, assume a spherical Earth; otherwise, use ellipsoidal (WGS 84).
        semi_major: Semi-major axis of the ellipsoid in meters; defaults to WGS 84 (6378137).
        semi_minor: Semi-minor axis of the ellipsoid in meters; defaults to WGS 84 (6356752.314245).

    Returns:
        Geodetic coordinates as a 2D array (N, 3); (latitude, longitude, altitude) for each input point.
    """
    x = np.asarray(x)
    y = np.asarray(y)
    z = np.asarray(z)

    xyz: NDArray = np.stack((x, y, z)).T
    coords = np.array(
        [
            list(
                ecef_to_geo(
                    c[0],
                    c[1],
                    c[2],
                    target_radius=target_radius,
                    perfect_sphere=perfect_sphere,
                    semi_major=semi_major,
                    semi_minor=semi_minor,
                )
            )
            for c in xyz
        ]
    )
    return coords

sequence_geo_to_ecef #

sequence_geo_to_ecef(
    lats: NDArray | list[SupportsFloat],
    lons: NDArray | list[SupportsFloat],
    alts: NDArray | list[SupportsFloat] | SupportsFloat | None = None,
    target_radius: float = 1.0,
    perfect_sphere: bool = True,
    semi_major: float = SEMI_MAJOR_AXIS_METERS,
    semi_minor: float = SEMI_MINOR_AXIS_METERS,
) -> NDArray

Converts sequences of geodetic coordinates to Earth-centered, Earth-fixed (ECEF) coordinates.

Parameters:

Name Type Description Default
lats NDArray | list[SupportsFloat]

Latitude values in degrees (positive north, negative south).

 required
lons NDArray | list[SupportsFloat]

Longitude values in degrees (positive east, negative west).

 required
alts NDArray | list[SupportsFloat] | SupportsFloat | None

Height above the ellipsoid in meters; None assumes ellipsoid surface.

 None
target_radius float

Target mean radius of the ECEF coordinate system; defaults to 1.0.

 1.0
perfect_sphere bool

If True, assume a spherical Earth; otherwise, use ellipsoidal (WGS 84).

 True
semi_major float

Semi-major axis of the ellipsoid in meters; defaults to WGS 84 (6378137).

 SEMI_MAJOR_AXIS_METERS
semi_minor float

Semi-minor axis of the ellipsoid in meters; defaults to WGS 84 (6356752.314245).

 SEMI_MINOR_AXIS_METERS

Returns:

Type Description
NDArray

ECEF coordinates as a 2D array (N, 3), where N is the number of input points.
Source code in earthcarekit/geo/convertsions.py 
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
def sequence_geo_to_ecef(
    lats: NDArray | list[SupportsFloat],
    lons: NDArray | list[SupportsFloat],
    alts: NDArray | list[SupportsFloat] | SupportsFloat | None = None,
    target_radius: float = 1.0,
    perfect_sphere: bool = True,
    semi_major: float = SEMI_MAJOR_AXIS_METERS,
    semi_minor: float = SEMI_MINOR_AXIS_METERS,
) -> NDArray:
    """Converts sequences of geodetic coordinates to Earth-centered, Earth-fixed (ECEF) coordinates.

    Args:
        lats: Latitude values in degrees (positive north, negative south).
        lons: Longitude values in degrees (positive east, negative west).
        alts: Height above the ellipsoid in meters; None assumes ellipsoid surface.
        target_radius: Target mean radius of the ECEF coordinate system; defaults to 1.0.
        perfect_sphere: If True, assume a spherical Earth; otherwise, use ellipsoidal (WGS 84).
        semi_major: Semi-major axis of the ellipsoid in meters; defaults to WGS 84 (6378137).
        semi_minor: Semi-minor axis of the ellipsoid in meters; defaults to WGS 84 (6356752.314245).

    Returns:
        ECEF coordinates as a 2D array (N, 3), where N is the number of input points.
    """
    lats = np.asarray(lats)
    lons = np.asarray(lons)
    if alts is None:
        alts = np.zeros(lats.shape)
    elif isinstance(alts, (float, int)):
        alts = np.repeat(float(alts), lats.shape)
    else:
        alts = np.asarray(alts)
    coords: NDArray = np.stack((lats, lons, alts)).T
    xyz = np.array(
        [
            list(
                geo_to_ecef(
                    c[0],
                    c[1],
                    c[2],
                    target_radius=target_radius,
                    perfect_sphere=perfect_sphere,
                    semi_major=semi_major,
                    semi_minor=semi_minor,
                )
            )
            for c in coords
        ]
    )
    return xyz

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