PyGeodetics is a pure-Python library for practical geodesy. It bundles coordinate conversions, map projections, geodesic problem solvers and local-frame utilities in a single, dependency-light package that works equally well with scalar inputs and NumPy arrays.
What's supported
- Global coordinate conversions — geodetic ↔ ECEF with three inverse solvers (Bowring, Vermeille, Heiskanen–Moritz iterative).
- Local tangent frames — full round-trip ENU / NED / ECEF / geodetic, plus Azimuth / Elevation / Range (AER) look-angle conversions.
- Geodesic problems on the ellipsoid — Vincenty forward (direct) and inverse problems, Vincenty geodesic distance, meridional (M) and prime-vertical (N) radii of curvature, radius of curvature for a given azimuth (Euler), mean radius.
- Map projections — Transverse Mercator (EPSG 9807), Mercator Variant C (EPSG 1044), Polar Stereographic (EPSG 9810), with the grid convergence and point-scale factor computed from either geographic or projected coordinates.
- Grid systems — turnkey UTM (
geodetic_to_utm/utm_to_geodeticwith automatic zone and hemisphere handling) and UPS (geodetic_to_ups/ups_to_geodetic) for the polar regions above 84°N / below 80°S. - MGRS — encoding and parsing across UTM and UPS at every standard precision level (100 km → 1 m).
- EPSG helpers — compute the EPSG code for any WGS84 / NAD83 UTM zone
or WGS84 UPS polar CRS from a simple arithmetic formula (no external
database or network call):
utm_epsg,ups_epsg. - Polygon utilities on the ellipsoid — perimeter, area (with an
optional
geodesic=Truemode for accurate continental / hemispheric polygons), centroid, bounding box, geodesic densification (insert vertices so no edge exceeds a given length) and geodesic interpolation between two points. - Multiple reference ellipsoids — WGS84 (default), GRS80,
International 1924, Clarke 1866, Bessel 1841 and Bessel Modified; any
computation accepts a custom
Ellipsoidinstance. - NumPy-friendly — every public function accepts scalars or arrays and returns the same shape, making the library usable as a drop-in for batch processing and notebooks.
pip install pygeodeticsCoordinate system conversions
Local tangent frames (ENU / NED / AER)
Geodesic problems on the ellipsoid
- Geodetic Inverse Problem on the GRS80 ellipsoid
- Geodetic Direct Problem on the GRS80 ellipsoid
- Distance between two points (Vincenty)
- Radius of curvature for a given azimuth (Euler)
- Mean radius of the International 1924 ellipsoid
- Meridional radius of curvature (M)
- Normal radius of curvature (N)
Map projections
Grid references and EPSG helpers
Polygon utilities
from pygeodetics import Geodetic
# Initialize Geodetic class with WGS84 ellipsoid
geod = Geodetic()
lat = 59.907072474276958 # Latitude in degrees
lon = 10.754482924017791 # Longitude in degrees
h = 63.8281 # Height in meters
X, Y, Z = geod.geod2ecef(lat, lon, h)
print(f"Geodetic to ECEF:\nX: {X:.4f} m\nY: {Y:.4f} m\nZ: {Z:.4f} m")from pygeodetics import Geodetic
X, Y, Z = 3149785.9652, 598260.8822, 5495348.4927
geod = Geodetic()
lat, lon, h = geod.ecef2geod(X, Y, Z, angle_unit='deg')
print(f"ECEF to Geodetic:\nLatitude: {lat:.6f}°\nLongitude: {lon:.6f}°\nHeight: {h:.3f} m")from pygeodetics import Geodetic
X, Y, Z = 3149785.9652, 598260.8822, 5495348.4927
lat0, lon0, h0 = 58.907072, 9.75448, 63.8281
e, n, u = Geodetic().ecef2enu(X, Y, Z, lat0, lon0, h0, radians=False)
print(f"ECEF to ENU:\nEast: {e:.6f} m\nNorth: {n:.6f} m\nUp: {u:.6f} m")from pygeodetics import Geodetic
X, Y, Z = 3149785.9652, 598260.8822, 5495348.4927
lat0, lon0, h0 = 58.907072, 9.75448, 63.8281
n, e, d = Geodetic().ecef2ned(X, Y, Z, lat0, lon0, h0)
print(f"ECEF to NED:\nNorth: {n:.6f} m\nEast: {e:.6f} m\nDown: {d:.6f} m")from pygeodetics import Geodetic
from pygeodetics.Ellipsoid import GRS80
geod = Geodetic(GRS80())
lat1, lon1 = 52.2296756, 21.0122287
lat2, lon2 = 41.8919300, 12.5113300
az1, az2, distance = geod.inverse_problem(lat1, lon1, lat2, lon2, quadrant_correction=False)
print(f"Geodetic Inverse Problem:\nForward Azimuth: {az1:.6f}°\nFinal Azimuth at P2: {az2:.6f}°\nDistance: {distance:.3f} m")from pygeodetics import Geodetic
from pygeodetics.Ellipsoid import GRS80
geod = Geodetic(GRS80())
lat1, lon1 = 52.2296756, 21.0122287
az1 = -147.4628043168
d = 1316208.08334
lat2, lon2, az2 = geod.direct_problem(lat1, lon1, az1, d, quadrant_correction=True)
print(f"Geodetic Direct Problem:\nDestination Latitude: {lat2:.6f}°\nDestination Longitude: {lon2:.6f}°\nFinal Azimuth at destination: {az2:.6f}°")from pygeodetics import Geodetic
lat = 45
azimuth = 30
radius = Geodetic().radius_of_curvature(lat, azimuth, radians=False)
print(f"Radius of Curvature:\n{radius:.3f} meters")from pygeodetics import Geodetic
from pygeodetics.Ellipsoid import International1924
geod = Geodetic(International1924())
mean_radius = geod.get_mean_radius()
print(f"Mean Radius of the Ellipsoid:\n{mean_radius:.3f} meters")from pygeodetics import Geodetic
# Define the coordinates of the first point
lat1 = 52.2296756
lon1 = 21.0122287
# Define the coordinates of the second point
lat2 = 41.8919300
lon2 = 12.5113300
distances = Geodetic().distance_between_two_points(lon1, lat1, lon2, lat2, radians=False)
print(f"Distances between the two points: {distances}")from pygeodetics import Geodetic
# Compute the mean radius of the ellipsoid at a given latitude
lat = 61.456121547 # Latitude in degrees
mradius = Geodetic().mrad(lat)
print(f"Mean Radius of the Ellipsoid at Latitude {lat}°: {mradius:.3f} meters")from pygeodetics import Geodetic
# Compute the normal radius of the ellipsoid at a given latitude
lat = 61.456121547 # Latitude in degrees
mradius = Geodetic().nrad(lat)
print(f"Normal Radius of the Ellipsoid at Latitude {lat}°:\n{mradius:.3f} meters")from pygeodetics import MercatorVariantC
import numpy as np
conv = MercatorVariantC(
a=6378245.0, f=1/298.3,
latSP1=np.deg2rad(42), latFO=np.deg2rad(42), lonFO=np.deg2rad(51),
EFO=0.0, NFO=0.0
)
lon, lat = 53, 53
E, N = conv.geog_to_projected([[lon, lat]], unit="deg").ravel()
rlon, rlat = conv.projected_to_geog([[E, N]]).ravel()
print(f"Easting = {E:.2f} m\nNorthing = {N:.2f} m")
print(f"Reversed lon = {rlon:.8f}°\nReversed lat = {rlat:.8f}°")import numpy as np
from pygeodetics import TransverseMercator
np.set_printoptions(precision=8, suppress=True)
# Projection parameters (radians for origins)
lat_origin = 0.0 # φ0 (rad)
lon_origin = np.radians(9.0) # λ0 (rad), e.g. UTM zone 32
scale_factor = 0.9996 # k0
false_easting = 500000.0 # FE (m)
false_northing = 0.0 # FN (m)
# Ellipsoid (WGS84)
a = 6378137.0
f = 1 / 298.257223563
tm = TransverseMercator(
lat_origin, lon_origin, scale_factor,
false_easting, false_northing, a, f
)
# Input coordinates (lon, lat, h)
coords = np.array([
[3.0, 60.0, 100.0],
[3.2, 61.0, 102.0],
])
# Perform forward projection
easting, northing, height = tm.geog_to_projected(coords, unit="deg")
results = np.vstack((easting, northing, height)).T
print(f"\nProjected Coordinates TM:\n{results}")
# Perform inverse projection
proj_coordinates = np.vstack([easting, northing, height]).T
lon_back, lat_back, height = tm.projected_to_geog(proj_coordinates, unit="deg")
results = np.vstack((lon_back, lat_back, height)).T
print(f"\nGeographic Coordinates TM:\n{results}")A turnkey wrapper around TransverseMercator that applies the standard UTM
defaults (scale factor 0.9996, false easting 500 000 m, false northing
0 m / 10 000 000 m for N / S hemisphere) and selects the UTM zone
automatically from the longitude. Supports scalar and NumPy array inputs.
from pygeodetics import geodetic_to_utm, utm_to_geodetic
# Forward: geodetic -> UTM (auto zone, auto hemisphere)
easting, northing, zone, band = geodetic_to_utm(lat=60.0, lon=10.75)
print(f"E={easting:.3f} N={northing:.3f} zone={zone}{band}")
# Inverse: UTM -> geodetic (zone + hemisphere required)
lat, lon = utm_to_geodetic(easting, northing, zone=zone, hemisphere="N")
print(f"lat={lat:.10f}° lon={lon:.10f}°")
# Southern hemisphere works the same way (false northing 10 000 000 m
# is applied automatically on the forward call).
e_s, n_s, zone_s, band_s = geodetic_to_utm(-33.8688, 151.2093)
# Override the auto-selected zone (useful near zone borders or for the
# Norway / Svalbard exceptions, which are not auto-applied).
e, n, _, _ = geodetic_to_utm(60.0, 10.75, force_zone=33)
# Custom ellipsoid
from pygeodetics.Ellipsoid import GRS80
geodetic_to_utm(60.0, 10.75, ellipsoid=GRS80())UPS covers the polar regions where UTM is undefined (latitudes north of 84°N and south of 80°S). It uses the standard UPS defaults: scale factor 0.994, false easting / northing 2 000 000 m, central meridian 0°. The hemisphere is selected automatically from the sign of the latitude.
from pygeodetics import geodetic_to_ups, ups_to_geodetic
# Forward (auto hemisphere)
e, n, hemi = geodetic_to_ups(85.0, 0.0)
print(f"E={e:.3f} N={n:.3f} hemi={hemi}")
# Inverse (hemisphere required)
lat, lon = ups_to_geodetic(e, n, hemisphere="N")
# Both poles project to the false origin (2 000 000, 2 000 000)
geodetic_to_ups(90.0, 0.0) # -> (2_000_000, 2_000_000, 'N')
geodetic_to_ups(-90.0, 0.0) # -> (2_000_000, 2_000_000, 'S')Encode any geodetic position as an MGRS string. The function automatically chooses UTM (latitudes -80°..84°) or UPS (polar regions) and supports all standard precision levels (0..5 → 100 km..1 m).
from pygeodetics import to_mgrs, from_mgrs
# Forward at 1 m precision
to_mgrs(60.0, 10.75) # '32VNM9760352702'
to_mgrs(-33.8688, 151.2093) # '56HLH3436850948' (southern hemisphere)
to_mgrs(85.0, 0.0) # 'ZAB0000044542' (UPS north)
to_mgrs(-85.0, 45.0) # 'BFR9276792767' (UPS south)
# Lower precision
to_mgrs(60.0, 10.75, precision=3) # '32VNM976527' (100 m)
to_mgrs(60.0, 10.75, precision=0) # '32VNM' (100 km square only)
# Parsing (returns the SW corner of the cell)
lat, lon = from_mgrs('32VNM9760352702')
lat, lon = from_mgrs('32V NM 97603 52702') # spaces & lowercase acceptedWhen you need to tag projected coordinates with a CRS identifier for
downstream GIS tools (GDAL, QGIS, GeoJSON, etc.), these functions
calculate the EPSG code on the fly using the standard EPSG numbering
schemes — WGS84 UTM North is 32600 + zone, South is 32700 + zone,
NAD83 UTM North is 26900 + zone, and UPS North / South are the fixed
codes 32661 / 32761. No database, file, or network access is involved.
from pygeodetics import utm_epsg, ups_epsg
utm_epsg(32, 'N') # 32632 (WGS84 / UTM zone 32N)
utm_epsg(56, 'S') # 32756 (WGS84 / UTM zone 56S)
utm_epsg(15, 'N', datum='NAD83') # 26915 (NAD83 / UTM zone 15N)
ups_epsg('N') # 32661 (WGS84 / UPS North)
ups_epsg('S') # 32761 (WGS84 / UPS South)Round-trip conversions between geodetic, ECEF, and the local tangent
frames ENU (East-North-Up) and NED (North-East-Down). Cross-validated
against pymap3d to ~1e-12 precision.
from pygeodetics import (
geodetic2enu, enu2geodetic,
geodetic2ned, ned2geodetic,
enu2ecef, ned2ecef,
)
# Local observer at Oslo
lat0, lon0, h0 = 59.91, 10.75, 100.0
# Forward: geodetic -> ENU
e, n, u = geodetic2enu(60.0, 10.9, 250.0, lat0, lon0, h0)
print(f"ENU = ({e:.3f}, {n:.3f}, {u:.3f})")
# Inverse: ENU -> geodetic (sub-millimetre round trip)
lat, lon, h = enu2geodetic(e, n, u, lat0, lon0, h0)
# NED works the same way
n_, e_, d_ = geodetic2ned(60.0, 10.9, 250.0, lat0, lon0, h0)
lat, lon, h = ned2geodetic(n_, e_, d_, lat0, lon0, h0)
# Arrays are supported transparently
import numpy as np
lats = np.array([60.0, 60.1, 60.2])
lons = np.array([10.8, 10.9, 11.0])
hs = np.array([100.0, 200.0, 300.0])
e, n, u = geodetic2enu(lats, lons, hs, lat0, lon0, h0)Convert between AER (azimuth clockwise from north in degrees, elevation above horizon in degrees, slant range in metres) and any of geodetic, ECEF, ENU or NED.
from pygeodetics import (
geodetic2aer, aer2geodetic,
enu2aer, aer2enu,
ned2aer, aer2ned,
ecef2aer, aer2ecef,
)
# Look angle from observer to a satellite ground projection
lat0, lon0, h0 = 59.91, 10.75, 100.0
az, el, rng = geodetic2aer(60.0, 10.9, 250.0, lat0, lon0, h0)
print(f"Azimuth = {az:.3f}°, Elevation = {el:.3f}°, Range = {rng:.1f} m")
# Inverse: place a target at a given look-angle
lat, lon, h = aer2geodetic(az, el, rng, lat0, lon0, h0)Compute perimeter, ellipsoidal area, centroid, bounding box, geodesic densification and geodesic interpolation for a polygon defined by latitude/longitude vertices. Polygons are auto-closed if not already.
from pygeodetics import (
polygon_perimeter, polygon_area, polygon_centroid,
polygon_bounds, polygon_densify, geodesic_interpolate,
)
# Polygon over southern Norway
lats = [60.0, 60.0, 60.5, 60.5]
lons = [10.0, 11.0, 11.0, 10.0]
print("Perimeter:", polygon_perimeter(lats, lons), "m") # Vincenty sum
print("Area :", polygon_area(lats, lons), "m^2") # ellipsoidal
print("Centroid :", polygon_centroid(lats, lons)) # (lat, lon)
print("Bounds :", polygon_bounds(lats, lons)) # (lat_min, lon_min, lat_max, lon_max)
# Densify: insert geodesic vertices so no segment exceeds 5 km. The
# original vertices are preserved exactly; new ones lie on the true
# geodesic of each edge. Useful before projecting a sparsely sampled
# polygon to a flat map (avoids the straight-cartographic-edge artefact).
dense_lats, dense_lons = polygon_densify(lats, lons, max_segment_length=5000.0)
print(f"Densified vertex count: {len(dense_lats)} (from {len(lats)})")
# Sample 50 evenly spaced points along a single geodesic
sample_lats, sample_lons = geodesic_interpolate(60.0, 10.0, 70.0, 25.0, n_points=50)
# True geodesic-polygon area for a continental-scale polygon
# (the rhumb-edge formula is off by ~1% at this size; geodesic mode
# matches GeographicLib to ~1 ppm).
us_lats = [49.0, 49.0, 25.0, 25.0]
us_lons = [-125.0, -67.0, -67.0, -125.0]
print("US bbox area:", polygon_area(us_lats, us_lons, geodesic=True), "m^2")Note on area — by default
polygon_areaintegrates the authalic latitude along rhumb-line edges (Green's theorem on the ellipsoid). It matchespyproj.Geod.polygon_area_perimeterto better than 1 ppm for small and medium polygons (degrees-scale). For continental-scale or hemispheric polygons the rhumb / geodesic discrepancy can grow to a fraction of a percent — passgeodesic=Trueto compute the true geodesic-polygon area instead. In that mode the polygon is internally densified (default cap 10 km per segment, controllable viamax_segment_length) before the area integral, which converges to the GeographicLib reference to better than 1 ppm even for hemispheric polygons.
This section provides the mathematical foundation for the computations performed in PyGeodetics.
The conversion from geodetic coordinates
where:
-
$N = \frac{a}{\sqrt{1 - e^2 \sin^2\phi}}$ is the prime vertical radius of curvature, -
$a$ is the semi-major axis, -
$b$ is the semi-minor axis, -
$e^2 = \frac{a^2 - b^2}{a^2}$ is the first eccentricity squared, -
$\phi$ is the geodetic latitude, -
$\lambda$ is the geodetic longitude, -
$h$ is the height above the ellipsoid.
PyGeodetics implements three solvers for the inverse problem
-
Compute the longitude:
$$\lambda=\arctan2(Y, X)$$ -
Compute the intermediate values:
$$p=\sqrt{X^2 + Y^2}, \quad \theta = \arctan2(Z,a,; p,b)$$ -
Compute the latitude:
$$\phi=\arctan2!\left(Z + e'^2 b \sin^3\theta,; p - e^2 a \cos^3\theta\right)$$ ,where$e'^2=\frac{a^2 - b^2}{b^2}$ is the second eccentricity squared. -
Compute the height:
$$h=\frac{p}{\cos\phi}-N$$
The geodetic inverse problem calculates the geodesic distance
-
Compute the reduced latitudes (also called parametric latitudes):
$$U_1 = \arctan\left((1 - f) \tan\phi_1\right), \quad U_2 = \arctan\left((1 - f) \tan\phi_2\right)$$ ,where$f = \frac{a - b}{a}$ is the flattening. -
Set
$L = \lambda_2 - \lambda_1$ and initialise$\lambda = L$ . Iterate$\lambda$ (longitude on the auxiliary sphere) together with the angular distance$\sigma$ :$$\sigma = \arctan2!\left(\sqrt{(\cos U_2 \sin\lambda)^2 + (\cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos\lambda)^2},; \sin U_1 \sin U_2 + \cos U_1 \cos U_2 \cos\lambda\right)$$ The full update equation for
$\lambda$ and the auxiliary quantities$\sin\alpha$ ,$\cos 2\sigma_m$ ,$C$ are given in Section 7. -
Once converged, compute the forward and reverse azimuths:
$$\alpha_1 = \arctan2\left(\cos U_2 \sin\lambda,; \cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos\lambda\right)$$ $$\alpha_2 = \arctan2\left(\cos U_1 \sin\lambda,; -\sin U_1 \cos U_2 + \cos U_1 \sin U_2 \cos\lambda\right)$$ -
Compute the geodesic distance using the series expansion in
$u^2 = \cos^2\alpha,(a^2 - b^2)/b^2$ :$$s = b,A,(\sigma - \Delta\sigma)$$ with
$A$ ,$B$ and$\Delta\sigma$ defined in Section 7. The simpler one-term form$s \approx b(\sigma - \tfrac{f}{4} \sin\sigma \cos(2\sigma_m + \sigma))$ is not used here — it loses several kilometres of accuracy on continental baselines.
The radius of curvature in the meridian (
The transformation from ECEF to local East-North-Up (ENU) coordinates is given by:
where
PyGeodetics returns the IUGG arithmetic mean radius:
This is one of several conventional definitions (others include the authalic and volumetric mean radii); it is the simple arithmetic mean of the three semi-axes
Vincenty's formula is used to calculate the geodesic distance between two points on an ellipsoid. The method iteratively solves for the distance
-
Reduced Latitude: Compute the reduced latitudes:
$$U_1=\arctan\left((1 - f) \cdot \tan\phi_1\right), \quad U_2=\arctan\left((1 - f) \cdot\tan\phi_2\right)$$ ,where
$f = \frac{a - b}{a}$ is the flattening. -
Longitude Difference: Compute the difference in longitudes:
$$L=\lambda_2 - \lambda_1$$ -
Iterative Solution: Initialize
$\lambda = L$ and iteratively solve for$\lambda$ using:$$\lambda=L + (1 - C) f \sin\alpha \left[\sigma + C \sin\sigma \left(\cos2\sigma_m + C \cos\sigma \left(-1 + 2 \cos^2 2\sigma_m\right)\right)\right]$$ where:
-
$\sigma = \arctan2\left(\sqrt{(\cos U_2 \sin\lambda)^2 + (\cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos\lambda)^2}, \sin U_1 \sin U_2 + \cos U_1 \cos U_2 \cos\lambda\right)$ , -
$\sin\alpha = \frac{\cos U_1 \cos U_2 \sin\lambda}{\sin\sigma}$ , -
$\cos2\sigma_m = \cos\sigma - \frac{2 \sin U_1 \sin U_2}{\cos^2\alpha}$ , -
$C = \frac{f}{16} \cos^2\alpha \left(4 + f \left(4 - 3 \cos^2\alpha\right)\right)$ .
The iteration stops when
$|\lambda - \lambda_{\text{prev}}| < \text{tolerance}$ . -
-
Geodesic Distance: Compute the geodesic distance:
$$s=b A \left[\sigma - \delta\sigma\right]$$ ,where:-
$u^2=\frac{\cos^2\alpha (a^2 - b^2)}{b^2}$ , -
$A=1 + \frac{u^2}{16384} \left(4096 + u^2 \left(-768 + u^2 (320 - 175 u^2)\right)\right)$ , -
$B=\frac{u^2}{1024} \left(256 + u^2 \left(-128 + u^2 (74 - 47 u^2)\right)\right)$ , -
$\delta\sigma=B \sin\sigma \left[\cos2\sigma_m + \frac{B}{4} \left(\cos\sigma \left(-1 + 2 \cos^2 2\sigma_m\right) - \frac{B}{6} \cos2\sigma_m \left(-3 + 4 \sin^2\sigma\right) \left(-3 + 4 \cos^2 2\sigma_m\right)\right)\right]$ .
-
-
Azimuths: Compute the forward and reverse azimuths:
$$\alpha_1 = \arctan2\left(\cos U_2 \sin\lambda, \cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos\lambda\right)$$ $$\alpha_2 = \arctan2\left(\cos U_1 \sin\lambda, -\sin U_1 \cos U_2 + \cos U_1 \sin U_2 \cos\lambda\right)$$
Vincenty's formula is highly accurate for most geodetic calculations but may fail to converge for nearly antipodal points.
The meridional radius of curvature, denoted as
where:
-
$a$ is the semi-major axis of the ellipsoid, -
$e^2 = \frac{a^2 - b^2}{a^2}$ is the first eccentricity squared, -
$\phi$ is the geodetic latitude.
This formula accounts for the flattening of the ellipsoid and the variation in curvature with latitude.
The normal radius of curvature, denoted as
where:
-
$a$ is the semi-major axis of the ellipsoid, -
$e^2 = \frac{a^2 - b^2}{a^2}$ is the first eccentricity squared, -
$\phi$ is the geodetic latitude.
The value of
The grid convergence, denoted as
The grid convergence
where:
-
$\Delta\lambda = \lambda - \lambda_0$ is the longitude difference from the central meridian, -
$\epsilon^2 = \frac{e^2}{1 - e^2} \cos^2\phi$ is the second eccentricity squared, -
$e^2 = \frac{a^2 - b^2}{a^2}$ is the first eccentricity squared, -
$a$ is the semi-major axis, -
$b$ is the semi-minor axis, -
$\phi$ is the geodetic latitude, -
$\lambda$ is the geodetic longitude, -
$\lambda_0$ is the central meridian.
The grid convergence
where:
-
$\phi_f$ is the footpoint latitude, computed iteratively, -
$N_f = \frac{a}{\sqrt{1 - e^2 \sin^2\phi_f}}$ is the normal radius of curvature at the footpoint latitude, -
$\epsilon_f^2 = \frac{e^2}{1 - e^2} \cos^2\phi_f$ is the second eccentricity squared at the footpoint latitude, -
$x$ is the easting coordinate (adjusted for false easting), -
$y$ is the northing coordinate.
The scale factor, denoted as
The scale factor
where:
-
$\Delta\lambda = \lambda - \lambda_0$ is the longitude difference from the central meridian, -
$\epsilon^2 = \frac{e^2}{1 - e^2} \cos^2\phi$ is the second eccentricity squared, -
$e^2 = \frac{a^2 - b^2}{a^2}$ is the first eccentricity squared, -
$a$ is the semi-major axis, -
$b$ is the semi-minor axis, -
$\phi$ is the geodetic latitude, -
$\lambda$ is the geodetic longitude, -
$\lambda_0$ is the central meridian.
The scale factor
where:
-
$M_f = \frac{a(1 - e^2)}{(1 - e^2 \sin^2\phi_f)^{3/2}}$ is the meridional radius of curvature at the footpoint latitude, -
$N_f = \frac{a}{\sqrt{1 - e^2 \sin^2\phi_f}}$ is the normal radius of curvature at the footpoint latitude, -
$\phi_f$ is the footpoint latitude, computed iteratively, -
$x$ is the easting coordinate (adjusted for false easting), -
$y$ is the northing coordinate.
For a spherical Earth model, the scale factor
where:
-
$x$ is the easting coordinate (adjusted for false easting), -
$R$ is the radius of the sphere.
This section explains the mathematical foundation for the projections method supported by the library. In this release supported methods are:
- Mercator Variant C
- Transverse Mercator projections.
The Mercator Variant C projection is a cylindrical map projection that preserves angles, making it conformal. This projection is a variant of the Mercator projection, with specific parameters defined for its implementation.
The following table summarizes the key parameters of the Mercator Variant C projection as defined by EPSG:
| Parameter Name | Parameter EPSG Code | Sign Reversal | Description |
|---|---|---|---|
| Latitude of 1st standard parallel | 8823 | No | For a conic projection with two standard parallels, this is the latitude of one of the parallels of intersection of the cone with the ellipsoid. Scale is true along this parallel. |
| Longitude of natural origin | 8802 | No | The longitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. |
| Latitude of false origin | 8821 | No | The latitude of the point which is not the natural origin and at which grid coordinate values false easting and false northing are defined. |
| Easting at false origin | 8826 | No | The easting value assigned to the false origin. |
| Northing at false origin | 8827 | No | The northing value assigned to the false origin. |
The forward projection transforms geographic coordinates
- Constants:
$k_0 = \frac{\cos \phi_1}{\sqrt{1 - e^2 \sin^2 \phi_1}}$ - All logarithms below are natural logarithms (base
$\mathrm{e}$ ). To avoid clashing with the eccentricity$e$ , the natural-log base is denoted$\mathrm{e}$ in this section.
-
Compute the projection constant
$k_0$ :$$k_0 = \frac{\cos \phi_1}{\sqrt{1 - e^2 \sin^2 \phi_1}}$$ where:
-
$\phi_1$ = latitude of the first standard parallel (absolute value, positive), -
$e^2 = \tfrac{a^2 - b^2}{a^2}$ = eccentricity squared, -
$a$ = semi-major axis,$b$ = semi-minor axis.
-
-
Compute the meridional arc at latitude
$\phi$ :$$M(\phi) = a k_0 ,\ln \left[ \tan!\left(\tfrac{\pi}{4} + \tfrac{\phi}{2}\right) \left(\frac{1 - e \sin\phi}{1 + e \sin\phi}\right)^{\tfrac{e}{2}} \right]$$ -
Compute Easting and Northing:
$$E = E_0 + a k_0 (\lambda - \lambda_0)$$ $$N = N_0 - M(\phi_0) + M(\phi)$$ where:
-
$(\lambda_0, \phi_0)$ = longitude/latitude of the false origin, -
$(E_0, N_0)$ = false easting and northing.
-
The inverse projection transforms projected coordinates
-
Longitude:
$$\lambda = \lambda_0 + \frac{E - E_0}{a k_0}$$ -
Latitude (iterative series):
Define:
$$\chi = \frac{\pi}{2} - 2 \arctan !\left( \exp!\left[-\frac{N - N_0 + M(\phi_0)}{a k_0}\right] \right)$$ Then:
$$\phi=\chi+\left(\tfrac{e^2}{2} + \tfrac{5 e^4}{24} + \tfrac{e^6}{12} + \tfrac{13 e^8}{360}\right) \sin(2\chi)+\left(\tfrac{7 e^4}{48} +\tfrac{29 e^6}{240} + \tfrac{811 e^8}{11520}\right) \sin(4\chi)+\left(\tfrac{7 e^6}{120} + \tfrac{81 e^8}{1120}\right) \sin(6\chi)+\left(\tfrac{4279 e^8}{161280}\right) \sin(8\chi)$$
The Transverse Mercator (TM) projection is a conformal map projection widely used for large-scale mapping, such as the Universal Transverse Mercator (UTM) system.
The following table summarizes the key parameters of the Transverse Mercator projection as defined by EPSG:
| Parameter Name | Parameter EPSG Code | Sign Reversal | Description |
|---|---|---|---|
| Latitude of natural origin | 8801 | No | The latitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. |
| Longitude of natural origin | 8802 | No | The longitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. |
| Scale factor at natural origin | 8805 | No | The factor by which the map grid is reduced or enlarged during the projection process, defined by its value at the natural origin. |
| False easting | 8806 | No | The value assigned to the abscissa (east or west) axis of the projection grid at the natural origin to avoid negative coordinates over parts of the mapped area. |
| False northing | 8807 | No | The value assigned to the ordinate (north or south) axis of the projection grid at the natural origin to avoid negative coordinates over parts of the mapped area. |
The forward projection transforms geographic coordinates
-
Calculate constants for the projection:
$$n = \frac{f}{2 - f}$$ $$B = \frac{a}{1 + n} \left(1 + \frac{n^2}{4} + \frac{n^4}{64}\right)$$ $$h_1=\frac{n}{2} - \frac{2}{3}n^2 + \frac{5}{16}n^3 + \frac{41}{180}n^4$$ $$h_2=\frac{13}{48}n^2 - \frac{3}{5}n^3 + \frac{557}{1440}n^4$$ $$h_3=\frac{61}{240}n^3 - \frac{103}{140}n^4$$ $$h_4=\frac{49561}{161280}n^4$$ -
Compute the meridional arc distance from the equator to the projection origin:
If
$\phi_0 = 0$ $$M_0 = 0$$ or if
$\phi_0 = 90^\circ , (\pi/2 , \text{radians})$ $$M_0 = B \cdot \frac{\pi}{2}$$ or if
$\phi_0 = -90^\circ , (-\pi/2 , \text{radians})$ $$M_0 = B \cdot \left(-\frac{\pi}{2}\right)$$ Otherwise:
$$Q_0 = \sinh^{-1}(\tan\phi_0) - e \tanh^{-1}(e \sin\phi_0)$$ $$\beta_0 = \tan^{-1}(\sinh Q_0)$$ $$\xi_{0,0} = \sin^{-1}(\sin\beta_0)$$ $$\xi_0 = \xi_{0,0} + h_1 \sin(2\xi_{0,0}) + h_2 \sin(4\xi_{0,0}) + h_3 \sin(6\xi_{0,0}) + h_4 \sin(8\xi_{0,0})$$ $$M_0 = B \xi_0$$ -
Compute intermediate values for the given latitude
$\phi$ :$$Q = \sinh^{-1}(\tan\phi) - e \tanh^{-1}(e \sin\phi)$$ $$\beta = \tan^{-1}(\sinh Q)$$ $$\eta_0 = \tanh^{-1}(\cos\beta \sin(\lambda - \lambda_0))$$ $$\xi_0 = \sin^{-1}(\sin\beta \cosh\eta_0)$$ $$\xi = \xi_0 + h_1 \sin(2\xi_0) \cosh(2\eta_0) + h_2 \sin(4\xi_0) \cosh(4\eta_0) + h_3 \sin(6\xi_0) \cosh(6\eta_0) + h_4 \sin(8\xi_0) \cosh(8\eta_0)$$ $$\eta = \eta_0 + h_1 \cos(2\xi_0) \sinh(2\eta_0) + h_2 \cos(4\xi_0) \sinh(4\eta_0) + h_3 \cos(6\xi_0) \sinh(6\eta_0) + h_4 \cos(8\xi_0) \sinh(8\eta_0)$$ -
Compute the easting and northing:
$$E = E_0 + k_0 B \eta$$ $$N = N_0 + k_0 \left(B \xi - M_0\right)$$
The inverse projection transforms projected coordinates
-
Compute intermediate values:
$$\eta' = \frac{E - E_0}{k_0 B}, \quad \xi' = \frac{(N - N_0) + k_0 M_0}{k_0 B}$$ -
Iteratively compute
$\xi_0'$ and$\eta_0'$ :$$\xi_0' = \xi' - \left(h_1 \sin(2\xi') \cosh(2\eta') + h_2 \sin(4\xi') \cosh(4\eta') + h_3 \sin(6\xi') \cosh(6\eta') + h_4 \sin(8\xi') \cosh(8\eta')\right)$$ $$\eta_0' = \eta' - \left(h_1 \cos(2\xi') \sinh(2\eta') + h_2 \cos(4\xi') \sinh(4\eta') + h_3 \cos(6\xi') \sinh(6\eta') + h_4 \cos(8\xi') \sinh(8\eta')\right)$$ -
Compute
$\beta'$ and$Q'$ :$$\beta' = \sin^{-1}\left(\frac{\sin\xi_0'}{\cosh\eta_0'}\right)$$ $$Q' = \sinh^{-1}(\tan\beta')$$ -
Iteratively compute latitude
$\phi$ :$$Q'' = Q' + e \tanh^{-1}(e \tanh Q')$$ Repeat until the change in
$Q''$ is insignificant. Then:$$\phi = \tan^{-1}(\sinh Q'')$$ -
Compute longitude
$\lambda$ :$$\lambda = \lambda_0 + \sin^{-1}\left(\frac{\tanh\eta_0'}{\cos\beta'}\right)$$
- These formulas follow EPSG Guidance Note 7-2.
This project is licensed under the MIT License.