Map Projection Forward and Inverse in Python

Every cadastral survey plane is a projected coordinate reference system, and every parcel corner written to it has passed through a map projection: the forward mapping that carries a geographic coordinate $(\varphi, \lambda)$ on the ellipsoid onto a planar easting and northing $(E, N)$ , and the inverse mapping that recovers the geographic coordinate from the grid. This page narrows Algorithmic Math & Geodetic Workflows to that single conformal operation, implemented so both directions round-trip to within a survey tolerance rather than merely “look right” on a map. Getting the projection wrong is not a cosmetic error — a mis-set central meridian or a forgotten false easting relocates a boundary by hundreds of kilometres, and a truncated inverse series leaves a residual that quietly exceeds the millimetre budget a legal plan is graded against. The two conformal families that dominate national cadastral grids are Transverse Mercator, used for narrow north–south zones, and Lambert Conformal Conic, used for wide east–west territories; both are built here from the ellipsoid up, with the scale factor and grid convergence that let a surveyor reduce measured distances and bearings onto the plane.

Figure — the forward projection carries the ellipsoid graticule onto the survey plane; the inverse returns it, and the round-trip residual is the acceptance gate.

Datum and CRS Context: What the Projection Assumes

A projection is defined only relative to a geodetic datum and its ellipsoid. Transverse Mercator on GRS80 and Transverse Mercator on the Airy 1830 ellipsoid are different operations that produce eastings differing by tens of metres for the same input, so the ellipsoid constants — semi-major axis $a$ and flattening $f$ — are part of the projection definition, not an afterthought. The projection also assumes its geographic input is already expressed in the same datum as the projected CRS; if a coordinate is still in a legacy datum, it must first be shifted, and the ellipsoidal-to-Cartesian lift that begins that shift is covered in geodetic conversion math: ellipsoid to Cartesian. The statutory framing of which projection a cadastral jurisdiction mandates, and why conformality (angle preservation) is the governing property for boundary work, is set out in projection math fundamentals for cadastral surveys.

Both families here are conformal: at every point the scale is identical in all directions, so a right angle on the ground remains a right angle on the grid. That is exactly the property a boundary survey needs, because bearings between monuments must be preserved. What is not preserved is scale with distance from the projection’s central line, and that residual distortion is captured by the point scale factor.

The Projection Core: Forward Equations, Scale Factor, and Convergence

From the ellipsoid parameters we derive the second eccentricity and the third flattening used by the modern Krüger series:

e^2 = 2f - f^2, \qquad e'^2 = \frac{e^2}{1-e^2}, \qquad n = \frac{f}{2-f}

The ellipsoidal Transverse Mercator forward mapping is cleanest in the Krüger $n$ -series form, which converts geodetic latitude to the conformal latitude on an auxiliary sphere and then applies a rapidly converging trigonometric series. With $\lambda_0$ the central meridian and $k_0$ the scale factor on it, the conformal quantities are

Q = \operatorname{asinh}(\tan\varphi) - e\,\operatorname{atanh}(e\sin\varphi), \qquad \beta = \arctan(\sinh Q)

\eta_0 = \operatorname{atanh}\!\big(\cos\beta \,\sin(\lambda-\lambda_0)\big), \qquad \xi_0 = \arcsin\!\big(\sin\beta \,\cosh\eta_0\big)

and the northing and easting follow from the series coefficients $h_1 \ldots h_4$ (functions of $n$ ):

\xi = \xi_0 + \sum_{j=1}^{4} h_j \sin(2j\,\xi_0)\cosh(2j\,\eta_0), \qquad \eta = \eta_0 + \sum_{j=1}^{4} h_j \cos(2j\,\xi_0)\sinh(2j\,\eta_0)

E = \mathrm{FE} + k_0 B\,\eta, \qquad N = \mathrm{FN} + k_0\,(B\,\xi - M_0)

where $B = \tfrac{a}{1+n}\left(1 + \tfrac{n^2}{4} + \tfrac{n^4}{64}\right)$ is the rectifying radius and $M_0$ is the meridional value at the latitude of origin $\varphi_0$ .

The point scale factor $k$ and the grid convergence $\gamma$ — the angle between grid north and true north — are what tie the plane back to the ground. To leading terms in $\omega = \lambda - \lambda_0$ , with $t=\tan\varphi$ and $\psi = 1 + e'^2\cos^2\varphi$ ,

