Metadata-Version: 2.4
Name: pygeodetics
Version: 1.1.0
Summary: A Python library for geodetic computations
Author-email: Per Helge Aarnes <per.helge.aarnes@gmail.com>
License: MIT License
        
        Copyright (c) 2025 Per Helge Aarnes
        
        Permission is hereby granted, free of charge, to any person obtaining a copy
        of this software and associated documentation files (the "Software"), to deal
        in the Software without restriction, including without limitation the rights
        to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
        copies of the Software, and to permit persons to whom the Software is
        furnished to do so, subject to the following conditions:
        
        The above copyright notice and this permission notice shall be included in all
        copies or substantial portions of the Software.
        
        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
        LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
        OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
        SOFTWARE.
        
Project-URL: Homepage, https://github.com/paarnes/pygeodetics
Project-URL: Repository, https://github.com/paarnes/pygeodetics
Project-URL: Issues, https://github.com/paarnes/pygeodetics/issues
Keywords: geodesy,coordinates,geodetic,ECEF,ENU,conversion,ellipsoid,geodetic-direct,geodetic-inverse,vincenty-formula,transverse-mercator,map-projections,scale-factor,grid-convergence
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Developers
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.21
Dynamic: license-file

<p align="center">
    <img src="https://github.com/paarnes/pygeodetics/blob/main/docs/icon/pygeodetics.png?raw=true" alt="PyGeodetics Logo" width="700">
</p>

<p align="center">
    <a href="https://github.com/paarnes/pygeodetics/actions/workflows/tests.yml"><img alt="Tests" src="https://github.com/paarnes/pygeodetics/actions/workflows/tests.yml/badge.svg"></a>
    <a href="https://pypi.org/project/pygeodetics/"><img alt="PyPI" src="https://img.shields.io/pypi/v/pygeodetics?color=blue"></a>
    <a href="https://pypi.org/project/pygeodetics/"><img alt="Python versions" src="https://img.shields.io/pypi/pyversions/pygeodetics"></a>
    <a href="https://github.com/paarnes/pygeodetics/blob/main/LICENSE"><img alt="License: MIT" src="https://img.shields.io/github/license/paarnes/pygeodetics"></a>
    <a href="https://github.com/paarnes/pygeodetics/issues"><img alt="GitHub issues" src="https://img.shields.io/github/issues/paarnes/pygeodetics"></a>
</p>

## Introduction

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_geodetic`
  with 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=True` mode 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 `Ellipsoid` instance.
- **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.

## Installation
```sh
pip install pygeodetics
```

## Table of contents — code examples

