CoordTransform Cadastral & High-Precision Python

Implementing Transverse Mercator Forward Projection

The forward Transverse Mercator mapping — carrying geodetic (φ,λ)(\varphi, \lambda) to grid easting and northing — is the single most-used projection in cadastral work, and this page derives one self-contained, float64 implementation of it using the classic Redfearn meridional-arc series, the direction-of-travel companion to the round-trip treatment in map projection forward and inverse in Python under Algorithmic Math & Geodetic Workflows. Where that reference builds the modern Krüger nn-series and its inverse, here the forward mapping is expanded explicitly in powers of the longitude difference ω=λ−λ0\omega = \lambda - \lambda_0 so every term — meridional arc, prime-vertical radius, the curvature ratio — is visible and separately testable.

Redfearn forward Transverse Mercator data flow Geodetic latitude and longitude feed a box that computes the meridional arc M of phi and the prime-vertical radius nu, then a series in powers of omega equal to lambda minus the central meridian. The northing branch adds k0 times M minus M0 plus the even-power terms; the easting branch is k0 times nu times the odd-power terms plus the false easting. Output is easting and northing in metres. φ, λ (degrees) ω = λ − λ₀ Redfearn series M(φ), ν, ψ = ν/ρ powers of ω N = FN + k₀(…) E = FE + k₀ν(…)

Figure — the Redfearn forward mapping splits into an even-power northing branch built on the meridional arc and an odd-power easting branch.

The Series, Term by Term

Transverse Mercator is conformal, so a single scale governs both axes at each point, but that scale grows with distance from the central meridian; the Redfearn expansion makes the growth explicit. Two ellipsoidal radii appear: the prime-vertical radius ν\nu and the meridional radius ρ\rho, whose ratio ψ=ν/ρ=1+e′2cos⁡2φ\psi = \nu/\rho = 1 + e'^2\cos^2\varphi controls the higher-order terms.

ν=a1−e2sin⁡2φ,e′2=e21−e2,ω=λ−λ0\nu = \frac{a}{\sqrt{1 - e^2\sin^2\varphi}}, \qquad e'^2 = \frac{e^2}{1-e^2}, \qquad \omega = \lambda - \lambda_0

The northing is anchored on the meridional arc M(φ)M(\varphi), the true distance along the ellipsoid meridian from the equator to latitude φ\varphi:

M(φ)=a[(1−e24−3e464−5e6256)φ−(3e28+3e432+45e61024)sin⁡2φ+(15e4256+45e61024)sin⁡4φ−35e63072sin⁡6φ]M(\varphi) = a\Big[\big(1 - \tfrac{e^2}{4} - \tfrac{3e^4}{64} - \tfrac{5e^6}{256}\big)\varphi - \big(\tfrac{3e^2}{8} + \tfrac{3e^4}{32} + \tfrac{45e^6}{1024}\big)\sin 2\varphi + \big(\tfrac{15e^4}{256} + \tfrac{45e^6}{1024}\big)\sin 4\varphi - \tfrac{35e^6}{3072}\sin 6\varphi\Big]

With t=tan⁡φt = \tan\varphi, the easting collects the odd powers of ω\omega and the northing the even powers, each scaled by k0k_0:

E=FE+k0 ν[ωcos⁡φ+ω36cos⁡3φ (ψ−t2)+ω5120cos⁡5φ (…)]E = \mathrm{FE} + k_0\,\nu\Big[\omega\cos\varphi + \tfrac{\omega^3}{6}\cos^3\varphi\,(\psi - t^2) + \tfrac{\omega^5}{120}\cos^5\varphi\,(\ldots)\Big]
N=FN+k0[M(φ)−M(φ0)+ν(ω22sin⁡φcos⁡φ+ω424sin⁡φcos⁡3φ (4ψ2+ψ−t2)+…)]N = \mathrm{FN} + k_0\Big[M(\varphi) - M(\varphi_0) + \nu\big(\tfrac{\omega^2}{2}\sin\varphi\cos\varphi + \tfrac{\omega^4}{24}\sin\varphi\cos^3\varphi\,(4\psi^2 + \psi - t^2) + \ldots\big)\Big]

Truncating after the fifth/sixth-order terms holds sub-millimetre agreement across a standard 6° zone — adequate for cadastral use, though the Krüger nn-series in the parent reference is marginally tighter near the zone edge.

Complete Runnable Implementation

The function below is self-contained: it takes the ellipsoid constants, the projection parameters, and a point, and returns easting, northing, and the point scale factor. It works entirely in float64 and converts degrees to radians internally so the caller never has to.

from __future__ import annotations
import math
from dataclasses import dataclass


@dataclass(frozen=True)
class TMParameters:
    """Transverse Mercator definition. Angles in degrees, lengths in metres."""
    a: float                # semi-major axis, e.g. 6378137.0 (WGS84)
    inv_f: float            # inverse flattening, e.g. 298.257223563
    lon_cm_deg: float       # central meridian (longitude of natural origin)
    lat_origin_deg: float   # latitude of natural origin
    k0: float               # scale factor on the central meridian
    false_easting: float    # metres
    false_northing: float   # metres


def _meridian_arc(phi: float, a: float, e2: float) -> float:
    """Meridional arc M(phi): ellipsoid distance from the equator to phi (radians)."""
    return a * (
        (1 - e2/4 - 3*e2**2/64 - 5*e2**3/256) * phi
        - (3*e2/8 + 3*e2**2/32 + 45*e2**3/1024) * math.sin(2*phi)
        + (15*e2**2/256 + 45*e2**3/1024) * math.sin(4*phi)
        - (35*e2**3/3072) * math.sin(6*phi)
    )


