CoordTransform Cadastral & High-Precision Python

Projection Math Fundamentals for Cadastral Surveys

Turning ellipsoidal survey observations into legally defensible plane coordinates is the projection-math sub-task at the heart of the Core Transformation Fundamentals & Standards reference: a parcel corner, a monument, or a boundary bearing only becomes a record-of-survey coordinate once a conformal map projection has mapped it from latitude and longitude onto a grid. This page narrows the standards framework to one operation — computing Transverse Mercator eastings and northings, the point scale factor, and the combined grid-to-ground factor to survey-grade precision — and shows the deterministic arithmetic that keeps that operation auditable under ISO 19111 coordinate-operation semantics. Unlike thematic mapping or web visualization, cadastral projection cannot tolerate truncated series, implicit ellipsoid defaults, or floating-point drift: a one-part-per-million scale error is roughly 1 mm per kilometre, which compounds into a metre across a municipal dataset and into a contested boundary in court.

Cadastral projection-math workflow Vertical flowchart. Geodetic input latitude, longitude and ellipsoidal height feeds a step that resolves the EPSG zone into explicit ellipsoid and conversion parameters a, one-over-f, k-zero and the central meridian. The meridional arcs M and M-zero are computed, then the Transverse Mercator forward series yields easting, northing and the point scale factor k. An elevation reduction R over R-plus-h is computed and multiplied with the grid scale to give the combined grid-to-ground factor. A tolerance gate tests whether the residual against published control is at or below twenty-five millimetres: if it passes the coordinate is committed with round-half-to-even rounding, otherwise it is rejected to manual review. Geodetic input φ, λ, h deg, deg, m · target datum Resolve EPSG zone (explicit) a, 1∕f, k₀, λ₀, FE, FN Meridional arc M, M₀ · e⁶ series kept TM forward series (A⁵∕A⁶) → E, N, point scale k Elevation reduction R ∕ (R + h) Combined grid-to-ground k_comb = k · R∕(R+h) residual ≤ 25 mm? Reject → manual review Commit E, N round-half-to-even pass fail

The arithmetic below is the foundation that the grid-shift and datum stages build on. Before a coordinate is projected it must already sit in the correct datum realization — the choice of correction method is covered in NADCON vs NTv2: choosing the right datum shift, and the explicit reference-frame resolution that pins down axis order and units is covered in setting up high-precision coordinate reference systems.

The Conformal Projection Specification

Transverse Mercator ™ and Lambert Conformal Conic dominate cadastral practice because they are conformal: they preserve local angles, so a measured bearing maps to a grid bearing with only a small, computable convergence correction. Conformality is purchased at the cost of scale distortion that grows with distance from the central meridian ™ or standard parallels (Lambert). For boundary retracement the surveyor therefore needs three quantities at every point, not just the coordinate pair: the easting EE, the northing NN, and the point scale factor kk that converts a grid distance into a true geodesic distance on the ellipsoid.

Point-scale distortion across a transverse Mercator zone Distortion in parts per million against distance from the central meridian, computed from the transverse Mercator point-scale series with a central-meridian scale factor of 0.9996. At the central meridian the grid is 400 parts per million short. The distortion passes through zero near 180 kilometres and reaches roughly plus 570 parts per million at the 350 kilometre zone edge, which is 0.57 metres in every kilometre of measured distance. -500 0 500 1000 0 50 100 150 200 250 300 350 distance from central meridian (km) ppm point scale (ppm)

Figure — point-scale distortion computed from the transverse Mercator series, k0 = 0.9996.

The exact TM forward mapping is a series in the ellipsoidal eccentricity. Writing ϕ\phi for latitude, λ\lambda for longitude, λ0\lambda_0 for the central meridian, k0k_0 for the central-meridian scale, aa for the semi-major axis and e2e^2 for the first eccentricity squared, the standard Snyder formulation introduces the auxiliary terms

N=a1−e2sin⁡2ϕ,T=tan⁡2ϕ,C=e′2cos⁡2ϕ,A=(λ−λ0)cos⁡ϕN = \frac{a}{\sqrt{1 - e^2\sin^2\phi}}, \quad T = \tan^2\phi, \quad C = e'^2\cos^2\phi, \quad A = (\lambda - \lambda_0)\cos\phi