**Coordinate system conversions**
- [Geodetic to ECEF](#geodetic-to-ecef)
- [ECEF to Geodetic](#ecef-to-geodetic)
- [ECEF to ENU](#ecef-to-enu)
- [ECEF to NED](#ecef-to-ned)

**Local tangent frames (ENU / NED / AER)**
- [Bidirectional ENU / NED conversions](#bidirectional-enu--ned-conversions)
- [Azimuth / Elevation / Range (AER) conversions](#azimuth--elevation--range-aer-conversions)

**Geodesic problems on the ellipsoid**
- [Geodetic Inverse Problem on the GRS80 ellipsoid](#geodetic-inverse-problem-on-the-grs80-ellipsoid)
- [Geodetic Direct Problem on the GRS80 ellipsoid](#geodetic-direct-problem-on-the-grs80-ellipsoid)
- [Distance between two points (Vincenty)](#calculate-the-distance-between-two-points-on-the-ellipsoid-vincenty-formula)
- [Radius of curvature for a given azimuth (Euler)](#radius-of-curvature-for-a-given-azimuth-using-eulers-equation)
- [Mean radius of the International 1924 ellipsoid](#calculate-the-mean-radius-of-the-international1924-ellipsoid)
- [Meridional radius of curvature (M)](#calculate-the-meridional-radius-of-curvature-m-at-a-given-latitude)
- [Normal radius of curvature (N)](#calculate-the-normal-radius-of-curvature-n-at-a-given-latitude)

**Map projections**
- [Mercator Variant C](#use-of-mercator-variant-c-projection)
- [Transverse Mercator](#use-of-transverse-mercator-projection)
- [UTM convenience API](#utm-convenience-api-geodetic_to_utm--utm_to_geodetic)
- [UPS (Universal Polar Stereographic)](#ups-universal-polar-stereographic-for-the-poles)

**Grid references and EPSG helpers**
- [MGRS encoding and parsing](#mgrs-encoding-and-parsing)
- [EPSG helpers](#epsg-helpers)

**Polygon utilities**
- [Geodesic polygon utilities](#geodesic-polygon-utilities)

## Usage Examples

### Coordinate system conversions

#### Geodetic to ECEF
```python
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")

```

#### ECEF to Geodetic
```python
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")

```

#### ECEF to ENU
```python
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")

```

#### ECEF to NED
```python
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")

```

### Geodesic problems on the ellipsoid

#### Geodetic Inverse Problem on the GRS80 ellipsoid
```python
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")

```

#### Geodetic Direct Problem on the GRS80 ellipsoid
```python
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}°")

```

#### Radius of Curvature for a given Azimuth using Euler's equation.
```python
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")

```

#### Calculate the mean Radius of the International1924 Ellipsoid
```python
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")

```

#### Calculate the distance between two points on the ellipsoid (Vincenty formula)
```python
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}")

```

#### Calculate the meridional radius of curvature (M) at a given latitude

```python
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")
```


#### Calculate the normal radius of curvature (N) at a given latitude.

```python
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")
```

### Map projections

#### Use of Mercator Variant C projection
```python

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}°")
```


#### Use of Transverse Mercator projection

```python
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}")

```


#### UTM convenience API (`geodetic_to_utm` / `utm_to_geodetic`)

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.

```python
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 (Universal Polar Stereographic) for the poles

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.

```python
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')
```


### Grid references and EPSG helpers

#### MGRS encoding and parsing

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).

```python
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 accepted
```


#### EPSG helpers

When 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.

```python
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)
```

### Local tangent frames (ENU / NED / AER)

#### Bidirectional ENU / NED conversions

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.

```python
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)
```

#### Azimuth / Elevation / Range (AER) conversions

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.

```python
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)
```

### Polygon utilities on the ellipsoid

#### Geodesic polygon utilities

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.

```python
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_area` integrates the authalic
> latitude along rhumb-line edges (Green's theorem on the ellipsoid).
> It matches `pyproj.Geod.polygon_area_perimeter` to 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 — pass `geodesic=True` to compute the true
> geodesic-polygon area instead. In that mode the polygon is internally
> densified (default cap 10 km per segment, controllable via
> `max_segment_length`) before the area integral, which converges to
> the GeographicLib reference to better than 1 ppm even for
> hemispheric polygons.


## Math and the Theory Basis

This section provides the mathematical foundation for the computations performed in PyGeodetics. 

### 1. Geodetic to ECEF Conversion

The conversion from geodetic coordinates $(\phi, \lambda, h)$ to Earth-Centered, Earth-Fixed (ECEF) coordinates $(X, Y, Z)$ is given by:

$$
\begin{aligned}
    X &= (N + h) \cos\phi \cos\lambda, \\
    Y &= (N + h) \cos\phi \sin\lambda, \\
    Z &= \left( \frac{b^2}{a^2} N + h \right) \sin\phi,
\end{aligned}
$$

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.

### 2. ECEF to Geodetic Conversion

PyGeodetics implements three solvers for the inverse problem $(X, Y, Z) \to (\phi, \lambda, h)$ — an iterative Heiskanen–Moritz solver, Bowring's closed-form method, and Vermeille's closed-form method. The recommended (default) solver is Bowring's method, summarised below:

1. Compute the longitude:

   $$\lambda=\arctan2(Y, X)$$

2. Compute the intermediate values:

   $$p=\sqrt{X^2 + Y^2}, \quad \theta = \arctan2(Z\,a,\; p\,b)$$

3. 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.

4. Compute the height:

   $$h=\frac{p}{\cos\phi}-N$$

### 3. Geodetic Inverse Problem

The geodetic inverse problem calculates the geodesic distance $s$ and azimuths $(\alpha_1, \alpha_2)$ between two points $(\phi_1, \lambda_1)$ and $(\phi_2, \lambda_2)$. PyGeodetics solves it with Vincenty's formulae, which iterate on the auxiliary sphere using the *reduced* latitudes:

1. 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.

2. 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](#7-vincentys-formula-for-geodesic-distance).

3. 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)$$

4. 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](#7-vincentys-formula-for-geodesic-distance). 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.

### 4. Radius of Curvature

The radius of curvature in the meridian ($M$) and the prime vertical ($N$) are given by:

$$
M = \frac{a(1 - e^2)}{(1 - e^2 \sin^2\phi)^{3/2}}, \quad N = \frac{a}{\sqrt{1 - e^2 \sin^2\phi}}.
$$

### 5. Local ENU Coordinates

The transformation from ECEF to local East-North-Up (ENU) coordinates is given by:

$$
\begin{bmatrix}
e \\
n \\
u
\end{bmatrix}=\begin{bmatrix}
-\sin\lambda & \cos\lambda & 0 \\
-\sin\phi\cos\lambda & -\sin\phi\sin\lambda & \cos\phi \\
\cos\phi\cos\lambda & \cos\phi\sin\lambda & \sin\phi
\end{bmatrix}
\begin{bmatrix}
X - X_0 \\
Y - Y_0 \\
Z - Z_0
\end{bmatrix}
$$



where $(X_0, Y_0, Z_0)$ are the ECEF coordinates of the reference point.

### 6. Mean Radius of the Ellipsoid

PyGeodetics returns the IUGG arithmetic mean radius:

$$
R_1 = \frac{2a + b}{3}.
$$

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 $(a, a, b)$.

### 7. Vincenty's Formula for Geodesic Distance

Vincenty's formula is used to calculate the geodesic distance between two points on an ellipsoid. The method iteratively solves for the distance $s$ and azimuths $(\alpha_1, \alpha_2)$ between two points $(\phi_1, \lambda_1)$ and $(\phi_2, \lambda_2)$.

#### Steps:

1. **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.

2. **Longitude Difference**:
   Compute the difference in longitudes:

   $$L=\lambda_2 - \lambda_1$$

3. **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}$.

4. **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]$.

5. **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.

### 8. Meridional Radius of Curvature (M)

The meridional radius of curvature, denoted as $M$, represents the radius of curvature in the north-south direction along a meridian. It is computed as:

$$M=\frac{a(1 - e^2)}{(1 - e^2 \sin^2\phi)^{3/2}}$$

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.

### 9. Normal Radius of Curvature (N)

The normal radius of curvature, denoted as $N$, represents the radius of curvature in the east-west direction perpendicular to the meridian. It is computed as:

$$
N = \frac{a}{\sqrt{1 - e^2 \sin^2\phi}},
$$

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 $N$ varies with latitude due to the ellipsoidal shape of the Earth, being largest at the equator and smallest at the poles.

### 10. Grid Convergence in the Transverse Mercator Projection

The grid convergence, denoted as $\gamma$, is the angular difference between grid north and true north in the Transverse Mercator (TM) projection. It can be computed using either geographic coordinates or projected coordinates.

#### 10.1 Grid Convergence Using Geographic Coordinates

The grid convergence $\gamma$ at a point $(\phi, \lambda)$ is given by:

$$
\gamma = \Delta\lambda \sin\phi + \frac{\Delta\lambda^3}{3} \sin\phi \cos^2\phi (1 + 3\epsilon^2 + 2\epsilon^4) + \frac{\Delta\lambda^5}{15} \sin\phi \cos^4\phi (2 - \tan^2\phi),
$$

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.

#### 10.2 Grid Convergence Using Projected Coordinates

The grid convergence $\gamma$ at a point $(x, y)$ in projected coordinates is given by:

$$
\gamma = \frac{x \tan\phi_f}{N_f} - \frac{x^3 \tan\phi_f}{3 N_f^3} \left(1 + \tan^2\phi_f - \epsilon_f^2 - 2\epsilon_f^4\right),
$$

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.

### 11. Scale Factor in the Transverse Mercator Projection

The scale factor, denoted as $k$, describes the distortion of distances in the Transverse Mercator projection. It can be computed using either geographic coordinates or projected coordinates.

#### 11.1 Scale Factor Using Geographic Coordinates

The scale factor $k$ at a point $(\phi, \lambda)$ is given by:

$$
k = 1 + \frac{\Delta\lambda^2}{2} \cos^2\phi (1 + \epsilon^2) + \frac{\Delta\lambda^4}{24} \cos^4\phi (5 + 4\tan^2\phi),
$$

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.

#### 11.2 Scale Factor Using Projected Coordinates

The scale factor $k$ at a point $(x, y)$ in projected coordinates is given by:

$$
k = 1 + \frac{x^2}{2 M_f N_f} + \frac{x^4}{24 N_f^4},
$$

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.

#### 11.3 Scale Factor for a Sphere

For a spherical Earth model, the scale factor $k$ is given by:

$$
k = \cosh\left(\frac{x}{R}\right),
$$

where:
- $x$ is the easting coordinate (adjusted for false easting),
- $R$ is the radius of the sphere.

### 12. Projections

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.

#### 12.1 Mercator Variant C Projection

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.

##### Parameters of the Mercator Variant C Projection

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.                                                                                                                                                                |

### 12.1.1 Forward Projection (Geographic → Projected)

The forward projection transforms geographic coordinates $(\lambda, \phi)$ (longitude, latitude) to projected coordinates $(E, N)$ (easting, northing).  

- 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.

1. **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.

2. **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]$$

3. **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.

---

### 12.1.2 Inverse Projection (Projected → Geographic)

The inverse projection transforms projected coordinates $(E, N)$ back to geographic $(\lambda, \phi)$.

1. **Longitude:**

   $$\lambda = \lambda_0 + \frac{E - E_0}{a k_0}$$

2. **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)$$

---


#### 12.2 Transverse Mercator Projection

The Transverse Mercator (TM) projection is a conformal map projection widely used for large-scale mapping, such as the Universal Transverse Mercator (UTM) system.

#### Parameters of the Transverse Mercator Projection

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.                                              |


##### 12.2.1 Forward Projection (Geographic to Projected)

The forward projection transforms geographic coordinates $(\lambda, \phi)$ to projected coordinates $(E, N)$ as follows:

1. **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$$

2. **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$$

3. **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)$$

4. **Compute the easting and northing**:

   $$E = E_0 + k_0 B \eta$$

   $$N = N_0 + k_0 \left(B \xi - M_0\right)$$

##### 12.2.2 Inverse Projection (Projected to Geographic)

The inverse projection transforms projected coordinates $(E, N)$ back to geographic coordinates $(\lambda, \phi)$ as follows:

1. **Compute intermediate values**:

   $$\eta' = \frac{E - E_0}{k_0 B}, \quad \xi' = \frac{(N - N_0) + k_0 M_0}{k_0 B}$$

2. **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)$$

3. **Compute $\beta'$ and $Q'$**:

   $$\beta' = \sin^{-1}\left(\frac{\sin\xi_0'}{\cosh\eta_0'}\right)$$

   $$Q' = \sinh^{-1}(\tan\beta')$$

4. **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'')$$

5. **Compute longitude $\lambda$**:

   $$\lambda = \lambda_0 + \sin^{-1}\left(\frac{\tanh\eta_0'}{\cos\beta'}\right)$$


### EPSG Reference

- These formulas follow **EPSG Guidance Note 7-2**.  

## License
This project is licensed under the MIT License.

