Skip to content

API reference

earthcarekit.geo

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

create_spherical_grid 

create_spherical_grid(
    nlat, nlon=None, reduced=False, lat_spacing="regular"
)

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/utils/geo/grid/_create_global_grid.py 
196
197
198
199
200
201
202
203
204
205
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
def create_spherical_grid(
    nlat: int,
    nlon: int | None = None,
    reduced: bool = False,
    lat_spacing: Literal["regular", "sinusoidal", "gaussian"] = "regular",
) -> 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)
    elif lat_spacing == "sinusoidal":
        lat_c, lat_b = _get_sinusoidal_latitudes(nlat)
    elif lat_spacing == "gaussian":
        lat_c, lat_b = _get_gaussian_latitudes(nlat)
    else:
        raise ValueError("grid_type must be 'regular', 'gaussian', or 'sinusoidal'")

    lon_b_list = _get_lon_bounds_per_lat(nlon, lat_c, reduced)

    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,
    y,
    z,
    target_radius=1.0,
    perfect_sphere=True,
    semi_major=SEMI_MAJOR_AXIS_METERS,
    semi_minor=SEMI_MINOR_AXIS_METERS,
)

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

Parameters:

Name Type Description Default
x, y, z float

Cartesian ECEF coordinates.

 required
target_radius float

Target mean radius of the Earth ellipsoid in the new cartesian coordinate system. Defaults to 1.

 1.0
perfect_sphere bool

If True, assume a spherical Earth, else ellipsoidal (WGS-84).

 True
semi_major float

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

 SEMI_MAJOR_AXIS_METERS
semi_minor float

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

 SEMI_MINOR_AXIS_METERS

Returns:

Name Type Description
coords tuple[float, float, float]
  • lat (float): Latitude in degrees
  • lon (float): Longitude in degrees
  • alt (float): Altitude above ellipsoid in meters
Source code in earthcarekit/utils/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
153
154
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)
    back to geodetic coordinates (latitude, longitude, altitude).

    Args:
        x, y, z (float): Cartesian ECEF coordinates.
        target_radius (float): Target mean radius of the Earth ellipsoid in the new cartesian coordinate system. Defaults to 1.
        perfect_sphere (bool): If True, assume a spherical Earth, else ellipsoidal (WGS-84).
        semi_major (float, optional): Semi-major axis of the Earth ellipsoid in meters. Defaults to 6378137 (WGS 84).
        semi_minor (float, optional): Semi-minor axis of the Earth ellipsoid in meters. Defaults to 6356752.314245 (WGS 84).

    Returns:
        coords (tuple[float, float, float]):

            - lat (float): Latitude in degrees
            - lon (float): Longitude in degrees
            - alt (float): Altitude above 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,
    lon,
    alt=None,
    target_radius=1.0,
    perfect_sphere=True,
    semi_major=SEMI_MAJOR_AXIS_METERS,
    semi_minor=SEMI_MINOR_AXIS_METERS,
)

Converts geodetic coordinates (i.e. latitude, longitude and altitude above ellipsoid) to Earth-centered, Earth-fixed (ECEF) coordinates (i.e. x, y and z in cartesian coordinates).

Parameters:

Name Type Description Default
lat float

Latitude angle north (positive) and south (negative) of the equator in degrees.

 required
lon float

Longitude angle east (positive) and west (negative) of the prime meridian in degrees.

 required
alt float

Height above above the Earth ellipsoid in meters.

 None
target_radius float

Target mean radius of the Earth ellipsoid in the new cartesian coordinate system. Defaults to 1.

 1.0
semi_major float

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

 SEMI_MAJOR_AXIS_METERS
semi_minor float

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

 SEMI_MINOR_AXIS_METERS

Returns:

Name Type Description
coords tuple[float, float, float]

3D coordinates in meters (ECEF: A right-handed cartesian coordinate system that has its origin at the Earth's center and is fixed with respect to the Earth's rotation).

  • x (float): Point along the axis passing through the equator at the prime meridian (i.e. latitude = 0, longitude = 0 degrees).
  • y (float): Point along the axis passing through the equator 90 degrees east of the Prime Meridian (i.e. latitude = 0, longitude = 90 degrees).
  • z (float): Point along the axis passing through the north pole (i.e. latitude = 90 degrees).