where e′2=e2/(1−e2)e'^2 = e^2/(1 - e^2) is the second eccentricity squared. The meridional arc length MM from the equator to latitude ϕ\phi is

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

The plane coordinates then follow, with M0M_0 the meridional arc at the latitude of origin ϕ0\phi_0 and FE\mathit{FE}, FN\mathit{FN} the false easting and northing:

E=FE+k0N[A+(1−T+C)A36+(5−18T+T2+72C−58e′2)A5120]E = \mathit{FE} + k_0 N\left[A + (1 - T + C)\tfrac{A^3}{6} + (5 - 18T + T^2 + 72C - 58e'^2)\tfrac{A^5}{120}\right]
Nplane=FN+k0{M−M0+Ntan⁡ϕ[A22+(5−T+9C+4C2)A424+(61−58T+T2+600C−330e′2)A6720]}N_{\text{plane}} = \mathit{FN} + k_0\left\{M - M_0 + N\tan\phi\left[\tfrac{A^2}{2} + (5 - T + 9C + 4C^2)\tfrac{A^4}{24} + (61 - 58T + T^2 + 600C - 330e'^2)\tfrac{A^6}{720}\right]\right\}

The point scale factor is the companion series

k=k0[1+(1+C)A22+(5−4T+42C+13C2−28e′2)A424+(61−148T+16T2)A6720]k = k_0\left[1 + (1 + C)\tfrac{A^2}{2} + (5 - 4T + 42C + 13C^2 - 28e'^2)\tfrac{A^4}{24} + (61 - 148T + 16T^2)\tfrac{A^6}{720}\right]

Finally, survey deliverables work at ground distance, so the grid scale must be combined with an elevation (sea-level) reduction. With ellipsoidal height hh and the geometric-mean radius of curvature R=MρNR = \sqrt{M_\rho N} — where Mρ=a(1−e2)/(1−e2sin⁡2ϕ)3/2M_\rho = a(1-e^2)/(1-e^2\sin^2\phi)^{3/2} is the meridian radius — the combined grid-to-ground factor is

kcombined=k×RR+hk_{\text{combined}} = k \times \frac{R}{R + h}

Retaining the A5A^5/A6A^6 terms keeps the projection error below 1 mm out to roughly 5° of longitude from the central meridian, which comfortably covers any state-plane or UTM zone. Dropping them — the truncation most “good enough” GIS code commits — reintroduces decimetre-level error at the zone edge.

Step-by-Step Implementation

The reference implementation keeps every ellipsoid and zone parameter explicit, performs the trigonometric series in IEEE 754 float64, and applies deterministic round-half-to-even only at the final commit so that no intermediate rounding biases the result. There are no hidden library defaults; the zone is declared in full, satisfying the ISO 19111 requirement that a coordinate operation name its source CRS, target CRS, and every method parameter.

What a cadastral projection call must return alongside the coordinates A projection call fans out into four outputs that must travel together. The easting and northing are the coordinates themselves. The point scale factor converts between ground distance and grid distance. The grid convergence converts between geodetic azimuth and grid bearing. The zone identifier and its EPSG code make the pair interpretable at all. Dropping any of the last three turns a defensible observation into an ambiguous number pair. project(lat, lon) one geodetic position Easting, northing metres on the grid Point scale k ground to grid Convergence gamma azimuth to bearing Zone + EPSG code makes it interpretable

Figure — the four quantities a cadastral projection routine must return together.

Step 1 — Declare the zone and ellipsoid explicitly 

from __future__ import annotations

import math
from dataclasses import dataclass


@dataclass(frozen=True)
class TMZone:
    """A fully specified Transverse Mercator zone (ISO 19111 conversion params).

    No implicit defaults: every parameter that affects the coordinate is named,
    so a pipeline can hash this object into its audit record.
    """

    a: float          # semi-major axis, metres (e.g. GRS80 = 6378137.0)
    inv_f: float      # inverse flattening 1/f (GRS80 = 298.257222101)
    lat0_deg: float   # latitude of origin, degrees
    lon0_deg: float   # central meridian, degrees
    k0: float         # scale factor on the central meridian
    false_easting: float
    false_northing: float

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

    @property
    def e2(self) -> float:
        # First eccentricity squared, derived — never assumed.
        return 2.0 * self.f - self.f * self.f

    @property
    def ep2(self) -> float:
        # Second eccentricity squared e'^2 = e^2 / (1 - e^2).
        return self.e2 / (1.0 - self.e2)

Step 2 — Compute the meridional arc 

def meridian_arc(a: float, e2: float, lat_rad: float) -> float:
    """Meridional arc length M from the equator to lat_rad (Snyder 3-21).

    The e^6 term is retained so the arc is accurate to sub-millimetre at any
    latitude — required for survey-grade northings.
    """
    e4 = e2 * e2
    e6 = e4 * e2
    return a * (
        (1 - e2 / 4 - 3 * e4 / 64 - 5 * e6 / 256) * lat_rad
        - (3 * e2 / 8 + 3 * e4 / 32 + 45 * e6 / 1024) * math.sin(2 * lat_rad)
        + (15 * e4 / 256 + 45 * e6 / 1024) * math.sin(4 * lat_rad)
        - (35 * e6 / 3072) * math.sin(6 * lat_rad)
    )

Step 3 — Forward Transverse Mercator projection 

def tm_forward(zone: TMZone, lat_deg: float, lon_deg: float) -> tuple[float, float, float]:
    """Project geodetic (lat, lon) to grid (E, N) and the point scale factor k.

    Returns (easting_m, northing_m, k_point). Implements Snyder 8-9..8-11 with
    the A^5/A^6 terms kept (sub-millimetre to ~5 deg from the central meridian).
    """
    lat = math.radians(lat_deg)
    lon = math.radians(lon_deg)
    lat0 = math.radians(zone.lat0_deg)
    lon0 = math.radians(zone.lon0_deg)

    e2 = zone.e2
    ep2 = zone.ep2
    sin_lat = math.sin(lat)
    cos_lat = math.cos(lat)
    tan_lat = math.tan(lat)

    nu = zone.a / math.sqrt(1.0 - e2 * sin_lat * sin_lat)  # prime-vertical radius
    t = tan_lat * tan_lat                                   # T = tan^2(phi)
    c = ep2 * cos_lat * cos_lat                             # C = e'^2 cos^2(phi)
    a_term = (lon - lon0) * cos_lat                         # A = dlon * cos(phi)

    m = meridian_arc(zone.a, e2, lat)
    m0 = meridian_arc(zone.a, e2, lat0)

    easting = zone.false_easting + zone.k0 * nu * (
        a_term
        + (1 - t + c) * a_term ** 3 / 6
        + (5 - 18 * t + t * t + 72 * c - 58 * ep2) * a_term ** 5 / 120
    )
    northing = zone.false_northing + zone.k0 * (
        m - m0
        + nu * tan_lat * (
            a_term ** 2 / 2
            + (5 - t + 9 * c + 4 * c * c) * a_term ** 4 / 24
            + (61 - 58 * t + t * t + 600 * c - 330 * ep2) * a_term ** 6 / 720
        )
    )
    k_point = zone.k0 * (
        1
        + (1 + c) * a_term ** 2 / 2
        + (5 - 4 * t + 42 * c + 13 * c * c - 28 * ep2) * a_term ** 4 / 24
        + (61 - 148 * t + 16 * t * t) * a_term ** 6 / 720
    )
    return easting, northing, k_point

Step 4 — Resolve the combined grid-to-ground scale factor 

def combined_scale_factor(
    zone: TMZone, lat_deg: float, ellipsoidal_height_m: float, k_point: float
) -> tuple[float, float]:
    """Combine the projection (grid) factor with the elevation reduction.

    Returns (k_combined, R) where R is the geometric-mean radius of curvature.
    k_combined converts a ground distance to a grid distance and back.
    """
    lat = math.radians(lat_deg)
    e2 = zone.e2
    sin2 = math.sin(lat) ** 2
    meridian_radius = zone.a * (1.0 - e2) / (1.0 - e2 * sin2) ** 1.5
    prime_vertical = zone.a / math.sqrt(1.0 - e2 * sin2)
    r_mean = math.sqrt(meridian_radius * prime_vertical)        # R = sqrt(M_rho * nu)
    k_elev = r_mean / (r_mean + ellipsoidal_height_m)           # sea-level reduction
    return k_point * k_elev, r_mean

Parameter and Return-Value Reference

Symbol / field Type Units Valid range Cadastral significance
a float metres 6 370 000–6 380 000 Ellipsoid semi-major axis; GRS80 = 6378137.0 for NAD83 work
inv_f float dimensionless 290–300 Inverse flattening; fixes ellipsoid shape (GRS80 = 298.257222101)
lat0_deg, lon0_deg float degrees −90…90 / −180…180 Origin and central meridian of the published zone
k0 float dimensionless 0.9996–1.0 Central-meridian scale; UTM = 0.9996, varies per SPCS zone
false_easting, false_northing float metres ≥ 0 Offsets that keep grid coordinates positive
lat_deg, lon_deg float degrees within zone extent The geodetic point to be projected
ellipsoidal_height_m float metres −500…9000 Height above the ellipsoid; drives the elevation reduction
easting_m, northing_m float metres zone-dependent Plane coordinates committed to the record of survey
k_point float dimensionless ≈ k0…1.0004 Grid scale; multiply geodesic distance to get grid distance
k_combined float dimensionless ≈ 0.9994…1.0004 Grid-to-ground factor; the value the field crew actually applies
r_mean (R) float metres ≈ 6 370 000–6 400 000 Geometric-mean radius used in the elevation reduction

Worked Example: a Vermont State Plane Monument

Vermont uses a single Transverse Mercator zone, EPSG:32145 (NAD83 / Vermont, GRS80), with central meridian −72.5°, latitude of origin 42.5°, k0=0.999964285714286k_0 = 0.999964285714286, false easting 500 000 m and false northing 0. Take a control point near Montpelier at latitude 44.26°, longitude −72.58°, ellipsoidal height 350 m.

vermont = TMZone(
    a=6378137.0,
    inv_f=298.257222101,
    lat0_deg=42.5,
    lon0_deg=-72.5,
    k0=0.999964285714286,
    false_easting=500000.0,
    false_northing=0.0,
)

easting, northing, k_point = tm_forward(vermont, 44.26, -72.58)
k_combined, r_mean = combined_scale_factor(vermont, 44.26, 350.0, k_point)

print(f"E = {easting:.4f} m")          # E = 493611.8194 m
print(f"N = {northing:.4f} m")         # N = 195532.3917 m
print(f"k_point   = {k_point:.9f}")    # k_point   = 0.999964787
print(f"k_combined = {k_combined:.9f}")# k_combined = 0.999909912

The point sits about 6.4 km west of the central meridian, so the grid scale (0.9999647870.999964787) barely differs from k0k_0. The elevation reduction does the heavier lifting: at 350 m it pulls the combined factor down to 0.9999099120.999909912, roughly 90 ppm. A 1 000 m ground distance measured at this monument therefore corresponds to 1000×0.999909912=999.9101000 \times 0.999909912 = 999.910 m on the grid — a 90 mm difference that a cadastral pipeline must apply, not ignore.

Verification and Residual Analysis

A projected coordinate is only trustworthy once it has been compared against an independently published value for the same monument. The verification step computes the horizontal residual, tests it against the survey-grade tolerance, and emits a structured record that an auditor can replay. Tolerance budgeting itself is treated in depth in tuning transformation thresholds for survey-grade work; here it is the gate that decides commit versus reject.

import json
from dataclasses import dataclass


@dataclass(frozen=True)
class Residual:
    d_east_m: float
    d_north_m: float
    horizontal_m: float
    within_tolerance: bool


def verify_against_control(
    computed_en: tuple[float, float],
    published_en: tuple[float, float],
    tolerance_m: float = 0.025,  # +/-25 mm, 95% CI: typical urban cadastral budget
) -> Residual:
    """Compare a computed (E, N) against a published control coordinate."""
    d_e = computed_en[0] - published_en[0]
    d_n = computed_en[1] - published_en[1]
    horizontal = math.hypot(d_e, d_n)
    return Residual(d_e, d_n, horizontal, horizontal <= tolerance_m)


def audit_record(epsg: int, residual: Residual) -> str:
    """Emit a deterministic JSON log line for the audit trail."""
    return json.dumps(
        {
            "epsg": epsg,
            "d_east_mm": round(residual.d_east_m * 1000, 2),
            "d_north_mm": round(residual.d_north_m * 1000, 2),
            "horizontal_mm": round(residual.horizontal_m * 1000, 2),
            "status": "PASS" if residual.within_tolerance else "FAIL",
        },
        sort_keys=True,
    )


# Published monument coordinate for the worked example above:
published = (493611.804, 195532.402)
res = verify_against_control((493611.8194, 195532.3917), published)
print(audit_record(32145, res))
# {"d_east_mm": 15.4, "d_north_mm": -10.3, "epsg": 32145, "horizontal_mm": 18.53, "status": "PASS"}

assert res.within_tolerance, "Residual exceeds survey-grade tolerance — reject and review"

The 18.5 mm horizontal residual clears the ±25 mm gate, so the coordinate commits with a PASS record. A failure here is never silently absorbed; it routes to manual review, exactly as a missing grid routes through the fallback routing strategies for missing grid files. Aggregating these residuals across a control network — RMSE, residual histograms, and outlier flagging — is the subject of validating datum alignment with control points.

Troubleshooting and Gotchas

  • Wrong ellipsoid for the datum. Projecting NAD83 coordinates with WGS84 parameters (or vice versa) looks harmless because aa differs by only ~0.1 mm, but pairing the wrong ellipsoid with the wrong datum realization injects up to a 2 m horizontal shift. Always derive e2 from the zone’s declared inv_f; never hard-code an eccentricity.
  • Truncating the projection series. Many lightweight projection routines drop the A5A^5/A6A^6 terms. Inside a few kilometres of the central meridian the error is invisible, but at the longitude edge of a state-plane or UTM zone it reaches decimetres — past any cadastral tolerance. Keep the full series shown in Step 3.
  • Confusing grid scale with combined scale. The point scale factor kk reduces to the ellipsoid only; it does not account for terrain height. Applying kk instead of kcombinedk_{\text{combined}} leaves the entire elevation reduction (90 ppm at 350 m, far more in the mountains) unmodelled. Field distances always use the combined factor.
  • Orthometric height vs ellipsoidal height. The elevation reduction needs height above the ellipsoid (hh), not the geoid (orthometric HH). Feeding a levelling-derived orthometric height skips the geoid undulation NN, biasing the combined factor by the local separation — often 20–30 m of height, i.e. several ppm of scale.
  • Projecting before the datum shift. Projection assumes the coordinate is already in the target datum. If the input is still in a legacy realization, project after the shift, not before; reconciling those realizations is covered in handling CRS mismatches in cadastral datasets.

Frequently Asked Questions

Why does a grid distance differ from the distance I measured on the ground?

Because a projected coordinate system is a flattened model of a curved surface, and flattening always distorts. The point scale factor tells you by how much at a given position: multiply a ground distance by the mean point scale along the line to get the grid distance, and divide to go back. On a transverse Mercator zone with a central-meridian scale factor of 0.9996 the grid is 400 parts per million short on the central meridian and roughly 570 parts per million long at the zone edge — a 40 centimetre disagreement over a kilometre in either direction. A cadastral traverse that reduces its distances without applying the scale factor will close beautifully against itself and disagree with every neighbouring survey.

Is the grid convergence the same as magnetic declination?

No, and conflating them is a classic source of rotated parcels. Grid convergence is the angle between grid north and geodetic (true) north at a point, and it is a purely geometric consequence of the projection — it depends only on your position relative to the central meridian and can be computed exactly. Magnetic declination is the angle between magnetic north and true north, is a property of the Earth’s magnetic field, changes measurably from year to year, and can only be looked up in a model. Converting a geodetic azimuth to a grid bearing needs convergence; converting a compass bearing needs both.

Can I just use a local plane approximation over a small site?

Over a site of a few hundred metres the curvature error is genuinely negligible, and a local plane is defensible for internal geometry. It stops being defensible the moment the coordinates leave the site: a local plane has no datum, no scale factor and no convergence, so nobody downstream can relate it to anything. If you work in a local plane, publish the transformation that maps it to the national grid alongside it — that is exactly the job the affine and similarity transformations in the workflows section exist to do.

Which ellipsoid constants should a projection routine use?

The ones belonging to the geographic CRS that the coordinates are actually in, read from the CRS definition rather than from a constant in your source file. GRS80 and WGS84 differ in the flattening term at about the tenth decimal place, which is under a millimetre and usually ignorable; Clarke 1866 and GRS80 differ by hundreds of metres in the semi-major axis and are absolutely not interchangeable. A projection routine that hard-codes its ellipsoid works perfectly until the day someone hands it NAD27 data.