k = k_0\left[1 + \tfrac{1}{2}\,\omega^2\cos^2\varphi\,\psi + \cdots\right], \qquad \gamma = \omega\sin\varphi\left[1 + \tfrac{1}{3}\,\omega^2\cos^2\varphi\,(2\psi^2-\psi) + \cdots\right]

A surveyor multiplies a measured horizontal distance by $k$ (and a small height correction) to reduce it to grid, and rotates a measured astronomic azimuth by $\gamma$ to obtain a grid bearing — so a projection library that returns only $(E,N)$ is insufficient for cadastral reduction; $k$ and $\gamma$ must be available at every point.

The Lambert Conformal Conic family replaces the auxiliary sphere with a cone. Its forward mapping uses the isometric parameter $t(\varphi)$ and the parallel radius $m(\varphi)$ :

t = \frac{\tan\!\left(\tfrac{\pi}{4}-\tfrac{\varphi}{2}\right)}{\left(\dfrac{1-e\sin\varphi}{1+e\sin\varphi}\right)^{e/2}}, \qquad m = \frac{\cos\varphi}{\sqrt{1-e^2\sin^2\varphi}}, \qquad n = \frac{\ln m_1 - \ln m_2}{\ln t_1 - \ln t_2}

with $n$ the cone constant fixed by the two standard parallels $\varphi_1, \varphi_2$ . The plane coordinates are then $E = \mathrm{FE} + r\sin\theta$ , $N = \mathrm{FN} + r_0 - r\cos\theta$ , where $r = aF t^{n}$ and $\theta = n(\lambda-\lambda_0)$ .

Step-by-Step Implementation

The module below builds the Transverse Mercator forward mapping, its series inverse, and the scale-factor/convergence pair, then adds Lambert Conformal Conic as a second conformal engine. Everything is float64, type-hinted, and free of external projection libraries so every constant is auditable. A pyproj cross-check appears in the worked example.

Step 1 — Define the ellipsoid and derive the series constants

from __future__ import annotations
import math
import json
import logging
from dataclasses import dataclass

logger = logging.getLogger("projection")


@dataclass(frozen=True)
class Ellipsoid:
    """Reference ellipsoid; defaults are GRS80 / WGS84 (EPSG:7019 / 7030)."""
    a: float = 6378137.0                     # semi-major axis (m)
    inv_f: float = 298.257222101             # inverse flattening (GRS80)

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

    @property
    def e2(self) -> float:
        return self.f * (2.0 - self.f)       # first eccentricity squared

    @property
    def e(self) -> float:
        return math.sqrt(self.e2)

    @property
    def n(self) -> float:
        return self.f / (2.0 - self.f)       # third flattening


def _kruger_forward_coeffs(n: float) -> list[float]:
    # EPSG Guidance Note 7-2, Transverse Mercator (method 9807): h1..h4.
    return [
        n/2.0 - (2.0/3.0)*n**2 + (5.0/16.0)*n**3 + (41.0/180.0)*n**4,
        (13.0/48.0)*n**2 - (3.0/5.0)*n**3 + (557.0/1440.0)*n**4,
        (61.0/240.0)*n**3 - (103.0/140.0)*n**4,
        (49561.0/161280.0)*n**4,
    ]


def _kruger_inverse_coeffs(n: float) -> list[float]:
    # Inverse (grid -> geographic) series coefficients h1'..h4'.
    return [
        n/2.0 - (2.0/3.0)*n**2 + (37.0/96.0)*n**3 - (1.0/360.0)*n**4,
        (1.0/48.0)*n**2 + (1.0/15.0)*n**3 - (437.0/1440.0)*n**4,
        (17.0/480.0)*n**3 - (37.0/840.0)*n**4,
        (4397.0/161280.0)*n**4,
    ]

Step 2 — Parameterise a Transverse Mercator projected CRS

@dataclass(frozen=True)
class TransverseMercator:
    """Ellipsoidal Transverse Mercator (EPSG method 9807).

    lat_origin / lon_origin are the natural origin (degrees); k0 the scale
    factor on the central meridian; false_easting / false_northing in metres.
    """
    ellipsoid: Ellipsoid
    lat_origin_deg: float
    lon_origin_deg: float
    k0: float
    false_easting: float
    false_northing: float

    @property
    def _B(self) -> float:
        n = self.ellipsoid.n
        return (self.ellipsoid.a / (1.0 + n)) * (1.0 + n**2/4.0 + n**4/64.0)

    def _M0(self) -> float:
        # Rectifying value at the latitude of origin (0 when lat_origin == 0).
        phi0 = math.radians(self.lat_origin_deg)
        if phi0 == 0.0:
            return 0.0
        e = self.ellipsoid.e
        q0 = math.asinh(math.tan(phi0)) - e * math.atanh(e * math.sin(phi0))
        beta0 = math.atan(math.sinh(q0))
        xi0 = math.asin(math.sin(beta0))     # eta = 0 on the meridian of origin
        h = _kruger_forward_coeffs(self.ellipsoid.n)
        xi = xi0 + sum(h[j-1] * math.sin(2*j*xi0) for j in range(1, 5))
        return self._B * xi