def tm_forward(p: TMParameters, lat_deg: float, lon_deg: float) -> tuple[float, float, float]:
    """Project geodetic (lat, lon) in degrees to (easting, northing, scale_factor).

    Ellipsoidal Redfearn series truncated at the 5th/6th order in the longitude
    difference; float64 throughout. Returns metres and a dimensionless scale factor.
    """
    f = 1.0 / p.inv_f
    e2 = f * (2.0 - f)                       # first eccentricity squared
    ep2 = e2 / (1.0 - e2)                    # second eccentricity squared

    phi = math.radians(lat_deg)
    phi0 = math.radians(p.lat_origin_deg)
    omega = math.radians(lon_deg - p.lon_cm_deg)   # degrees -> radians here only

    sin_phi, cos_phi = math.sin(phi), math.cos(phi)
    t = math.tan(phi)
    t2, t4 = t * t, t**4
    nu = p.a / math.sqrt(1.0 - e2 * sin_phi * sin_phi)   # prime-vertical radius
    psi = 1.0 + ep2 * cos_phi * cos_phi                  # nu / rho

    c = cos_phi
    easting = p.false_easting + p.k0 * nu * (
        omega * c
        + (omega**3 / 6.0) * c**3 * (psi - t2)
        + (omega**5 / 120.0) * c**5 * (
            4.0*psi**3*(1 - 6*t2) + psi**2*(1 + 8*t2) - psi*2*t2 + t4)
    )

    m_phi = _meridian_arc(phi, p.a, e2)
    m_phi0 = _meridian_arc(phi0, p.a, e2)
    northing = p.false_northing + p.k0 * (
        m_phi - m_phi0
        + nu * (
            (omega**2 / 2.0) * sin_phi * c
            + (omega**4 / 24.0) * sin_phi * c**3 * (4.0*psi**2 + psi - t2)
            + (omega**6 / 720.0) * sin_phi * c**5 * (
                8.0*psi**4*(11 - 24*t2) - 28*psi**3*(1 - 6*t2)
                + psi**2*(1 - 32*t2) - psi*2*t2 + t4)
        )
    )

    # Point scale factor: k0 grown by the distance from the central meridian.
    k = p.k0 * (
        1.0
        + 0.5 * omega**2 * c**2 * psi
        + (omega**4 / 24.0) * c**4 * (4.0*psi**3*(1 - 6*t2) + psi**2*(1 + 24*t2) - 4*psi*t2)
    )
    return easting, northing, k

Parameter Reference

Name Type Units Valid range Notes
a float metres 6378137.0 (WGS84/GRS80) semi-major axis of the CRS ellipsoid
inv_f float — 298.2572… inverse flattening; sets the eccentricity
lon_cm_deg float degrees zone central meridian subtract from lon before the series; wrong value shifts every easting
lat_origin_deg float degrees 0 for UTM enters only through M(φ0)M(\varphi_0) in the northing
k0 float — 0.9996 (UTM) scale factor on the central meridian
false_easting / false_northing float metres 500000 / 0 (UTM N) added last to keep coordinates positive
lat_deg, lon_deg float degrees within the zone the point to project; degrees, not radians
return tuple[float, float, float] m, m, — — easting, northing, point scale factor

Worked Example

Projecting a point in Vienna onto UTM zone 33N (EPSG:32633: central meridian 15°E, k0=0.9996k_0 = 0.9996, false easting 500 000 m) reproduces the published grid coordinate.

utm33 = TMParameters(
    a=6378137.0, inv_f=298.257223563, lon_cm_deg=15.0, lat_origin_deg=0.0,
    k0=0.9996, false_easting=500000.0, false_northing=0.0,
)
E, N, k = tm_forward(utm33, 48.208889, 16.372778)
print(round(E, 3), round(N, 3), round(k, 9))
# 601987.912 5340428.817 0.999727803

Validation Check

The forward output is gated against an independent oracle — pyproj, driving the same EPSG:32633 definition — and accepted only if the two agree within the millimetre budget cadastral deliverables are graded against.

from pyproj import Transformer

xf = Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)
px, py = xf.transform(16.372778, 48.208889)     # always_xy => (lon, lat)
E, N, _ = tm_forward(utm33, 48.208889, 16.372778)

separation = math.hypot(E - px, N - py)
assert separation < 1e-3, f"Redfearn forward disagrees with PROJ by {separation*1e3:.3f} mm"

The recovery of (φ,λ)(\varphi, \lambda) from these eastings and northings — the inverse that must round-trip this result — is built in the parent map projection forward and inverse in Python, and the covariance handling that turns a projection residual into a reported uncertainty is developed in error distribution modeling in Python.

Common Mistakes

Using the spherical Transverse Mercator formulas for cadastral accuracy
Dropping the eccentricity terms — treating the Earth as a sphere — introduces errors of hundreds of metres at mid-latitudes because it ignores the meridional arc's ellipsoidal shape. The _meridian_arc series and the psi = nu/rho ratio are exactly the ellipsoidal corrections a survey grid depends on; never substitute R * phi for M(phi).
Forgetting the false easting and northing
The series computes the signed distance from the central meridian and the latitude of origin; those are added to false_easting and false_northing only at the very end. Omit them and points west of the central meridian come out negative, and every coordinate is offset by 500 000 m in easting for a UTM zone — a defect that no downstream check will silently repair.
Wrong central meridian, or feeding degrees where radians are expected
Two failure modes that look similar. Using a neighbouring zone's central meridian offsets every easting by the meridian difference; passing lon_deg already in radians (or forgetting to subtract lon_cm_deg before converting) inflates omega and produces nonsense. This implementation converts degrees to radians in exactly one place — inside tm_forward — so callers must always pass plain degrees.