CoordTransform Cadastral & High-Precision Python

Geodetic Conversion Math: Ellipsoid to Cartesian for Survey-Grade Pipelines

The forward conversion from geodetic coordinates (latitude, longitude, ellipsoidal height) to Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates (X, Y, Z) is the first deterministic stage of nearly every workflow in Algorithmic Math & Geodetic Workflows: it lifts a position off the curved reference ellipsoid into a rectangular frame where datum shifts and least-squares adjustments become linear algebra. Unlike a map projection, which introduces zone-dependent scale distortion, the ellipsoid-to-Cartesian mapping is closed-form, non-iterative, and preserves the exact geometric definition of the ellipsoid. For cadastral registration, national mapping frameworks, and high-precision GIS infrastructure, this step must enforce strict numerical determinism, explicit precision boundaries, and full transformation traceability so that epoch consistency and metadata integrity survive across distributed compute nodes.

Forward geodetic-to-ECEF conversion pipeline Flowchart: geodetic input latitude, longitude and ellipsoidal height is cast to float64 radians, then a bounds-and-finite gate either routes invalid rows to a NaN sentinel or continues; the ellipsoid constant e-squared and the prime-vertical radius N are computed, the X, Y, Z triple is assembled in float64, a Bowring inverse round-trip is run, and a residual-versus-tolerance gate either flags the row for review or commits the ECEF triple. Geodetic input φ, λ, h deg, deg, m · EPSG:4979 Cast float64 → radians in bounds & finite? Route to NaN sentinel Ellipsoid constants e² = 2f − f² Prime-vertical radius N = a ⁄ √(1 − e² sin²φ) Assemble ECEF X, Y, Z (float64) Bowring inverse round-trip check residual ≤ tol? Flag for review Commit ECEF triple yes no pass fail

Specification: The Forward Geodetic-to-ECEF Definition

The conversion operates on a reference ellipsoid parameterized by semi-major axis aa (metres) and inverse flattening 1/f1/f. Flattening and the first eccentricity squared e2e^2 are derived constants, computed once rather than per coordinate:

f=11/f,e2=2f−f2 f = \frac{1}{1/f}, \qquad e^2 = 2f - f^2

Given geodetic latitude ϕ\phi (radians), longitude λ\lambda (radians), and ellipsoidal height hh (metres), the prime vertical radius of curvature NN — the distance along the ellipsoidal normal from the surface to the minor axis — resolves as:

N(ϕ)=a1−e2sin⁡2ϕ N(\phi) = \frac{a}{\sqrt{1 - e^2 \sin^2 \phi}}

The ECEF coordinates then follow directly, with no convergence loop:

X=(N+h)cos⁡ϕcos⁡λ X = (N + h)\cos\phi\cos\lambda
Y=(N+h)cos⁡ϕsin⁡λ Y = (N + h)\cos\phi\sin\lambda
Z=(N(1−e2)+h)sin⁡ϕ Z = \left(N(1 - e^2) + h\right)\sin\phi

Because the formulation is strictly forward, it carries none of the convergence uncertainty inherent in the inverse (ECEF-to-geodetic) direction, which is handled separately by Bowring’s closed-form solution in the verification section below. Determinism is preserved by three rules that the implementation enforces: evaluate every term in IEEE 754 double precision (float64), convert degrees to radians explicitly before any trigonometric call, and bound inputs to ϕ∈[−π/2,π/2]\phi \in [-\pi/2, \pi/2] and λ∈[−π,π]\lambda \in [-\pi, \pi]. One geodetic prerequisite governs the vertical: orthometric (levelled) heights must be converted to ellipsoidal heights through a geoid separation model before this step runs, or the resulting Z carries a systematic vertical bias of tens of metres. The reference-frame resolution that fixes which aa, 1/f1/f, and epoch apply is covered in setting up high-precision coordinate reference systems, and the broader compliance context lives in core transformation fundamentals and standards. Authoritative parameter definitions and transformation method codes trace to the ISO 19111 coordinate operation model.

Step-by-Step Implementation

The implementation builds in four ordered stages. Each stage is a runnable, type-hinted function in the same module; later blocks reuse the names defined earlier, so run them in sequence. Every function enforces numpy.float64 precision and avoids implicit float promotion.

Step 1 — Pin the ellipsoid constants

Compute ff and e2e^2 once from the published defining parameters. Re-deriving them inside the per-coordinate loop is a common source of low-order floating-point drift.

from __future__ import annotations

import logging
import math
from dataclasses import dataclass
from typing import Union

import numpy as np

logger = logging.getLogger("geodetic.ecef")