Step 3 — Implement the forward mapping (φ,λ → E,N)

    def forward(self, lat_deg: float, lon_deg: float) -> tuple[float, float]:
        """Project geodetic (lat, lon) in degrees to (easting, northing) in metres."""
        phi = math.radians(lat_deg)
        lam = math.radians(lon_deg)
        lam0 = math.radians(self.lon_origin_deg)
        e = self.ellipsoid.e
        h = _kruger_forward_coeffs(self.ellipsoid.n)

        q = math.asinh(math.tan(phi)) - e * math.atanh(e * math.sin(phi))
        beta = math.atan(math.sinh(q))
        # atanh input must stay < 1: it does inside a valid TM zone (|dlon| small).
        eta0 = math.atanh(math.cos(beta) * math.sin(lam - lam0))
        xi0 = math.asin(math.sin(beta) * math.cosh(eta0))

        xi = xi0 + sum(h[j-1] * math.sin(2*j*xi0) * math.cosh(2*j*eta0)
                       for j in range(1, 5))
        eta = eta0 + sum(h[j-1] * math.cos(2*j*xi0) * math.sinh(2*j*eta0)
                         for j in range(1, 5))

        easting = self.false_easting + self.k0 * self._B * eta
        northing = self.false_northing + self.k0 * (self._B * xi - self._M0())
        return easting, northing

Step 4 — Implement the series inverse (E,N → φ,λ)

The inverse reverses the same series, then recovers geodetic latitude from the conformal latitude by a short fixed-point iteration. That iteration is why the inverse is not a closed-form negation of the forward: the ellipsoidal correction $Q'' = Q' + e\,\operatorname{atanh}(e\tanh Q'')$ must be iterated until it stops moving.

    def inverse(self, easting: float, northing: float) -> tuple[float, float]:
        """Recover geodetic (lat, lon) in degrees from (easting, northing)."""
        lam0 = math.radians(self.lon_origin_deg)
        e = self.ellipsoid.e
        hp = _kruger_inverse_coeffs(self.ellipsoid.n)

        eta = (easting - self.false_easting) / (self.k0 * self._B)
        xi = ((northing - self.false_northing) / self.k0 + self._M0()) / self._B

        xi0 = xi - sum(hp[j-1] * math.sin(2*j*xi) * math.cosh(2*j*eta)
                       for j in range(1, 5))
        eta0 = eta - sum(hp[j-1] * math.cos(2*j*xi) * math.sinh(2*j*eta)
                         for j in range(1, 5))

        beta = math.asin(math.sin(xi0) / math.cosh(eta0))
        q = math.asinh(math.tan(beta))
        q_iter = q
        for _ in range(10):                  # converges in ~3 passes at survey lat
            q_iter = q + e * math.atanh(e * math.tanh(q_iter))
        lat = math.degrees(math.atan(math.sinh(q_iter)))
        lon = math.degrees(lam0 + math.asin(math.tanh(eta0) / math.cos(beta)))
        return lat, lon

Step 5 — Point scale factor and grid convergence

    def scale_and_convergence(self, lat_deg: float, lon_deg: float) -> tuple[float, float]:
        """Return (point scale factor k, grid convergence gamma in degrees)."""
        phi = math.radians(lat_deg)
        w = math.radians(lon_deg - self.lon_origin_deg)   # omega = lon - CM
        ep2 = self.ellipsoid.e2 / (1.0 - self.ellipsoid.e2)
        cosp, sinp, t = math.cos(phi), math.sin(phi), math.tan(phi)
        psi = 1.0 + ep2 * cosp**2
        t2 = t * t
        k = self.k0 * (
            1.0
            + 0.5 * w**2 * cosp**2 * psi
            + (w**4 / 24.0) * cosp**4 * (4.0*psi**3*(1-6*t2) + psi**2*(1+24*t2) - 4*psi*t2)
        )
        gamma = w * sinp * (
            1.0 + (w**2 / 3.0) * cosp**2 * (2.0*psi**2 - psi) + (w**4 / 15.0) * cosp**4
        )
        return k, math.degrees(gamma)

