"""
Map Projections for Tracking and Navigation.
This module provides map projection functions commonly used in tracking,
navigation, and geospatial applications. All projections convert between
geodetic coordinates (latitude, longitude) and planar map coordinates (x, y).
Projections included:
- Mercator: Conformal cylindrical projection, good for small-scale maps
- Transverse Mercator: Conformal projection for UTM zones
- UTM: Universal Transverse Mercator with zone handling
- Stereographic: Conformal azimuthal projection, good for polar regions
- Lambert Conformal Conic: Conformal conic projection for mid-latitudes
- Azimuthal Equidistant: Preserves distances from center point
All angles are in radians unless otherwise noted.
References
----------
- Snyder, J. P. "Map Projections: A Working Manual." U.S. Geological
Survey Professional Paper 1395, 1987.
- NIMA Technical Report 8350.2, "Department of Defense World Geodetic
System 1984," Third Edition, 2000.
- Karney, C. F. F. "Transverse Mercator with an accuracy of a few
nanometers." Journal of Geodesy 85.8 (2011): 475-485.
"""
from typing import Any, NamedTuple, Optional, Tuple, Union, cast
import numpy as np
from numpy.typing import NDArray
from pytcl.core.array_utils import wrap_to_pi
# WGS84 ellipsoid parameters
WGS84_A = 6378137.0 # Semi-major axis (m)
WGS84_F = 1.0 / 298.257223563 # Flattening
WGS84_B = WGS84_A * (1 - WGS84_F) # Semi-minor axis (m)
WGS84_E = np.sqrt(2 * WGS84_F - WGS84_F**2) # First eccentricity
WGS84_E2 = WGS84_E**2 # First eccentricity squared
WGS84_EP2 = WGS84_E2 / (1 - WGS84_E2) # Second eccentricity squared
[docs]
class ProjectionResult(NamedTuple):
"""Result of a map projection.
Parameters
----------
x : float
Easting coordinate in meters.
y : float
Northing coordinate in meters.
scale : float
Scale factor at the point (ratio of map distance to ground distance).
convergence : float
Grid convergence angle in radians (angle from true north to grid north).
"""
x: float
y: float
scale: float
convergence: float
[docs]
class UTMResult(NamedTuple):
"""Result of UTM projection.
Parameters
----------
easting : float
UTM easting coordinate in meters.
northing : float
UTM northing coordinate in meters.
zone : int
UTM zone number (1-60).
hemisphere : str
'N' for northern hemisphere, 'S' for southern.
scale : float
Scale factor at the point.
convergence : float
Grid convergence angle in radians.
"""
easting: float
northing: float
zone: int
hemisphere: str
scale: float
convergence: float
# =============================================================================
# Mercator Projection
# =============================================================================
def _wrap_longitude_difference(
lon: Union[float, NDArray[np.floating[Any]]],
lon0: Union[float, NDArray[np.floating[Any]]],
) -> Union[float, NDArray[np.floating[Any]]]:
"""Longitude difference ``lon - lon0``, wrapped into [-pi, pi).
Cylindrical and transverse projections compute their easting and
series terms from this difference; without wrapping, a central
meridian near +-180 degrees and a point on the far side of the
antimeridian produce a difference near +-360 degrees instead of near
zero (gh-25 follow-up).
``[()]`` rather than ``float()``: indexing a 0-d array with the empty
tuple returns a ``np.float64`` scalar (a ``float`` subclass) without
forcing array input through ``float()``, which raises for anything but
a 0-d array. An array `lon` therefore still comes back as an array.
"""
return wrap_to_pi(lon - lon0)[()]
[docs]
def mercator(
lat: float,
lon: float,
lon0: float = 0.0,
a: float = WGS84_A,
e: float = WGS84_E,
) -> ProjectionResult:
"""
Ellipsoidal Mercator projection (forward).
The Mercator projection is a conformal cylindrical projection where
rhumb lines (lines of constant bearing) appear as straight lines.
Scale increases with latitude, becoming infinite at the poles.
Parameters
----------
lat : float
Geodetic latitude in radians.
lon : float
Longitude in radians.
lon0 : float, optional
Central meridian in radians. Default is 0.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e : float, optional
First eccentricity. Default is WGS84.
Returns
-------
ProjectionResult
Projected coordinates with scale and convergence.
Notes
-----
The Mercator projection is conformal, meaning it preserves local angles.
However, it distorts areas, especially at high latitudes.
Examples
--------
>>> import numpy as np
>>> result = mercator(np.radians(45), np.radians(-75))
>>> print(f"x={result.x:.1f}, y={result.y:.1f}")
x=-8348961.8, y=5591295.9
"""
# Easting
x = a * _wrap_longitude_difference(lon, lon0)
# Northing using isometric latitude
sin_lat = np.sin(lat)
y = a * np.log(
np.tan(np.pi / 4 + lat / 2) * ((1 - e * sin_lat) / (1 + e * sin_lat)) ** (e / 2)
)
# Scale factor. Uses the caller's eccentricity: this read the global
# WGS84_E2 regardless of the `e` argument, so a custom ellipsoid was
# honored in the projected coordinates and ignored in the scale (gh-25).
cos_lat = np.cos(lat)
w = np.sqrt(1 - e**2 * sin_lat**2)
scale = w / cos_lat if cos_lat > 1e-10 else np.inf
# Convergence is zero for Mercator
convergence = 0.0
return ProjectionResult(x, y, scale, convergence)
[docs]
def mercator_inverse(
x: float,
y: float,
lon0: float = 0.0,
a: float = WGS84_A,
e: float = WGS84_E,
tol: float = 1e-12,
max_iter: int = 10,
) -> Tuple[float, float]:
"""
Ellipsoidal Mercator projection (inverse).
Parameters
----------
x : float
Easting coordinate in meters.
y : float
Northing coordinate in meters.
lon0 : float, optional
Central meridian in radians. Default is 0.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e : float, optional
First eccentricity. Default is WGS84.
tol : float, optional
Convergence tolerance. Default is 1e-12.
max_iter : int, optional
Maximum iterations. Default is 10.
Returns
-------
Tuple[float, float]
(latitude, longitude) in radians.
Examples
--------
>>> import numpy as np
>>> lat, lon = mercator_inverse(1000000, 5000000)
"""
# Longitude, canonicalized into [-pi, pi): lon0 + offset can otherwise
# land just past the antimeridian for a central meridian near +-180
# degrees (gh-25 follow-up). Cast: this function's public signature is
# scalar-only (`Tuple[float, float]`, unchanged here); the array
# support `_wrap_longitude_difference` restores is for other callers
# in this module, so its wider Union return is narrowed back for the
# type checker without touching this function's own declared type.
lon = cast(float, _wrap_longitude_difference(x / a + lon0, 0.0))
# Latitude using iterative solution
t = np.exp(-y / a)
lat = np.pi / 2 - 2 * np.arctan(t)
for _ in range(max_iter):
sin_lat = np.sin(lat)
lat_new = np.pi / 2 - 2 * np.arctan(
t * ((1 - e * sin_lat) / (1 + e * sin_lat)) ** (e / 2)
)
if abs(lat_new - lat) < tol:
break
lat = lat_new
return lat, lon
# =============================================================================
# Transverse Mercator Projection
# =============================================================================
[docs]
def transverse_mercator(
lat: float,
lon: float,
lat0: float = 0.0,
lon0: float = 0.0,
k0: float = 1.0,
a: float = WGS84_A,
e2: float = WGS84_E2,
) -> ProjectionResult:
"""
Transverse Mercator projection (forward).
The Transverse Mercator is a conformal projection that rotates the
cylinder to be tangent along a meridian. It's the basis for UTM.
Parameters
----------
lat : float
Geodetic latitude in radians.
lon : float
Longitude in radians.
lat0 : float, optional
Origin latitude in radians. Default is 0.
lon0 : float, optional
Central meridian in radians. Default is 0.
k0 : float, optional
Scale factor at central meridian. Default is 1.0.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e2 : float, optional
First eccentricity squared. Default is WGS84.
Returns
-------
ProjectionResult
Projected coordinates with scale and convergence.
Notes
-----
This implementation uses the Redfearn series expansion. Measured against
PROJ's ``+proj=tmerc`` on WGS84, worst case over latitudes 0-75 degrees:
================== ==========
From the meridian Difference
================== ==========
2 degrees < 0.1 mm
3 degrees 0.1 mm
4 degrees 0.7 mm
5 degrees 3.3 mm
6 degrees 11.5 mm
================== ==========
The docstring previously said "about 1 meter within 4 degrees", which
understated it by three orders of magnitude and could have led a caller to
reject the function as too coarse (gh-25 follow-up).
Examples
--------
>>> import numpy as np
>>> result = transverse_mercator(np.radians(45), np.radians(-75),
... lon0=np.radians(-75))
"""
# Derived constants
ep2 = e2 / (1 - e2) # Second eccentricity squared
# Third flattening, from the caller's own ellipsoid. The semi-minor axis
# was taken from the global WGS84_B even when a custom `a` and `e2` were
# supplied, which silently mixed two ellipsoids (gh-25). The inverse
# sibling below (`transverse_mercator_inverse`) had the identical bug,
# missed by the original gh-25 fix; it is now fixed the same way.
b = a * np.sqrt(1 - e2)
n = (a - b) / (a + b)
# Compute meridian arc length from equator (Helmert series in
# geodetic latitude: M = A*(lat - a1*sin2lat + a2*sin4lat - a3*sin6lat))
A = a / (1 + n) * (1 + n**2 / 4 + n**4 / 64)
alpha = [
0, # alpha[0] unused
3 * n / 2 - 9 * n**3 / 16,
-(15 * n**2 / 16 - 15 * n**4 / 32),
35 * n**3 / 48,
]
# Conformal latitude
sin_lat = np.sin(lat)
cos_lat = np.cos(lat)
tan_lat = np.tan(lat)
eta2 = ep2 * cos_lat**2
# Longitude difference, wrapped so a central meridian near +-180 degrees
# doesn't see a point on the far side of the antimeridian as ~360
# degrees away (gh-25 follow-up).
dlon = _wrap_longitude_difference(lon, lon0)
# Compute projection using series
t = tan_lat
t2 = t**2
t4 = t**4
c = eta2
c2 = c**2
cos2 = cos_lat**2
cos4 = cos_lat**4
cos6 = cos_lat**6
dl = dlon
dl2 = dl**2
dl4 = dl**4
dl6 = dl**6
# Meridian arc length
N = a / np.sqrt(1 - e2 * sin_lat**2) # Radius of curvature in prime vertical
# Arc length from equator to latitude. `sigma = sigma - ...` (not `-=`)
# so the series is computed out of place: `sigma = lat` aliases the
# caller's array, and an in-place `-=` would mutate it through that
# alias across successive loop iterations (gh-25 follow-up).
sigma = lat
for i in range(1, 4):
sigma = sigma - alpha[i] * np.sin(2 * i * lat)
M = A * sigma
# Arc length from origin latitude (same aliasing hazard as above).
sigma0 = lat0
for i in range(1, 4):
sigma0 = sigma0 - alpha[i] * np.sin(2 * i * lat0)
M0 = A * sigma0
# Easting
x = (
k0
* N
* (
dl * cos_lat
+ dl**3 * cos_lat**3 / 6 * (1 - t2 + c)
+ dl**5 * cos_lat**5 / 120 * (5 - 18 * t2 + t4 + 14 * c - 58 * t2 * c)
)
)
# Northing
y = k0 * (
(M - M0)
+ N * tan_lat * (dl2 * cos2 / 2 + dl4 * cos4 / 24 * (5 - t2 + 9 * c + 4 * c2))
+ N * tan_lat * dl6 * cos6 / 720 * (61 - 58 * t2 + t4 + 270 * c - 330 * t2 * c)
)
# Scale factor
scale = k0 * (
1
+ dl2 * cos2 / 2 * (1 + c)
+ dl4 * cos4 / 24 * (5 - 4 * t2 + 14 * c + 13 * c2 - 28 * t2 * c)
)
# Convergence
convergence = (
dl * sin_lat
+ dl**3 * sin_lat * cos2 / 3 * (1 + 3 * c + 2 * c2)
+ dl**5 * sin_lat * cos4 / 15 * (2 - t2)
)
return ProjectionResult(x, y, scale, convergence)
[docs]
def transverse_mercator_inverse(
x: float,
y: float,
lat0: float = 0.0,
lon0: float = 0.0,
k0: float = 1.0,
a: float = WGS84_A,
e2: float = WGS84_E2,
) -> Tuple[float, float]:
"""
Transverse Mercator projection (inverse).
Parameters
----------
x : float
Easting coordinate in meters.
y : float
Northing coordinate in meters.
lat0 : float, optional
Origin latitude in radians. Default is 0.
lon0 : float, optional
Central meridian in radians. Default is 0.
k0 : float, optional
Scale factor at central meridian. Default is 1.0.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e2 : float, optional
First eccentricity squared. Default is WGS84.
Returns
-------
Tuple[float, float]
(latitude, longitude) in radians.
"""
# Derived constants
ep2 = e2 / (1 - e2)
# Third flattening, from the caller's own ellipsoid (mirrors the forward
# function's gh-25 fix): b must come from the caller's a/e2, not the
# global WGS84_B, or a non-WGS84 caller silently mixes two ellipsoids.
b = a * np.sqrt(1 - e2)
n = (a - b) / (a + b)
A = a / (1 + n) * (1 + n**2 / 4 + n**4 / 64)
# Rectifying -> geodetic (footpoint) latitude series
beta = [
0,
3 * n / 2 - 27 * n**3 / 32,
21 * n**2 / 16 - 55 * n**4 / 32,
151 * n**3 / 96,
]
# Arc length from origin latitude (must match the forward series).
# `sigma0 = sigma0 - ...` (not `-=`): `sigma0 = lat0` aliases the
# caller's array, and an in-place `-=` would mutate it through that
# alias (gh-25 follow-up).
sigma0 = lat0
alpha = [
0,
3 * n / 2 - 9 * n**3 / 16,
-(15 * n**2 / 16 - 15 * n**4 / 32),
35 * n**3 / 48,
]
for i in range(1, 4):
sigma0 = sigma0 - alpha[i] * np.sin(2 * i * lat0)
M0 = A * sigma0
# Footprint latitude
M = M0 + y / k0
mu = M / A
# `lat_fp = lat_fp + ...` (not `+=`): `lat_fp = mu` aliases mu, and the
# loop reads mu on every iteration, so in-place accumulation feeds each
# term a mu already corrupted by the previous one for array input.
lat_fp = mu
for i in range(1, 4):
lat_fp = lat_fp + beta[i] * np.sin(2 * i * mu)
# Parameters at footprint latitude
sin_fp = np.sin(lat_fp)
cos_fp = np.cos(lat_fp)
tan_fp = np.tan(lat_fp)
N_fp = a / np.sqrt(1 - e2 * sin_fp**2)
R_fp = a * (1 - e2) / (1 - e2 * sin_fp**2) ** 1.5
t = tan_fp
t2 = t**2
t4 = t**4
c = ep2 * cos_fp**2
c2 = c**2
d = x / (k0 * N_fp)
d2 = d**2
d4 = d**4
d6 = d**6
# Latitude
lat = lat_fp - N_fp * tan_fp / R_fp * (
d2 / 2
- d4 / 24 * (5 + 3 * t2 + 10 * c - 4 * c2 - 9 * ep2)
+ d6 / 720 * (61 + 90 * t2 + 298 * c + 45 * t4 - 252 * ep2 - 3 * c2)
)
# Longitude, canonicalized into [-pi, pi): for a central meridian near
# +-180 degrees (a forced UTM zone across the antimeridian), lon0 +
# offset can otherwise land just past the antimeridian (gh-25
# follow-up). Cast: narrows the helper's wider Union return back to
# this function's own unchanged scalar-only declared type.
lon = cast(
float,
_wrap_longitude_difference(
lon0
+ (
d
- d**3 / 6 * (1 + 2 * t2 + c)
+ d**5 / 120 * (5 - 2 * c + 28 * t2 - 3 * c2 + 8 * ep2 + 24 * t4)
)
/ cos_fp,
0.0,
),
)
return lat, lon
# =============================================================================
# UTM Projection
# =============================================================================
[docs]
def utm_zone(lon: float, lat: float = 0.0) -> int:
"""
Determine UTM zone number from longitude.
Parameters
----------
lon : float
Longitude in radians.
lat : float, optional
Latitude in radians (used for Norway/Svalbard exceptions).
Returns
-------
int
UTM zone number (1-60).
Notes
-----
Standard zones are 6 degrees wide. Special zones exist for Norway
and Svalbard.
"""
lon_deg = np.degrees(lon)
lat_deg = np.degrees(lat)
# Normalize longitude to -180 to 180
while lon_deg < -180:
lon_deg += 360
while lon_deg >= 180:
lon_deg -= 360
# Standard zone calculation
zone = int((lon_deg + 180) / 6) + 1
# Norway exception (zone 32 extended)
if 56 <= lat_deg < 64 and 3 <= lon_deg < 12:
zone = 32
# Svalbard exceptions
if 72 <= lat_deg < 84:
if 0 <= lon_deg < 9:
zone = 31
elif 9 <= lon_deg < 21:
zone = 33
elif 21 <= lon_deg < 33:
zone = 35
elif 33 <= lon_deg < 42:
zone = 37
return zone
[docs]
def utm_central_meridian(zone: int) -> float:
"""
Get central meridian for UTM zone.
Parameters
----------
zone : int
UTM zone number (1-60).
Returns
-------
float
Central meridian in radians.
"""
return np.radians((zone - 1) * 6 - 180 + 3)
[docs]
def geodetic2utm(lat: float, lon: float, zone: Optional[int] = None) -> UTMResult:
"""
Convert geodetic coordinates to UTM.
Parameters
----------
lat : float
Geodetic latitude in radians.
lon : float
Longitude in radians.
zone : int, optional
Force specific UTM zone. If None, computed from longitude.
Returns
-------
UTMResult
UTM coordinates with zone information.
Examples
--------
>>> import numpy as np
>>> result = geodetic2utm(np.radians(45.0), np.radians(-75.5))
>>> print(f"Zone {result.zone}{result.hemisphere}: "
... f"E={result.easting:.1f}, N={result.northing:.1f}")
Zone 18N: E=460592.4, N=4983072.0
"""
# Determine zone
if zone is None:
zone = utm_zone(lon, lat)
# Central meridian
lon0 = utm_central_meridian(zone)
# UTM scale factor
k0 = 0.9996
# Project
result = transverse_mercator(lat, lon, lat0=0.0, lon0=lon0, k0=k0)
# Apply false easting/northing
easting = result.x + 500000.0
if lat >= 0:
northing = result.y
hemisphere = "N"
else:
northing = result.y + 10000000.0
hemisphere = "S"
return UTMResult(
easting, northing, zone, hemisphere, result.scale, result.convergence
)
[docs]
def utm2geodetic(
easting: float,
northing: float,
zone: int,
hemisphere: str = "N",
) -> Tuple[float, float]:
"""
Convert UTM coordinates to geodetic.
Parameters
----------
easting : float
UTM easting in meters.
northing : float
UTM northing in meters.
zone : int
UTM zone number (1-60).
hemisphere : str, optional
'N' for northern, 'S' for southern hemisphere.
Returns
-------
Tuple[float, float]
(latitude, longitude) in radians.
Examples
--------
>>> lat, lon = utm2geodetic(500000, 5000000, 18, 'N')
>>> print(f"Lat: {np.degrees(lat):.4f}, Lon: {np.degrees(lon):.4f}")
Lat: 45.1535, Lon: -75.0000
"""
# Remove false easting/northing
x = easting - 500000.0
if hemisphere.upper() == "S":
y = northing - 10000000.0
else:
y = northing
# Central meridian
lon0 = utm_central_meridian(zone)
# UTM scale factor
k0 = 0.9996
return transverse_mercator_inverse(x, y, lat0=0.0, lon0=lon0, k0=k0)
# =============================================================================
# Stereographic Projection
# =============================================================================
[docs]
def stereographic(
lat: float,
lon: float,
lat0: float,
lon0: float,
k0: float = 1.0,
a: float = WGS84_A,
e2: float = WGS84_E2,
) -> ProjectionResult:
"""
Oblique stereographic projection (forward).
The stereographic projection is conformal and azimuthal. It's commonly
used for polar regions and local surveys.
Parameters
----------
lat : float
Geodetic latitude in radians.
lon : float
Longitude in radians.
lat0 : float
Center latitude in radians.
lon0 : float
Center longitude in radians.
k0 : float, optional
Scale factor at center. Default is 1.0.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e2 : float, optional
First eccentricity squared. Default is WGS84.
Returns
-------
ProjectionResult
Projected coordinates with scale and convergence.
Notes
-----
Do not use ``lat0 = +-pi/2`` for polar work: this function's conformal
latitude substitution diverges from PROJ's polar stereographic by 5.7 km
at 85 degrees latitude and 34.7 km at 60 degrees (measured against
``+proj=stere`` on WGS84). Use :func:`polar_stereographic` instead,
which matches PROJ/UPS to sub-nanometer precision. For the ellipsoidal
oblique case (any other ``lat0``), use :func:`oblique_stereographic`
(EPSG method 9809), which matches PROJ's ``+proj=sterea`` to
sub-micrometer precision; this function's own oblique accuracy is
documented below and is generally worse.
Examples
--------
>>> import numpy as np
>>> # Polar stereographic centered at North Pole
>>> result = stereographic(np.radians(85), np.radians(45),
... np.radians(90), 0)
Notes
-----
The oblique case maps onto a conformal sphere of *Gaussian* radius,
``sqrt(M0 * N0)`` at the origin latitude. PROJ's ``+proj=sterea`` uses a
sphere of radius ``N0`` instead, so the two disagree away from the center
(gh-25). **How much depends strongly on the origin latitude**, because the
two radii coincide where the meridional and prime-vertical curvatures do.
Measured against ``+proj=sterea`` on WGS84:
========== ============= ==============
``lat0`` At 400 km At 1,738 km
========== ============= ==============
0 deg 1.34 km 5.9 km
30 deg 0.03 km 0.73 km
45 deg 1.38 km 6.7 km
60 deg 2.73 km 12.5 km
75 deg 3.71 km 16.5 km
========== ============= ==============
An earlier version of this note quoted a single pair of figures, 2.5 km and
8.7 km, as though they were properties of the projection. They are one
origin latitude: the disagreement varies eightyfold across the table above.
Both are conformal and each is self-consistent with its own inverse. But
coordinates from this function are not interchangeable with ``+proj=sterea``
output, and it is not survey-grade far from the projection center. For
coordinates that match PROJ (to sub-micrometer), use
:func:`oblique_stereographic` -- the EPSG method 9809 implementation.
"""
e = np.sqrt(e2)
# Conformal latitudes
def conformal_lat(phi: float) -> float:
sin_phi = np.sin(phi)
return (
2
* np.arctan(
np.tan(np.pi / 4 + phi / 2)
* ((1 - e * sin_phi) / (1 + e * sin_phi)) ** (e / 2)
)
- np.pi / 2
)
chi = conformal_lat(lat)
chi0 = conformal_lat(lat0)
# Radius of curvature
sin_lat0 = np.sin(lat0)
N0 = a / np.sqrt(1 - e2 * sin_lat0**2)
R = N0 * np.sqrt(1 - e2) / (1 - e2 * sin_lat0**2)
# Projection
sin_chi = np.sin(chi)
cos_chi = np.cos(chi)
sin_chi0 = np.sin(chi0)
cos_chi0 = np.cos(chi0)
dlon = lon - lon0
cos_dlon = np.cos(dlon)
sin_dlon = np.sin(dlon)
k_denom = 1 + sin_chi0 * sin_chi + cos_chi0 * cos_chi * cos_dlon
k = 2 * k0 / k_denom
x = k * R * cos_chi * sin_dlon
y = k * R * (cos_chi0 * sin_chi - sin_chi0 * cos_chi * cos_dlon)
# Scale factor
scale = k
# Grid convergence: the angle from true north to grid north. The
# stereographic projection is conformal, so this is exactly the
# difference between the grid and true bearings of any one direction;
# taking the direction from the point back to the projection centre,
# whose grid bearing is ``azimuth + pi`` because that ray is straight in
# the map, gives ``back_azimuth - azimuth - pi``. Both azimuths are on
# the conformal sphere, matching the projection above.
azimuth = np.arctan2(
cos_chi * sin_dlon,
cos_chi0 * sin_chi - sin_chi0 * cos_chi * cos_dlon,
)
back_azimuth = np.arctan2(
-cos_chi0 * sin_dlon,
cos_chi * sin_chi0 - sin_chi * cos_chi0 * cos_dlon,
)
convergence = float(wrap_to_pi(back_azimuth - azimuth - np.pi))
return ProjectionResult(x, y, scale, convergence)
[docs]
def stereographic_inverse(
x: float,
y: float,
lat0: float,
lon0: float,
k0: float = 1.0,
a: float = WGS84_A,
e2: float = WGS84_E2,
tol: float = 1e-12,
max_iter: int = 10,
) -> Tuple[float, float]:
"""
Oblique stereographic projection (inverse).
Parameters
----------
x : float
Easting in meters.
y : float
Northing in meters.
lat0 : float
Center latitude in radians.
lon0 : float
Center longitude in radians.
k0 : float, optional
Scale factor at center. Default is 1.0.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e2 : float, optional
First eccentricity squared. Default is WGS84.
tol : float, optional
Convergence tolerance. Default is 1e-12.
max_iter : int, optional
Maximum iterations. Default is 10.
Returns
-------
Tuple[float, float]
(latitude, longitude) in radians.
"""
e = np.sqrt(e2)
# Conformal latitude of center
sin_lat0 = np.sin(lat0)
chi0 = (
2
* np.arctan(
np.tan(np.pi / 4 + lat0 / 2)
* ((1 - e * sin_lat0) / (1 + e * sin_lat0)) ** (e / 2)
)
- np.pi / 2
)
sin_chi0 = np.sin(chi0)
cos_chi0 = np.cos(chi0)
# Radius
N0 = a / np.sqrt(1 - e2 * sin_lat0**2)
R = N0 * np.sqrt(1 - e2) / (1 - e2 * sin_lat0**2)
# Distance from center
rho = np.sqrt(x**2 + y**2)
if rho < 1e-10:
return lat0, lon0
c = 2 * np.arctan(rho / (2 * k0 * R))
sin_c = np.sin(c)
cos_c = np.cos(c)
# Conformal latitude
chi = np.arcsin(cos_c * sin_chi0 + y * sin_c * cos_chi0 / rho)
# Longitude, canonicalized into [-pi, pi): lon0 + offset can otherwise
# land just past the antimeridian for a centre near +-180 degrees
# (gh-25 follow-up). Cast: narrows the helper's wider Union return
# back to this function's own unchanged scalar-only declared type.
lon = cast(
float,
_wrap_longitude_difference(
lon0 + np.arctan2(x * sin_c, rho * cos_chi0 * cos_c - y * sin_chi0 * sin_c),
0.0,
),
)
# Invert conformal latitude
lat = chi
for _ in range(max_iter):
sin_lat = np.sin(lat)
lat_new = (
2
* np.arctan(
np.tan(np.pi / 4 + chi / 2)
* ((1 + e * sin_lat) / (1 - e * sin_lat)) ** (e / 2)
)
- np.pi / 2
)
if abs(lat_new - lat) < tol:
break
lat = lat_new
return lat, lon
[docs]
def polar_stereographic(
lat: float,
lon: float,
north: bool = True,
k0: float = 0.994,
a: float = WGS84_A,
e2: float = WGS84_E2,
) -> ProjectionResult:
"""
Polar stereographic projection (forward).
Standard polar stereographic used for polar regions.
Parameters
----------
lat : float
Geodetic latitude in radians.
lon : float
Longitude in radians.
north : bool, optional
True for North Pole center, False for South Pole. Default is True.
k0 : float, optional
Scale factor at pole. Default is 0.994 (UPS standard).
a : float, optional
Semi-major axis in meters. Default is WGS84.
e2 : float, optional
First eccentricity squared. Default is WGS84.
Returns
-------
ProjectionResult
Projected coordinates with scale and convergence.
Notes
-----
Uses the exact ellipsoidal polar stereographic equations
(Snyder 1987, eqs. 15-9 and 21-33/21-34), matching standard
UPS/EPSG implementations.
Examples
--------
>>> import numpy as np
>>> # Arctic location
>>> result = polar_stereographic(np.radians(80), np.radians(45))
"""
e = np.sqrt(e2)
# Work in the north-polar aspect; mirror for the south
if north:
phi, lam = lat, lon
else:
phi, lam = -lat, -lon
sin_phi = np.sin(phi)
t = np.tan(np.pi / 4 - phi / 2) * ((1 + e * sin_phi) / (1 - e * sin_phi)) ** (e / 2)
rho = 2 * a * k0 * t / np.sqrt((1 + e) ** (1 + e) * (1 - e) ** (1 - e))
x = rho * np.sin(lam)
y = -rho * np.cos(lam)
if not north:
x, y = -x, -y
# Scale factor k = rho / (a * m); at the pole the limit is k0
cos_phi = np.cos(phi)
if cos_phi > 1e-12:
m = cos_phi / np.sqrt(1 - e2 * sin_phi**2)
scale = rho / (a * m)
else:
scale = k0
# Grid convergence: +dlon (north aspect), -dlon (south aspect)
convergence = lam
return ProjectionResult(x, y, scale, convergence)
# =============================================================================
# Lambert Conformal Conic Projection
# =============================================================================
# =============================================================================
# Azimuthal Equidistant Projection
# =============================================================================
def _authalic_radius(a: float, e2: float) -> float:
"""Radius of the sphere with the same surface area as the ellipsoid."""
if e2 == 0:
return a
e = np.sqrt(e2)
q_p = 1 + (1 - e2) / (2 * e) * np.log((1 + e) / (1 - e))
return a * np.sqrt(q_p / 2)
[docs]
def azimuthal_equidistant(
lat: float,
lon: float,
lat0: float,
lon0: float,
a: float = WGS84_A,
e2: float = WGS84_E2,
) -> ProjectionResult:
"""
Azimuthal equidistant projection (forward).
Distances from the center point are preserved. Useful for showing
distances from a specific location (e.g., radio coverage).
Parameters
----------
lat : float
Geodetic latitude in radians.
lon : float
Longitude in radians.
lat0 : float
Center latitude in radians.
lon0 : float
Center longitude in radians.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e2 : float, optional
First eccentricity squared. Default is WGS84.
Returns
-------
ProjectionResult
Projected coordinates with scale and convergence.
Notes
-----
This implementation uses the spherical approximation for simplicity.
For exact ellipsoidal distances (matching PROJ's ``+proj=aeqd``), use
:func:`azimuthal_equidistant_exact`.
This projection is not conformal, so there is no single scale factor.
The radial scale is exactly 1 by construction and is what
``ProjectionResult.scale`` reports; the tangential scale is
``c / sin(c)`` for angular distance ``c`` from the centre, reaching
about 1.005 at 1,000 km and 1.11 at 5,000 km. ``convergence`` accounts
for this and is the true grid convergence, not the conformal
approximation.
Examples
--------
>>> import numpy as np
>>> result = azimuthal_equidistant(np.radians(40), np.radians(-75),
... np.radians(38.9), np.radians(-77))
"""
# Use spherical approximation with authalic radius
R = _authalic_radius(a, e2)
sin_lat = np.sin(lat)
cos_lat = np.cos(lat)
sin_lat0 = np.sin(lat0)
cos_lat0 = np.cos(lat0)
dlon = lon - lon0
cos_dlon = np.cos(dlon)
sin_dlon = np.sin(dlon)
# Angular distance from center
cos_c = sin_lat0 * sin_lat + cos_lat0 * cos_lat * cos_dlon
cos_c = np.clip(cos_c, -1, 1)
c = np.arccos(cos_c)
if c < 1e-10:
return ProjectionResult(0.0, 0.0, 1.0, 0.0)
# Azimuth from center
k = R * c / np.sin(c)
x = k * cos_lat * sin_dlon
y = k * (cos_lat0 * sin_lat - sin_lat0 * cos_lat * cos_dlon)
# Radial scale is exactly 1 by construction; this projection is not
# conformal, so there is no single scale factor -- the tangential scale
# is c / sin(c), which grows without bound away from the centre (about
# 1.11 at 5,000 km). ``ProjectionResult.scale`` reports the radial one.
scale = 1.0
# Grid convergence: the angle from true north to grid north. Unlike the
# conformal projections above, the projected meridian here is not simply
# rotated -- decomposed in the map's polar frame at this point it has a
# radial component -cos(back_azimuth) and a tangential component
# (c / sin c) * sin(back_azimuth), and the convergence follows from their
# ratio. The two agree in the small-c limit, where c / sin c -> 1.
azimuth = np.arctan2(
cos_lat * sin_dlon,
cos_lat0 * sin_lat - sin_lat0 * cos_lat * cos_dlon,
)
back_azimuth = np.arctan2(
-cos_lat0 * sin_dlon,
cos_lat * sin_lat0 - sin_lat * cos_lat0 * cos_dlon,
)
convergence = float(
wrap_to_pi(
-azimuth
- np.arctan2(
(c / np.sin(c)) * np.sin(back_azimuth),
-np.cos(back_azimuth),
)
)
)
return ProjectionResult(x, y, scale, convergence)
[docs]
def azimuthal_equidistant_inverse(
x: float,
y: float,
lat0: float,
lon0: float,
a: float = WGS84_A,
e2: float = WGS84_E2,
) -> Tuple[float, float]:
"""
Azimuthal equidistant projection (inverse).
Parameters
----------
x : float
Easting in meters.
y : float
Northing in meters.
lat0 : float
Center latitude in radians.
lon0 : float
Center longitude in radians.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e2 : float, optional
First eccentricity squared. Default is WGS84.
Returns
-------
Tuple[float, float]
(latitude, longitude) in radians.
"""
R = _authalic_radius(a, e2)
rho = np.sqrt(x**2 + y**2)
if rho < 1e-10:
return lat0, lon0
c = rho / R
sin_c = np.sin(c)
cos_c = np.cos(c)
sin_lat0 = np.sin(lat0)
cos_lat0 = np.cos(lat0)
lat = np.arcsin(cos_c * sin_lat0 + y * sin_c * cos_lat0 / rho)
# Canonicalize into [-pi, pi): lon0 + offset can otherwise land just
# past the antimeridian for a centre near +-180 degrees (gh-25
# follow-up). Cast: narrows the helper's wider Union return back to
# this function's own unchanged scalar-only declared type.
lon = cast(
float,
_wrap_longitude_difference(
lon0 + np.arctan2(x * sin_c, rho * cos_lat0 * cos_c - y * sin_lat0 * sin_c),
0.0,
),
)
return lat, lon
def _sterea_constants(
lat0: float, k0: float, a: float, e2: float
) -> Tuple[float, float, float, float, float]:
"""Conformal-sphere constants of EPSG method 9809 at the origin.
Returns (R, n, c, chi0, e) where R is the conformal sphere radius,
n and c the latitude-mapping constants, and chi0 the conformal
latitude of the origin on that sphere.
"""
e = np.sqrt(e2)
sin_lat0 = np.sin(lat0)
rho0 = a * (1 - e2) / (1 - e2 * sin_lat0**2) ** 1.5
nu0 = a / np.sqrt(1 - e2 * sin_lat0**2)
R = np.sqrt(rho0 * nu0)
n = np.sqrt(1 + e2 * np.cos(lat0) ** 4 / (1 - e2))
s1 = (1 + sin_lat0) / (1 - sin_lat0)
s2 = (1 - e * sin_lat0) / (1 + e * sin_lat0)
w1 = (s1 * s2**e) ** n
sin_chi0_prov = (w1 - 1) / (w1 + 1)
c = (n + sin_lat0) * (1 - sin_chi0_prov) / ((n - sin_lat0) * (1 + sin_chi0_prov))
w2 = c * w1
chi0 = np.arcsin((w2 - 1) / (w2 + 1))
return R, n, c, chi0, e
[docs]
def oblique_stereographic(
lat: float,
lon: float,
lat0: float,
lon0: float,
k0: float = 1.0,
a: float = WGS84_A,
e2: float = WGS84_E2,
) -> ProjectionResult:
"""
Oblique (double) stereographic projection, EPSG method 9809 (forward).
The ellipsoid is first mapped conformally onto a sphere of Gaussian
radius at the origin using the full EPSG latitude mapping (constants
``n`` and ``c``), then that sphere is projected stereographically.
This is the algorithm PROJ implements as ``+proj=sterea`` (the
projection of the Dutch RD grid, EPSG:28992), so unlike
:func:`stereographic` -- whose simpler conformal-latitude substitution
diverges from ``sterea`` by kilometers far from the origin (gh-25) --
coordinates from this function are interchangeable with PROJ output.
Parameters
----------
lat : float
Geodetic latitude in radians.
lon : float
Longitude in radians.
lat0 : float
Origin latitude in radians.
lon0 : float
Origin longitude in radians.
k0 : float, optional
Scale factor at the origin. Default is 1.0.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e2 : float, optional
First eccentricity squared. Default is WGS84.
Returns
-------
ProjectionResult
Projected coordinates with scale and convergence.
References
----------
- IOGP Publication 373-7-2, "Coordinate Conversions and
Transformations including Formulas", EPSG method code 9809.
Examples
--------
>>> import numpy as np
>>> res = oblique_stereographic(np.radians(53.0), np.radians(6.0),
... np.radians(52.156161), np.radians(5.387639),
... k0=0.9999079)
>>> bool(abs(res.scale - 0.9999079) < 1e-3)
True
"""
R, n, c, chi0, e = _sterea_constants(lat0, k0, a, e2)
sin_lat = np.sin(lat)
sa = (1 + sin_lat) / (1 - sin_lat)
sb = (1 - e * sin_lat) / (1 + e * sin_lat)
w = c * (sa * sb**e) ** n
chi = np.arcsin((w - 1) / (w + 1))
# Wrap the longitude difference into [-pi, pi) before scaling by n: for
# n != 1 (any real ellipsoid), sin/cos of n * (lon - lon0) is not
# 2*pi-periodic in the unwrapped difference, so an origin near +-180
# degrees and a point on the far side of the antimeridian produced a
# scaled angle far from the true small-angle value (gh-25 follow-up).
dlam = n * _wrap_longitude_difference(lon, lon0)
sin_chi = np.sin(chi)
cos_chi = np.cos(chi)
sin_chi0 = np.sin(chi0)
cos_chi0 = np.cos(chi0)
cos_dlam = np.cos(dlam)
sin_dlam = np.sin(dlam)
B = 1 + sin_chi * sin_chi0 + cos_chi * cos_chi0 * cos_dlam
x = 2 * R * k0 * cos_chi * sin_dlam / B
y = 2 * R * k0 * (sin_chi * cos_chi0 - cos_chi * sin_chi0 * cos_dlam) / B
# Composite conformal scale: ellipsoid -> conformal sphere
# (n R cos(chi) / (nu cos(lat))) times sphere -> plane (2 k0 / B).
nu = a / np.sqrt(1 - e2 * sin_lat**2)
scale = float(2 * n * R * k0 * cos_chi / (B * nu * np.cos(lat)))
# Grid convergence via the ray back to the projection centre, which is
# straight in the map; conformality transfers sphere bearings to
# ellipsoid bearings unchanged (meridians map to meridians).
azimuth = np.arctan2(
cos_chi * sin_dlam,
cos_chi0 * sin_chi - sin_chi0 * cos_chi * cos_dlam,
)
back_azimuth = np.arctan2(
-cos_chi0 * sin_dlam,
cos_chi * sin_chi0 - sin_chi * cos_chi0 * cos_dlam,
)
convergence = float(wrap_to_pi(back_azimuth - azimuth - np.pi))
return ProjectionResult(float(x), float(y), scale, convergence)
[docs]
def oblique_stereographic_inverse(
x: float,
y: float,
lat0: float,
lon0: float,
k0: float = 1.0,
a: float = WGS84_A,
e2: float = WGS84_E2,
) -> Tuple[float, float]:
"""
Oblique (double) stereographic projection, EPSG method 9809 (inverse).
Parameters
----------
x : float
Easting in meters.
y : float
Northing in meters.
lat0 : float
Origin latitude in radians.
lon0 : float
Origin longitude in radians.
k0 : float, optional
Scale factor at the origin. Default is 1.0.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e2 : float, optional
First eccentricity squared. Default is WGS84.
Returns
-------
Tuple[float, float]
(latitude, longitude) in radians.
Examples
--------
>>> import numpy as np
>>> lat0, lon0 = np.radians(52.156161), np.radians(5.387639)
>>> fwd = oblique_stereographic(np.radians(53.0), np.radians(6.0),
... lat0, lon0, k0=0.9999079)
>>> lat, lon = oblique_stereographic_inverse(fwd.x, fwd.y, lat0, lon0,
... k0=0.9999079)
>>> bool(abs(np.degrees(lat) - 53.0) < 1e-9)
True
"""
R, n, c, chi0, e = _sterea_constants(lat0, k0, a, e2)
g = 2 * R * k0 * np.tan(np.pi / 4 - chi0 / 2)
h = 4 * R * k0 * np.tan(chi0) + g
i = np.arctan2(x, h + y)
j = np.arctan2(x, g - y) - i
chi = chi0 + 2 * np.arctan((y - x * np.tan(j / 2)) / (2 * R * k0))
dlam = j + 2 * i
# Canonicalize into [-pi, pi): lon0 + offset can otherwise land just
# past the antimeridian for an origin near +-180 degrees (gh-25
# follow-up).
lon = _wrap_longitude_difference(dlam / n + lon0, 0.0)
# Invert the conformal-latitude mapping for the geodetic latitude:
# psi is the sphere's isometric latitude pulled back through n and c.
sin_chi = np.sin(chi)
psi = 0.5 * np.log((1 + sin_chi) / (c * (1 - sin_chi))) / n
lat = 2 * np.arctan(np.exp(psi)) - np.pi / 2
for _ in range(25):
sin_lat = np.sin(lat)
psi_i = np.log(
np.tan(lat / 2 + np.pi / 4)
* ((1 - e * sin_lat) / (1 + e * sin_lat)) ** (e / 2)
)
lat_new = lat - (psi_i - psi) * np.cos(lat) * (1 - e2 * sin_lat**2) / (1 - e2)
if abs(lat_new - lat) < 1e-14:
lat = lat_new
break
lat = lat_new
return float(lat), float(lon)
[docs]
def azimuthal_equidistant_exact(
lat: float,
lon: float,
lat0: float,
lon0: float,
a: float = WGS84_A,
e2: float = WGS84_E2,
) -> ProjectionResult:
"""
Exact ellipsoidal azimuthal equidistant projection (forward).
The map coordinates are the true geodesic distance and azimuth from
the projection centre: ``x = s12 sin(az1)``, ``y = s12 cos(az1)``,
computed with Karney's algorithm (geographiclib). This is what PROJ's
``+proj=aeqd`` computes on an ellipsoid, so unlike
:func:`azimuthal_equidistant` -- a spherical approximation on the
authalic sphere -- distances from the centre are exact everywhere,
antipodes included.
Requires the ``geodesy`` extra (geographiclib).
Parameters
----------
lat : float
Geodetic latitude in radians.
lon : float
Longitude in radians.
lat0 : float
Centre latitude in radians.
lon0 : float
Centre longitude in radians.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e2 : float, optional
First eccentricity squared. Default is WGS84.
Returns
-------
ProjectionResult
Projected coordinates. ``scale`` is the radial scale, exactly 1
by construction; ``convergence`` is the true grid convergence,
computed from the geodesic azimuths and reduced length.
Raises
------
DependencyError
If geographiclib is not installed.
Examples
--------
>>> import numpy as np
>>> res = azimuthal_equidistant_exact(np.radians(40), np.radians(-75),
... np.radians(38.9), np.radians(-77))
>>> bool(np.hypot(res.x, res.y) > 0)
True
>>> res.scale
1.0
"""
geod = _karney_geodesic(a, e2)
from geographiclib.geodesic import Geodesic
g = geod.Inverse(
np.degrees(lat0),
np.degrees(lon0),
np.degrees(lat),
np.degrees(lon),
Geodesic.STANDARD | Geodesic.REDUCEDLENGTH,
)
s12 = g["s12"]
if s12 < 1e-9:
return ProjectionResult(0.0, 0.0, 1.0, 0.0)
az1 = np.radians(g["azi1"])
az2 = np.radians(g["azi2"])
m12 = g["m12"]
x = s12 * np.sin(az1)
y = s12 * np.cos(az1)
# Radial scale is exactly 1 by construction. The tangential scale is
# s12 / m12 (map arc over the geodesic's reduced length), which the
# convergence must account for, exactly as the spherical version
# accounts for c / sin(c) -- its special case with m12 = R sin(s/R).
back_azimuth = wrap_to_pi(az2 + np.pi)
convergence = float(
wrap_to_pi(
-az1 - np.arctan2((s12 / m12) * np.sin(back_azimuth), -np.cos(back_azimuth))
)
)
return ProjectionResult(float(x), float(y), 1.0, convergence)
[docs]
def azimuthal_equidistant_exact_inverse(
x: float,
y: float,
lat0: float,
lon0: float,
a: float = WGS84_A,
e2: float = WGS84_E2,
) -> Tuple[float, float]:
"""
Exact ellipsoidal azimuthal equidistant projection (inverse).
Solves the direct geodesic problem from the projection centre with
azimuth ``atan2(x, y)`` and distance ``hypot(x, y)`` (Karney's
algorithm, geographiclib). Requires the ``geodesy`` extra.
Parameters
----------
x : float
Easting in meters.
y : float
Northing in meters.
lat0 : float
Centre latitude in radians.
lon0 : float
Centre longitude in radians.
a : float, optional
Semi-major axis in meters. Default is WGS84.
e2 : float, optional
First eccentricity squared. Default is WGS84.
Returns
-------
Tuple[float, float]
(latitude, longitude) in radians.
Raises
------
DependencyError
If geographiclib is not installed.
Examples
--------
>>> import numpy as np
>>> lat0, lon0 = np.radians(38.9), np.radians(-77)
>>> fwd = azimuthal_equidistant_exact(np.radians(40), np.radians(-75),
... lat0, lon0)
>>> lat, lon = azimuthal_equidistant_exact_inverse(fwd.x, fwd.y,
... lat0, lon0)
>>> bool(abs(np.degrees(lat) - 40.0) < 1e-9)
True
"""
geod = _karney_geodesic(a, e2)
s = float(np.hypot(x, y))
if s < 1e-9:
return float(lat0), float(lon0)
az = np.degrees(np.arctan2(x, y))
g = geod.Direct(np.degrees(lat0), np.degrees(lon0), az, s)
return float(np.radians(g["lat2"])), float(np.radians(g["lon2"]))
def _karney_geodesic(a: float, e2: float):
"""A geographiclib Geodesic for (a, e2), or a loud DependencyError."""
try:
from geographiclib.geodesic import Geodesic
except ImportError:
from pytcl.core.exceptions import DependencyError
raise DependencyError(
"azimuthal_equidistant_exact requires geographiclib. "
"Install the `geodesy` extra: pip install nrl-tracker[geodesy]."
) from None
f = 1.0 - np.sqrt(1.0 - e2)
return Geodesic(a, f)
# =============================================================================
# Batch Operations
# =============================================================================
[docs]
def geodetic2utm_batch(
lats: NDArray[np.floating],
lons: NDArray[np.floating],
zone: Optional[int] = None,
) -> tuple[NDArray[np.floating], NDArray[np.floating], NDArray[np.intp], NDArray[Any]]:
"""
Batch convert geodetic coordinates to UTM.
Parameters
----------
lats : ndarray
Geodetic latitudes in radians.
lons : ndarray
Longitudes in radians.
zone : int, optional
Force specific UTM zone. If None, computed per point.
Returns
-------
eastings : ndarray
UTM eastings in meters.
northings : ndarray
UTM northings in meters.
zones : ndarray
UTM zone numbers.
hemispheres : ndarray
Hemisphere characters ('N' or 'S').
"""
n = len(lats)
eastings = np.zeros(n)
northings = np.zeros(n)
zones = np.zeros(n, dtype=np.int_)
hemispheres = np.empty(n, dtype="U1")
for i in range(n):
result = geodetic2utm(lats[i], lons[i], zone)
eastings[i] = result.easting
northings[i] = result.northing
zones[i] = result.zone
hemispheres[i] = result.hemisphere
return eastings, northings, zones, hemispheres
__all__ = [
# Constants
"WGS84_A",
"WGS84_B",
"WGS84_F",
"WGS84_E",
"WGS84_E2",
"WGS84_EP2",
# Result types
"ProjectionResult",
"UTMResult",
# Mercator
"mercator",
"mercator_inverse",
# Transverse Mercator
"transverse_mercator",
"transverse_mercator_inverse",
# UTM
"utm_zone",
"utm_central_meridian",
"geodetic2utm",
"utm2geodetic",
"geodetic2utm_batch",
# Stereographic
"stereographic",
"stereographic_inverse",
"oblique_stereographic",
"oblique_stereographic_inverse",
"polar_stereographic",
# Lambert Conformal Conic
"lambert_conformal_conic",
"lambert_conformal_conic_inverse",
# Azimuthal Equidistant
"azimuthal_equidistant",
"azimuthal_equidistant_inverse",
"azimuthal_equidistant_exact",
"azimuthal_equidistant_exact_inverse",
]