"""Spherical geometry helpers using SLERP (Spherical Linear intERPolation).
See `Spherical linear interpolation
<https://en.wikipedia.org/wiki/Spherical_linear_interpolation>`_.
These are rougher than pyproj in that they model Earth as a perfect sphere of radius
:data:`pycontrails.physics.constants.radius_earth` rather than a WGS84 ellipsoid but
significantly faster for vectorized numpy operations.
The pycontrails library already includes :func:`pycontrails.physics.geo.haversine` and
:func:`pycontrails.physics.geo.azimuth` functions for computing great circle distances
and azimuths. The functions here extend those utilities with forward projection and
interpolation capabilities.
"""
import numpy as np
import numpy.typing as npt
from pycontrails.physics import constants
[docs]
def spherical_fwd(
lon: npt.NDArray[np.floating],
lat: npt.NDArray[np.floating],
az: npt.NDArray[np.floating],
dist: npt.NDArray[np.floating],
) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]]:
r"""Project from ``(lon, lat)`` along azimuth ``az`` by ``dist`` meters on a sphere.
Equivalent to :meth:`pyproj.Geod.fwd`, which returns ``lon2, lat2, _``.
Parameters
----------
lon : npt.NDArray[np.floating]
Longitude of starting point, [:math:`\deg`].
lat : npt.NDArray[np.floating]
Latitude of starting point, [:math:`\deg`].
az : npt.NDArray[np.floating]
Forward azimuth, [:math:`\deg`].
dist : npt.NDArray[np.floating]
Distance to project, [:math:`m`].
Returns
-------
lon2 : npt.NDArray[np.floating]
Longitude of projected point, [:math:`\deg`].
lat2 : npt.NDArray[np.floating]
Latitude of projected point, [:math:`\deg`].
"""
lonr = np.deg2rad(lon)
latr = np.deg2rad(lat)
azr = np.deg2rad(az)
d = dist / constants.radius_earth # angular distance in radians
sin_d = np.sin(d)
cos_d = np.cos(d)
sin_lat = np.sin(latr)
cos_lat = np.cos(latr)
lat2 = np.arcsin(sin_lat * cos_d + cos_lat * sin_d * np.cos(azr))
lon2 = lonr + np.arctan2(
np.sin(azr) * sin_d * cos_lat,
cos_d - sin_lat * np.sin(lat2),
)
return np.rad2deg(lon2), np.rad2deg(lat2)
[docs]
def gc_interp(
lon1: npt.NDArray[np.floating],
lat1: npt.NDArray[np.floating],
lon2: npt.NDArray[np.floating],
lat2: npt.NDArray[np.floating],
frac: npt.NDArray[np.floating],
) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]]:
r"""Interpolate along great circles via SLERP.
Equivalent to :meth:`pyproj.Geod.fwd` evaluated at ``frac * dist``, but without
the separate :meth:`pyproj.Geod.inv` call to obtain ``az`` and ``dist``.
Parameters
----------
lon1 : npt.NDArray[np.floating]
Longitude of source, [:math:`\deg`].
lat1 : npt.NDArray[np.floating]
Latitude of source, [:math:`\deg`].
lon2 : npt.NDArray[np.floating]
Longitude of destination, [:math:`\deg`].
lat2 : npt.NDArray[np.floating]
Latitude of destination, [:math:`\deg`].
frac : npt.NDArray[np.floating]
Fractional position along the arc, 0 = source, 1 = dest.
Returns
-------
lon : npt.NDArray[np.floating]
Interpolated longitude, [:math:`\deg`].
lat : npt.NDArray[np.floating]
Interpolated latitude, [:math:`\deg`].
"""
lon1r = np.deg2rad(lon1)
lat1r = np.deg2rad(lat1)
lon2r = np.deg2rad(lon2)
lat2r = np.deg2rad(lat2)
cos_lat1 = np.cos(lat1r)
cos_lat2 = np.cos(lat2r)
# Central angle via haversine
dlat = lat2r - lat1r
dlon = lon2r - lon1r
a = np.sin(dlat / 2) ** 2 + cos_lat1 * cos_lat2 * np.sin(dlon / 2) ** 2
sigma = 2.0 * np.arctan2(np.sqrt(a), np.sqrt(1.0 - a))
sin_sigma = np.sin(sigma)
# Guard against zero-length arcs: as sigma -> 0, we have A -> 1 - frac and B -> frac
zero_arc = sin_sigma < 1e-12
safe_sin_sigma = np.where(zero_arc, 1.0, sin_sigma) # avoid division by zero
A = np.where(zero_arc, 1.0 - frac, np.sin((1.0 - frac) * sigma) / safe_sin_sigma)
B = np.where(zero_arc, frac, np.sin(frac * sigma) / safe_sin_sigma)
# Cartesian SLERP
x = A * cos_lat1 * np.cos(lon1r) + B * cos_lat2 * np.cos(lon2r)
y = A * cos_lat1 * np.sin(lon1r) + B * cos_lat2 * np.sin(lon2r)
z = A * np.sin(lat1r) + B * np.sin(lat2r)
lat_out = np.rad2deg(np.arctan2(z, np.sqrt(x**2 + y**2)))
lon_out = np.rad2deg(np.arctan2(y, x))
return lon_out, lat_out
[docs]
def gc_npts(
lon1: float, lat1: float, lon2: float, lat2: float, n: int
) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
r"""Return n equally-spaced intermediate points along a great circle.
Excludes the endpoints themselves. Equivalent to :meth:`pyproj.Geod.npts`.
Parameters
----------
lon1 : float
Longitude of source, [:math:`\deg`].
lat1 : float
Latitude of source, [:math:`\deg`].
lon2 : float
Longitude of destination, [:math:`\deg`].
lat2 : float
Latitude of destination, [:math:`\deg`].
n : int
Number of intermediate points.
Returns
-------
lon : npt.NDArray[np.float64]
Longitude of intermediate points, [:math:`\deg`]. Always ``np.float64``.
lat : npt.NDArray[np.float64]
Latitude of intermediate points, [:math:`\deg`]. Always ``np.float64``.
"""
frac = np.linspace(0.0, 1.0, n + 2, dtype=np.float64)[1:-1]
return gc_interp(lon1, lat1, lon2, lat2, frac)