Step 6 — A second conformal engine: Lambert Conformal Conic (2SP)

@dataclass(frozen=True)
class LambertConformalConic:
    """Lambert Conformal Conic, two standard parallels (EPSG method 9802)."""
    ellipsoid: Ellipsoid
    lat1_deg: float
    lat2_deg: float
    lat_origin_deg: float
    lon_origin_deg: float
    false_easting: float
    false_northing: float

    def _params(self) -> tuple[float, float, float]:
        e = self.ellipsoid.e
        def m(p: float) -> float:
            return math.cos(p) / math.sqrt(1.0 - self.ellipsoid.e2 * math.sin(p)**2)
        def t(p: float) -> float:
            return (math.tan(math.pi/4 - p/2)
                    / (((1 - e*math.sin(p)) / (1 + e*math.sin(p))) ** (e/2)))
        p1, p2 = math.radians(self.lat1_deg), math.radians(self.lat2_deg)
        n = (math.log(m(p1)) - math.log(m(p2))) / (math.log(t(p1)) - math.log(t(p2)))
        F = m(p1) / (n * t(p1)**n)
        r0 = self.ellipsoid.a * F * t(math.radians(self.lat_origin_deg))**n
        return n, F, r0

    def forward(self, lat_deg: float, lon_deg: float) -> tuple[float, float]:
        e = self.ellipsoid.e
        n, F, r0 = self._params()
        p = math.radians(lat_deg)
        t = (math.tan(math.pi/4 - p/2)
             / (((1 - e*math.sin(p)) / (1 + e*math.sin(p))) ** (e/2)))
        r = self.ellipsoid.a * F * t**n
        theta = n * math.radians(lon_deg - self.lon_origin_deg)
        return (self.false_easting + r*math.sin(theta),
                self.false_northing + r0 - r*math.cos(theta))

    def inverse(self, easting: float, northing: float) -> tuple[float, float]:
        e = self.ellipsoid.e
        n, F, r0 = self._params()
        de = easting - self.false_easting
        dn = r0 - (northing - self.false_northing)
        rp = math.copysign(math.hypot(de, dn), n)
        tp = (rp / (self.ellipsoid.a * F)) ** (1.0 / n)
        theta = math.atan2(de, dn)
        lon = math.degrees(theta / n + math.radians(self.lon_origin_deg))
        phi = math.pi/2 - 2*math.atan(tp)
        for _ in range(15):                   # iterate the isometric latitude
            phi = math.pi/2 - 2*math.atan(
                tp * (((1 - e*math.sin(phi)) / (1 + e*math.sin(phi))) ** (e/2)))
        return math.degrees(phi), lon

The parallel between the two classes is deliberate: each exposes forward and inverse, so a pipeline can select the engine from the projected CRS metadata and treat them interchangeably. Choosing which family a jurisdiction actually mandates — and reconciling parcels delivered in a foreign grid — is the wider concern of handling CRS mismatches in cadastral datasets.

Parameter and Return-Value Reference

Name	Type	Units	Valid range	Cadastral significance
`lat_deg`, `lon_deg`	`float`	degrees	−90…90, −180…180	forward input; must be in the projected CRS’s datum
`lon_origin_deg`	`float`	degrees	zone central meridian	a wrong CM shifts every easting by the meridian offset
`lat_origin_deg`	`float`	degrees	natural origin	sets the northing datum via $M_0$
`k0`	`float`	—	0.9996 (UTM), 1.0 (many state planes)	scales the whole grid; misreading it is a systematic error
`false_easting` / `false_northing`	`float`	metres	e.g. 500000 / 0	keeps eastings positive; omitting it yields negative or absurd $E$
`easting` / `northing` (`E`, `N`)	`float`	metres	within zone extent	the projected plane coordinate written to the cadastral record
`k` (return of `scale_and_convergence`)	`float`	—	≈ 0.9996…1.001	multiply measured distance by $k$ to reduce to grid
`gamma`	`float`	degrees	small, signed	rotate true azimuth by $\gamma$ to obtain grid bearing
`inverse()` →	`tuple[float, float]`	degrees	recovered $(\varphi,\lambda)$	must round-trip to the forward input within tolerance