ArrayLike = Union[float, np.ndarray]


@dataclass(frozen=True)
class Ellipsoid:
    """Defining parameters of a reference ellipsoid (ISO 19111 datum ensemble member)."""

    a: float  # semi-major axis, metres
    inv_f: float  # inverse flattening, dimensionless

    @property
    def f(self) -> float:  # flattening: f = 1 / (1/f)
        return 1.0 / self.inv_f

    @property
    def e2(self) -> float:  # first eccentricity squared: e^2 = 2f - f^2 (ISO 19111 4.1.3)
        return 2.0 * self.f - self.f**2

    @property
    def b(self) -> float:  # semi-minor axis: b = a(1 - f)
        return self.a * (1.0 - self.f)


# EPSG:7030 — WGS 84 defining ellipsoid (a, 1/f). Used by EPSG:4979 / EPSG:4978.
WGS84 = Ellipsoid(a=6378137.0, inv_f=298.257223563)

Step 2 — Convert to radians and validate bounds

Latitude and longitude arrive in degrees from most data sources, so cast to float64 and convert before any trigonometric evaluation. Out-of-bounds or non-finite inputs are routed to a deterministic NaN sentinel rather than silently wrapping, satisfying the explicit-validity requirement of the ISO 19111 coordinate operation model.

def _to_radians_validated(
    lat_deg: ArrayLike, lon_deg: ArrayLike, h_m: ArrayLike
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Cast to float64, convert to radians, and NaN-route invalid coordinates."""
    # Explicit float64 cast prevents implicit float32 promotion downstream.
    lat = np.asarray(lat_deg, dtype=np.float64) * (math.pi / 180.0)
    lon = np.asarray(lon_deg, dtype=np.float64) * (math.pi / 180.0)
    h = np.asarray(h_m, dtype=np.float64)

    # Bounds per the geographic CRS domain of validity (phi in [-pi/2, pi/2]).
    valid = (
        (lat >= -math.pi / 2) & (lat <= math.pi / 2)
        & (lon >= -math.pi) & (lon <= math.pi)
        & np.isfinite(lat) & np.isfinite(lon) & np.isfinite(h)
    )
    if not np.all(valid):
        bad = ~valid
        logger.warning("Coordinate bounds violation on %d input(s); routing to NaN.", int(bad.sum()))
        lat = np.where(bad, np.nan, lat)
        lon = np.where(bad, np.nan, lon)
        h = np.where(bad, np.nan, h)
    return lat, lon, h

Step 3 — Evaluate the prime vertical radius of curvature

Compute N(ϕ)N(\phi) with an explicit guard against a zero denominator. The denominator only reaches zero for a non-physical eccentricity, but the guard keeps the function total and prevents a ZeroDivisionError from ever propagating into a cadastral pipeline.

def prime_vertical_radius(lat_rad: np.ndarray, ell: Ellipsoid = WGS84) -> np.ndarray:
    """N(phi) = a / sqrt(1 - e^2 sin^2 phi)  (ISO 19111 geographic-to-geocentric)."""
    sin2 = np.sin(lat_rad) ** 2
    denom = np.sqrt(1.0 - ell.e2 * sin2)
    # Replace an exact zero with the smallest positive float64 to stay total.
    denom = np.where(denom == 0.0, np.finfo(np.float64).tiny, denom)
    return ell.a / denom

Step 4 — Assemble the ECEF vector

Combine the validated angles, the height, and NN into the X, Y, Z triple. The function returns three float64 arrays so that no downstream consumer can re-introduce a narrower dtype.

def geodetic_to_ecef(
    lat_deg: ArrayLike,
    lon_deg: ArrayLike,
    h_m: ArrayLike,
    ell: Ellipsoid = WGS84,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Deterministic forward transform: geographic (EPSG:4979) -> geocentric (EPSG:4978).

    Method: EPSG:9602 Geographic/geocentric conversions. Enforces float64,
    explicit bounds checking, and NaN fallback routing.
    """
    lat, lon, h = _to_radians_validated(lat_deg, lon_deg, h_m)
    n = prime_vertical_radius(lat, ell)

    cos_phi, sin_phi = np.cos(lat), np.sin(lat)
    cos_lam, sin_lam = np.cos(lon), np.sin(lon)

    x = (n + h) * cos_phi * cos_lam
    y = (n + h) * cos_phi * sin_lam
    z = (n * (1.0 - ell.e2) + h) * sin_phi  # note the (1 - e^2) factor on Z only
    return x.astype(np.float64), y.astype(np.float64), z.astype(np.float64)

ISO 19111 requires explicit Coordinate Reference System identification and method registration on every operation. In practice each batch is annotated with the EPSG URNs urn:ogc:def:crs:EPSG::4979 (WGS 84 3D geographic) and urn:ogc:def:crs:EPSG::4978 (WGS 84 geocentric), the method code EPSG:9602, and the coordinate epoch. Where the ECEF output feeds a localized mapping frame, downstream alignment usually requires implementing affine transformations for local grids to absorb crustal deformation and legacy datum offsets.

Parameter and Return-Value Reference

Every key input and output, with the type, unit, valid range, and why it matters for survey-grade work. The table scrolls horizontally on narrow screens.

Name Type Units Valid range Cadastral significance
lat_deg float / np.ndarray degrees [−90,90][-90, 90] Geodetic (not geocentric) latitude; a swap with longitude relocates a parcel by hundreds of km
lon_deg float / np.ndarray degrees [−180,180][-180, 180] East-positive; sign error mirrors the position across the prime meridian
h_m float / np.ndarray metres unbounded (practical [−500,9000][-500, 9000]) Must be ellipsoidal height; an orthometric height injects a geoid-separation bias up to ~100 m
ell.a float metres 6378137.0 for WGS 84 Semi-major axis fixes the absolute scale of X, Y, Z
ell.inv_f float dimensionless 298.257223563 for WGS 84 Inverse flattening drives e2e^2 and the N(1−e2)N(1-e^2) term on Z
X np.ndarray (float64) metres ≈±6.4×106\approx \pm 6.4\times10^6 Geocentric X toward ϕ=0,λ=0\phi=0,\lambda=0
Y np.ndarray (float64) metres ≈±6.4×106\approx \pm 6.4\times10^6 Geocentric Y toward ϕ=0,λ=90∘E\phi=0,\lambda=90^\circ\text{E}
Z np.ndarray (float64) metres ≈±6.4×106\approx \pm 6.4\times10^6 Geocentric Z toward the north pole

Worked Example with Real-World Coordinates

Two checks pin the implementation to values that are exact by definition of the WGS 84 ellipsoid, then a realistic monument confirms behaviour off the axes. At the equator/prime-meridian intersection (ϕ=0,λ=0,h=0\phi=0,\lambda=0,h=0) the result must equal the semi-major axis aa; at the north pole (ϕ=90∘,h=0\phi=90^\circ,h=0) the Z component must equal the semi-minor axis b=a(1−f)≈6356752.314 mb = a(1-f) \approx 6356752.314\,\text{m}.

# Definitional control points — exact for the WGS 84 ellipsoid (EPSG:7030).
x0, y0, z0 = geodetic_to_ecef(0.0, 0.0, 0.0)
assert np.isclose(x0, WGS84.a, atol=1e-6), x0          # X = a at (0, 0, 0)
assert np.isclose(y0, 0.0, atol=1e-6) and np.isclose(z0, 0.0, atol=1e-6)

_, _, zp = geodetic_to_ecef(90.0, 0.0, 0.0)
assert np.isclose(zp, WGS84.b, atol=1e-6), zp          # Z = b at the pole

# Realistic monument: Boulder, CO area, WGS 84 geographic (EPSG:4979).
lat, lon, h = 40.01500000, -105.27055556, 1655.000
X, Y, Z = geodetic_to_ecef(lat, lon, h)
print(f"X={X:.3f}  Y={Y:.3f}  Z={Z:.3f}")
# X=-1288680.076  Y=-4720150.098  Z=4080325.441  (metres, EPSG:4978)

The printed triple is the geocentric position the rest of the survey-grade chain consumes: a datum shift, then a projection, then a residual gate against independent control. Treat the printed values as illustrative to three decimals; the verification step below is what actually proves correctness to sub-millimetre tolerance, independent of any hand-copied number.

Verification and Residual Analysis

The robust check for a forward conversion is a round trip: invert the ECEF vector back to geodetic coordinates with an independent algorithm and confirm the residual is below the survey-grade tolerance. Bowring’s closed-form inverse is used here so the round trip does not simply re-run the forward math in reverse.

def ecef_to_geodetic(
    x: ArrayLike, y: ArrayLike, z: ArrayLike, ell: Ellipsoid = WGS84
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Bowring closed-form inverse: geocentric (EPSG:4978) -> geographic (EPSG:4979)."""
    x = np.asarray(x, dtype=np.float64)
    y = np.asarray(y, dtype=np.float64)
    z = np.asarray(z, dtype=np.float64)

    e2, a, b = ell.e2, ell.a, ell.b
    ep2 = (a**2 - b**2) / b**2  # second eccentricity squared
    p = np.sqrt(x**2 + y**2)
    theta = np.arctan2(z * a, p * b)
    lat = np.arctan2(
        z + ep2 * b * np.sin(theta) ** 3,
        p - e2 * a * np.cos(theta) ** 3,
    )
    lon = np.arctan2(y, x)
    n = a / np.sqrt(1.0 - e2 * np.sin(lat) ** 2)
    h = p / np.cos(lat) - n
    return np.degrees(lat), np.degrees(lon), h


def assert_round_trip(
    lat_deg: float, lon_deg: float, h_m: float, tol_m: float = 1e-4
) -> dict[str, float]:
    """Forward then inverse; emit a structured audit record and enforce tolerance."""
    X, Y, Z = geodetic_to_ecef(lat_deg, lon_deg, h_m)
    lat2, lon2, h2 = ecef_to_geodetic(X, Y, Z)

    # Horizontal residual via local metres-per-degree at this latitude.
    m_per_deg_lat = 111_132.0
    m_per_deg_lon = 111_320.0 * math.cos(math.radians(lat_deg))
    d_north = float((lat2 - lat_deg) * m_per_deg_lat)
    d_east = float((lon2 - lon_deg) * m_per_deg_lon)
    d_up = float(h2 - h_m)
    residual_3d = math.sqrt(d_north**2 + d_east**2 + d_up**2)

    record = {
        "source_crs": "EPSG:4979",
        "target_crs": "EPSG:4978",
        "method": "EPSG:9602",
        "residual_3d_m": residual_3d,
        "tolerance_m": tol_m,
    }
    logger.info("ecef_roundtrip %s", record)
    assert residual_3d <= tol_m, f"round-trip residual {residual_3d:.2e} m exceeds {tol_m} m"
    return record


# Confirm the Boulder monument survives a forward/inverse round trip within 0.1 mm.
report = assert_round_trip(40.01500000, -105.27055556, 1655.000, tol_m=1e-4)
assert report["residual_3d_m"] < 1e-4

When this conversion is the front of a longer chain, residuals must be tracked across stages rather than only at the end. Systematic (non-zero-mean) residuals point to a wrong reference frame or a missing epoch and are formalized by error distribution modeling in Python; residuals from regional datum reconciliation are absorbed by polynomial shift algorithms for regional adjustments and resolved rigorously by a least squares adjustment for control networks. The control-point comparison itself, applied to a full datum alignment, is detailed in validating datum alignment with control points, and step-by-step isolation of model artefacts from genuine discrepancies is covered in debugging polynomial shift residuals in GIS.

Troubleshooting and Gotchas

The conversion is short, but a handful of specific mistakes turn an exact closed-form mapping into a silently biased one.

My Z values are tens of metres off but X and Y look fine — what happened?
Almost always an orthometric height was passed in where an ellipsoidal height is required. The Z term carries the height most directly, so a geoid–ellipsoid separation (the geoid undulation, typically -100 m to +85 m globally) appears mainly there. Convert levelled heights through a geoid model before calling the forward conversion, and record which geoid model was used in the audit metadata.
Why did re-deriving e² inside a vectorized loop change my results in the last digits?
Recomputing $e^2 = 2f - f^2$ per coordinate lets the compiler reorder the subtraction differently each call, perturbing the low-order bits of a value used in every term. Compute the ellipsoid constants once (as the frozen `Ellipsoid` dataclass does) so $e^2$ is bit-for-bit identical across the whole batch — a precondition for reproducible audit hashes.
Everything was correct in a notebook but the parcels shifted hundreds of kilometres in production.
This is the classic latitude/longitude axis-order swap. EPSG defines geographic CRSs as latitude-first, while many libraries and CSV exports are longitude-first. Assert the axis order from the source WKT2:2019 definition before calling `geodetic_to_ecef`, and never infer it from the magnitude of the numbers.
I mixed float32 arrays into the pipeline and lost sub-centimetre agreement.
State Plane and ECEF coordinates routinely exceed $10^6$ metres, and `float32` carries only ~7 significant digits — about a metre of resolution at that magnitude. The forward functions cast every input to `float64` at ingress and return `float64`; keep every downstream stage in double precision and difference against a local origin before accumulating sums.
A coordinate at the exact pole or with a NaN crashed the batch.
The implementation guards the $N$ denominator with `np.finfo(np.float64).tiny` and routes out-of-bounds or non-finite inputs to a NaN sentinel via `_to_radians_validated`, so a single bad row degrades to NaN instead of raising. Filter or flag NaN outputs explicitly before committing — never let them propagate silently into a deliverable.