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.

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 $E$ , the northing $N$ , and the point scale factor $k$ that converts a grid distance into a true geodesic distance on the ellipsoid.

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

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 = e^2/(1 - e^2)$ is the second eccentricity squared. The meridional arc length $M$ from the equator to latitude $\phi$ is

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 $M_0$ the meridional arc at the latitude of origin $\phi_0$ and $\mathit{FE}$ , $\mathit{FN}$ the false easting and northing:

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]

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 = 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 $h$ and the geometric-mean radius of curvature $R = \sqrt{M_\rho N}$ — where $M_\rho = a(1-e^2)/(1-e^2\sin^2\phi)^{3/2}$ is the meridian radius — the combined grid-to-ground factor is

k_{\text{combined}} = k \times \frac{R}{R + h}

Retaining the $A^5$ / $A^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.

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°, $k_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.999964787$ ) barely differs from $k_0$ . The elevation reduction does the heavier lifting: at 350 m it pulls the combined factor down to $0.999909912$ , roughly 90 ppm. A 1 000 m ground distance measured at this monument therefore corresponds to $1000 \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 $a$ 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 $A^5$ / $A^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 $k$ reduces to the ellipsoid only; it does not account for terrain height. Applying $k$ instead of $k_{\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 ( $h$ ), not the geoid (orthometric $H$ ). Feeding a levelling-derived orthometric height skips the geoid undulation $N$ , 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.

Core Transformation Fundamentals & Standards — parent reference for the full ISO 19111 transformation framework.
Handling CRS mismatches in cadastral datasets — resolving legacy state-plane and UTM realizations before projection.
NADCON vs NTv2: choosing the right datum shift — selecting the correction kernel that precedes the projection step.
Setting up high-precision coordinate reference systems — pinning down axis order, units, and zone parameters.
Validating datum alignment with control points — network-level residual analysis and RMSE reporting.
Geodetic conversion math: ellipsoid to Cartesian — the companion forward/inverse geodetic transforms.

Projection Math Fundamentals for Cadastral Surveys

The Conformal Projection Specification #

Step-by-Step Implementation #

Step 1 — Declare the zone and ellipsoid explicitly #

Step 2 — Compute the meridional arc #

Step 3 — Forward Transverse Mercator projection #

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

Parameter and Return-Value Reference #

Worked Example: a Vermont State Plane Monument #

Verification and Residual Analysis #

Troubleshooting and Gotchas #

Related #