Worked Example: UTM Zone 33N (EPSG:32633)

UTM zone 33N is Transverse Mercator on WGS84 with central meridian $\lambda_0 = 15^{\circ}$ E, latitude of origin $0^{\circ}$ , scale factor $k_0 = 0.9996$ , false easting 500 000 m and false northing 0 m. Projecting a point in central Vienna at $\varphi = 48.208889^{\circ}$ N, $\lambda = 16.372778^{\circ}$ E gives the published grid coordinate, which the code reproduces and which pyproj confirms to well under a millimetre.

wgs84 = Ellipsoid(a=6378137.0, inv_f=298.257223563)   # WGS84 for UTM
utm33 = TransverseMercator(
    ellipsoid=wgs84, lat_origin_deg=0.0, lon_origin_deg=15.0,
    k0=0.9996, false_easting=500000.0, false_northing=0.0,
)

lat, lon = 48.208889, 16.372778
E, N = utm33.forward(lat, lon)
k, gamma = utm33.scale_and_convergence(lat, lon)
print(round(E, 3), round(N, 3))     # 601987.912 5340428.817
print(round(k, 9), round(gamma, 6)) # 0.999727803 1.023603

back = utm33.inverse(E, N)
print(round(back[0], 9), round(back[1], 9))  # 48.208889 16.372778

The scale factor $k = 0.99973$ is above $k_0$ because the point lies $1.37^{\circ}$ east of the central meridian, out where the cylinder cuts the ellipsoid at a slightly larger radius; the grid convergence of $+1.02^{\circ}$ is the angle a surveyor must apply to turn a true bearing into a grid bearing here.

Verification and Round-Trip Residual Analysis

The forward mapping alone is untrustworthy until its inverse recovers the input. The acceptance test is the round-trip residual: project forward, invert, and measure the ground separation between the original and recovered position. Converting the degree difference to metres uses the local radii so the tolerance is expressed where the survey standard lives — in millimetres on the ground.

r_{\text{trip}} = \sqrt{\big(\Delta\varphi \cdot M\big)^2 + \big(\Delta\lambda \cdot \nu\cos\varphi\big)^2}

with $M$ the meridional radius and $\nu$ the prime-vertical radius at $\varphi$ .

def roundtrip_residual_m(proj: TransverseMercator, lat_deg: float, lon_deg: float,
                         tolerance_m: float = 0.001) -> dict[str, object]:
    """Forward-then-inverse a point and emit an audit record of the ground residual."""
    e2 = proj.ellipsoid.e2
    a = proj.ellipsoid.a
    E, N = proj.forward(lat_deg, lon_deg)
    lat_back, lon_back = proj.inverse(E, N)

    phi = math.radians(lat_deg)
    w = 1.0 - e2 * math.sin(phi) ** 2
    nu = a / math.sqrt(w)                       # prime-vertical radius
    M = a * (1.0 - e2) / w ** 1.5               # meridional radius
    d_north = math.radians(lat_back - lat_deg) * M
    d_east = math.radians(lon_back - lon_deg) * nu * math.cos(phi)
    residual = math.hypot(d_north, d_east)

    record = {
        "crs": "EPSG:32633",
        "input_deg": [lat_deg, lon_deg],
        "grid_m": [round(E, 4), round(N, 4)],
        "roundtrip_residual_mm": round(residual * 1000.0, 6),
        "tolerance_mm": tolerance_m * 1000.0,
        "passed": residual <= tolerance_m,
    }
    logger.info("projection_roundtrip %s", json.dumps(record))
    if not record["passed"]:
        raise ValueError(f"Round-trip residual {residual*1e3:.4f} mm exceeds tolerance")
    return record


rec = roundtrip_residual_m(utm33, 48.208889, 16.372778)
assert rec["passed"]           # residual is nanometre-scale for the Kruger series

For the Krüger $n$ -series the round-trip residual is nanometre-scale — far below any survey tolerance — which is exactly the guarantee a certifiable pipeline needs before it commits a projected coordinate. To cross-check against an independent implementation, compare the forward output to pyproj:

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 order
E, N = utm33.forward(48.208889, 16.372778)
assert math.hypot(E - px, N - py) < 1e-3        # agree to sub-millimetre

That millimetre-level agreement with PROJ is the same discipline applied when tuning transformation thresholds for survey-grade work: an independent oracle, a numeric residual, and a hard gate. The step-by-step derivation of the forward Transverse Mercator series on its own — including the meridional-arc form used by the classic Redfearn implementation — is worked in implementing Transverse Mercator forward projection.