Source code in earthcarekit/utils/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 (i.e. latitude, longitude and altitude above ellipsoid)
    to Earth-centered, Earth-fixed (ECEF) coordinates (i.e. x, y and z in cartesian coordinates).

    Args:
        lat (float): Latitude angle north (positive) and south (negative) of the equator in degrees.
        lon (float): Longitude angle east (positive) and west (negative) of the prime meridian in degrees.
        alt (float, optional): Height above above the Earth ellipsoid in meters.
        target_radius (float, optional): Target mean radius of the Earth ellipsoid in the new cartesian coordinate system. Defaults to 1.
        semi_major (float, optional): Semi-major axis of the Earth ellipsoid in meters. Defaults to 6378137 (WGS 84).
        semi_minor (float, optional): Semi-minor axis of the Earth ellipsoid in meters. Defaults to 6356752.314245 (WGS 84).

    Returns:
        coords (tuple[float, float, float]): 3D coordinates in meters (ECEF: A right-handed cartesian coordinate system that has its origin at the Earth's center and is fixed with respect to the Earth's rotation).

            - x (float): Point along the axis passing through the equator at the prime meridian (i.e. latitude = 0, longitude = 0 degrees).
            - y (float): Point along the axis passing through the equator 90 degrees east of the Prime Meridian (i.e. latitude = 0, longitude = 90 degrees).
            - z (float): Point along the axis passing through the north pole (i.e. latitude = 90 degrees).
    """
    lat = float(lat)
    lon = float(lon)
    if alt is None:
        alt = 0.0
    else:
        alt = float(alt)

    sin = lambda x: math.sin(x)
    cos = lambda x: math.cos(x)
    sqrt = lambda x: 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, b, units="km", tolerance=1e-12, max_iterations=10
)

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/utils/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
150
151
152
153
154
155
156
157
158
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

get_central_coords 

get_central_coords(latitude, longitude)

Calculates the central lat/lon coordinates.

Source code in earthcarekit/utils/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)

Calculates the central latitude coordinate.

Source code in earthcarekit/utils/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)

Calculates the central longitude coordinate.

Source code in earthcarekit/utils/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, coord2, f=0.5)

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:

Name Type Description
NDArray NDArray

A 2-element numpy.ndarray representing the interpolated lat/lon point.
Source code in earthcarekit/utils/geo/interpolate.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
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 (ArrayLike): The first lat/lon point.
        coord2 (ArrayLike): The second lat/lon point.
        f (float): A fractional value between 0 and 1. Defaults to 0.5, i.e., the mid point between coord1 and coord2.

    Returns:
        NDArray: 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(f"coord1 must be a 2-element sequence (lat, lon)")

    if coord2.shape != (2,):
        raise ValueError(f"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,
    *,
    lat_var=TRACK_LAT_VAR,
    lon_var=TRACK_LON_VAR,
    flatten=False
)

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. Defaults to TRACK_LAT_VAR.

 TRACK_LAT_VAR
lon_var str

Name of the longitude variable. Defaults to TRACK_LON_VAR.

 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

numpy.array: The extracted lat/lon coordinates.
Source code in earthcarekit/utils/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 (str, optional): Name of the latitude variable. Defaults to TRACK_LAT_VAR.
        lon_var (str, optional): Name of the longitude variable. Defaults to TRACK_LON_VAR.
        flatten (bool, optional):
            If True, the coordinates will be flattened to a 2D array

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

    Returns:
        numpy.array: 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, b, units="km", radius_m=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.

 required

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/utils/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_m: 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."
        )

    radius: float = radius_m
    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]

    hav = lambda theta: (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, lon1, lat2, lon2, f)

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/utils/geo/interpolate.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
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

vincenty 

vincenty(
    a, b, units="km", tolerance=1e-12, max_iterations=10
)

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/utils/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
150
151
152
153
154
155
156
157
158
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