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.

Contribution of each series term to the easting Bar chart of the contribution of each term of the transverse Mercator easting series at three hundred kilometres from the central meridian, in metres: the first term contributes 299 733, the third 265.8, the fifth 0.42 and the seventh 0.0003. Truncating after the fifth term leaves a sub-millimetre error; truncating after the third leaves nearly half a metre, which is why the term count is part of the specification and not an optimisation. 0.0001 0.001 0.01 0.1 1 10 100 1000 10000 100000 1e+06 m 299733 term 1 265.8 term 3 0.42 term 5 0.0003 term 7

Figure — how much each term of the series is worth at 300 km from the central meridian.

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

Constants a transverse Mercator forward projection requires Seven constants must be fixed before the series runs. The semi-major axis and flattening come from the ellipsoid of the source CRS, not from a global default. The central meridian and the latitude of origin fix the zone. The scale factor at the central meridian is 0.9996 for UTM but not for every national grid. The false easting and false northing keep coordinates positive. Taking any of these from the wrong grid definition produces coordinates that look entirely reasonable. a semi-major axis metres from the CRS ellipsoid f flattening ratio from the CRS ellipsoid lambda0 central meridian degrees fixes the zone phi0 latitude of origin degrees usually 0 k0 scale at CM ratio 0.9996 for UTM only FE / FN false origin metres keeps values positive

Figure — the seven constants a forward projection needs before it can start.

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.

Frequently Asked Questions

How many terms of the series do I need?

Through the eighth order — commonly written as terms one, three, five and seven of the easting series — is the standard for survey work, and it holds sub-millimetre accuracy across a normal three-degree zone. Truncating after the fifth term costs under a millimetre at the zone edge; truncating after the third costs nearly half a metre there while looking perfect near the central meridian, which is the worst kind of error because testing near the origin will not catch it.

Why does my implementation disagree with PROJ by a few millimetres far from the zone?

Almost certainly the series order, or the meridian arc formula. PROJ uses an expansion accurate well beyond the standard zone width; a classical eighth-order series diverges from it as you move out. Agreement should be tested at the zone edge and beyond, not at the origin. If the disagreement is constant rather than growing, suspect the false easting or the scale factor instead.

Should the projection routine also return the scale factor and convergence?

Yes. They are computed from quantities the projection already has, so they are nearly free, and without them the easting and northing cannot be used for reductions or bearings. A routine that returns only the coordinate pair pushes that work — and the chance of getting it wrong — onto every caller.