Troubleshooting and Gotchas

Eastings come out around 100 km too small or negative

The false easting was dropped. UTM and most cadastral grids add a false_easting of 500 000 m so points west of the central meridian still have positive coordinates. Compute k0 * B * eta and then add the false easting; skipping it leaves the raw distance from the central meridian, which is negative on the west side.

Every coordinate is off by a fixed longitude-sized shift

The central meridian lon_origin_deg is wrong — often the neighbouring zone's. UTM zones are 6° wide with central meridians at 3°, 9°, 15°, …; using zone 32's 9°E for a zone 33 point displaces every easting by the meridian offset times the local radius. Confirm the zone from the EPSG code before projecting.

Northings are consistent but biased by a constant

The latitude of origin or its meridional value M0 is misset. For UTM the latitude of origin is the equator, so M0 = 0; for many state-plane and national grids it is a non-zero parallel and M0 must be subtracted inside the northing term. A constant northing bias with correct eastings is the signature of a wrong M0.

The inverse latitude is close but never quite exact

The fixed-point iteration Q'' = Q' + e·atanh(e·tanh Q'') was run too few times, or not at all. It converges in about three passes at survey latitudes, but truncating it after one leaves a sub-metre latitude error. Iterate until the change falls below your tolerance, or run a fixed ten passes as a deterministic budget.

Distances measured on the ground don't match grid distances

That is the scale factor doing its job, not a bug. A ground distance must be multiplied by the point scale factor k (and reduced for height) before it matches the grid distance between the same two eastings and northings. Near a zone edge k departs from k0 by roughly 400 ppm, which is 0.4 m per kilometre — well outside cadastral tolerance if ignored.

Frequently Asked Questions

Why not just use the spherical Transverse Mercator formulas?

The spherical form ignores the flattening of the Earth and is wrong by hundreds of metres at cadastral latitudes. Every national grid is ellipsoidal, so the projection must use the ellipsoid's eccentricity through the Krüger or Redfearn series. The spherical version is only ever acceptable for small-scale thematic maps, never for boundary work.

How wide can a Transverse Mercator zone be before the series breaks down?

The four-term Krüger series holds sub-millimetre accuracy out to roughly ±4° of longitude from the central meridian, comfortably covering a 6° UTM zone. Beyond that the atanh argument approaches its singularity and the truncated series loses precision; wide territories use Lambert Conformal Conic instead, which is why the conic engine is included here.

Is the point scale factor the same as the zone's k0?

No. k0 is the scale factor fixed on the central meridian (0.9996 for UTM); the point scale factor k varies with distance from that meridian and equals k0 only on it. Distance reduction uses the point value k, computed per coordinate, not the constant k0.

Which ellipsoid should the projection use?

The one bound to the projected CRS, not a default. UTM zones use WGS84; RGF93 / Lambert-93 uses GRS80; the British National Grid uses Airy 1830. Using WGS84 constants for an Airy-based grid introduces a tens-of-metres error that no scale factor will absorb.

Implementing Transverse Mercator forward projection — the forward series derived on its own with the meridional-arc Redfearn form.
Geodetic conversion math: ellipsoid to Cartesian — the ECEF lift that precedes projection when a datum shift is required.
Projection math fundamentals for cadastral surveys — why conformality and grid reduction govern statutory boundary work.
Handling CRS mismatches in cadastral datasets — reconciling parcels delivered in a foreign projected grid.
Algorithmic Math & Geodetic Workflows — the parent reference on the deterministic transformation chain.

Map Projection Forward and Inverse in Python

Datum and CRS Context: What the Projection Assumes #

The Projection Core: Forward Equations, Scale Factor, and Convergence #

Step-by-Step Implementation #

Step 1 — Define the ellipsoid and derive the series constants #

Step 2 — Parameterise a Transverse Mercator projected CRS #

Step 3 — Implement the forward mapping (φ,λ → E,N) #

Step 4 — Implement the series inverse (E,N → φ,λ) #

Step 5 — Point scale factor and grid convergence #

Step 6 — A second conformal engine: Lambert Conformal Conic (2SP) #

Parameter and Return-Value Reference #

Worked Example: UTM Zone 33N (EPSG:32633) #

Verification and Round-Trip Residual Analysis #

Troubleshooting and Gotchas #

Frequently Asked Questions #

Related #