Asserting CRS Round-Trip Accuracy in pytest

9 min read

A round-trip assertion transforms a coordinate into a target reference system and back, then bounds how far it moved. It is the cheapest numerical check available on a spatial pipeline and the most commonly written badly, because the obvious version — one point, one threshold, one maximum — passes in almost every situation including the broken ones. This guide sits beneath coordinate reference system testing and shows how to build the version that actually discriminates: a fixture set chosen to span the projection domain, a recorded residual distribution rather than a scalar, and a failure message that names which of the three signature defects occurred.

Root cause: why a single point proves nothing

Projection error is not uniform across a coordinate reference system. It is smallest at the projection’s origin and grows towards the edges of the valid area, because that is what a projection is — a mapping that trades distortion in one place for accuracy in another. A round trip measured at the centre of a zone therefore returns a near-zero residual for almost any transformation, correct or not, and a suite anchored there is measuring the projection’s design rather than the pipeline’s behaviour.

The second reason a naive round trip is weak is structural rather than geographic. The assertion compares T−1(T(p))T^{-1}(T(p)) against pp, which validates the transformation against its own inverse. Two errors that cancel — as they do whenever the forward and inverse both omit the same datum grid — produce a residual of exactly zero. The check is not merely insensitive in that case; it is blind by construction.

r(p)=∥T−1(T(p))−p∥,gate: max⁡p∈Pr(p)≤τr(p) = \lVert T^{-1}(T(p)) - p \rVert, \qquad \text{gate: } \max_{p \in P} r(p) \le \tau

Everything useful about the assertion is in how PP is chosen and what is recorded besides the maximum.

Residual across a projection zone, and where the fixture points sit A curve of round-trip residual plotted across the width of a projection zone. The curve is close to zero in the centre, at the projection origin, and rises steeply towards each edge of the valid area. A horizontal dashed line marks the tolerance. A single fixture point placed at the centre sits well below the tolerance and reports a comfortable pass. Two further points placed near the left and right edges sit above the tolerance line and would have failed the same assertion. A caption states that the placement of the fixture points, not the threshold, decides whether the assertion can discriminate. residual position across the zone west edge origin east edge tolerance the single test point passes comfortably would have failed would have failed The placement of the fixture points, not the threshold, decides whether this assertion can discriminate at all. A centroid drawn from the data lands in the flat region almost every time.

Parameter reference

Element Choice that discriminates Choice that does not
Fixture points Origin, edge, outside, zone boundary, monument A centroid from the data
Reference for comparison At least one externally-published coordinate Only the transformation’s own inverse
Recorded statistic Mean, p95, max, and the worst coordinate The maximum alone
Tolerance unit Metres, after projecting Degrees, as stored
Failure message Rule, point, measured value, threshold “round trip failed”
Transformer construction always_xy=True, errcheck=True Defaults

Step-by-step implementation

The gate below targets pyproj 3.6+, pytest 7+ and a projected CRS in metres.

Step 1 — Declare the fixture set by role

Name each point for what it exercises. The names appear in the failure output, which is what makes a red build interpretable.

import pytest

# Each point is chosen for what it exercises, not sampled from the data.
FIXTURE_POINTS = {
    "origin":        (-2.0, 49.0),      # projection origin — residual should be ~0
    "zone_edge":     (1.76, 55.8),      # maximum scale distortion inside the area
    "outside_area":  (12.0, 55.8),      # must fail loudly, not return a plausible value
    "zone_boundary": (-6.0, 55.0),      # where PROJ may switch operations
}

# Published in both frames — the only external reference in the set.
MONUMENT = {"wgs84": (-1.542324, 53.797416), "target": (429157.19, 434005.51)}

TOLERANCE_M = 0.01

Step 2 — Measure the residual, do not assert yet

Separating measurement from assertion is what allows the distribution to be recorded. A function that returns numbers is also far easier to use from a diagnostic script.

from pyproj import Transformer

FWD = Transformer.from_crs("EPSG:4326", "EPSG:27700", always_xy=True, errcheck=True)
INV = Transformer.from_crs("EPSG:27700", "EPSG:4326", always_xy=True, errcheck=True)

def round_trip_residual_m(lon: float, lat: float) -> float:
    """Displacement in metres after a forward-and-back transform."""
    east, north = FWD.transform(lon, lat)
    lon2, lat2 = INV.transform(east, north)
    # Compare in the projected frame so the residual is already in metres.
    e2, n2 = FWD.transform(lon2, lat2)
    return ((e2 - east) ** 2 + (n2 - north) ** 2) ** 0.5

Step 3 — Assert per point, with the role in the message

Parametrising by name gives one test per role, so a failure identifies which property broke without any further investigation.

@pytest.mark.parametrize("role", ["origin", "zone_edge", "zone_boundary"])
def test_round_trip_within_budget(role):
    lon, lat = FIXTURE_POINTS[role]
    residual = round_trip_residual_m(lon, lat)
    assert residual <= TOLERANCE_M, (
        f"{role}: round-trip residual {residual:.4f} m exceeds {TOLERANCE_M} m "
        f"at ({lon}, {lat})"
    )

def test_outside_valid_area_raises():
    lon, lat = FIXTURE_POINTS["outside_area"]
    with pytest.raises(Exception):
        FWD.transform(lon, lat)     # errcheck=True turns infinity into a raise

Step 4 — Anchor against external truth

The round trip cannot see errors its inverse shares, so one assertion must compare against a value the transformation did not produce.

def test_monument_matches_published_position():
    east, north = FWD.transform(*MONUMENT["wgs84"])
    ex, ny = MONUMENT["target"]
    drift = ((east - ex) ** 2 + (north - ny) ** 2) ** 0.5
    assert drift <= 0.05, (
        f"monument sits {drift:.3f} m from its published position — "
        f"a compensating error the round trip cannot detect"
    )
Each fixture point and the defect only it can catch Five rows pairing a fixture point with the defect class it is uniquely able to detect. The origin point catches structurally wrong transformations such as a swapped axis order or an entirely incorrect target system, because a correct transformation has essentially no residual there. The zone-edge point catches a wrong zone selection or an approximate transformation, because distortion is greatest at the edge. The outside-area point catches a transformation that returns a plausible coordinate instead of failing for input beyond its valid domain. The zone-boundary point catches PROJ selecting different operations for adjacent features, producing a discontinuity mid-dataset. The published monument catches compensating errors shared by the forward and inverse transformations, which no round trip can see. A footer records that dropping any single row removes an entire defect class from coverage. FIXTURE POINT THE DEFECT ONLY IT CATCHES origin structurally wrong transform — swapped axis, wrong target zone edge wrong zone, or an approximate operation outside the valid area a plausible value returned instead of a failure zone boundary two operations applied within one dataset published monument compensating errors — nothing else can see these

Verify the fix

Run the CRS suite on its own, first, so a failure here is not buried among geometry results:

pytest -q tests/test_crs_roundtrip.py -v

Each parametrised case reports its role, so a failing run names the property that broke — zone_edge failing alone points at zone selection, while every case failing together points at the operation or the environment.

Recording the distribution, not the maximum

A gate that emits one number cannot distinguish the two numerical defect signatures, and they have opposite remedies. Record four values per run and the diagnosis becomes immediate: the mean residual, the 95th percentile, the maximum, and the coordinate at which the maximum occurred.

When the mean sits close to the maximum, every point moved by roughly the same amount — a systematic offset, almost always a missing datum grid, fixed by installing the grid package rather than by widening the budget. When the mean is far below the maximum and the worst point is at the edge of the domain, the error scales with position — magnitude-scaled drift, meaning precision loss or a wrong zone, and the only signature where a tolerance conversation is legitimate.

Same maximum, opposite diagnoses Two bar charts of per-point residuals sharing an identical maximum value. In the first, labelled systematic, all five bars are nearly the same height, so the mean sits just below the maximum; the diagnosis given is a missing datum grid and the remedy is to install the grid package. In the second, labelled magnitude-scaled, three bars are near zero and two are tall, so the mean sits far below the maximum; the diagnosis is precision loss or an incorrect zone and the remedy is to partition by zone before considering the threshold. A closing note states that a gate reporting only the maximum presents these two situations identically. Systematic — mean ≈ max max every point moved the same distance → missing datum grid · install it Magnitude-scaled — mean ≪ max max only the edge points moved → precision or zone · partition, then re-measure A gate that reports only the maximum presents these two identically — and their remedies are an environment fix and a data-partitioning decision. Emit the distribution as a run artefact; it costs four numbers.

Where the gate belongs and how often it runs

The whole fixture set transforms in single-digit milliseconds, so cost is not a reason to defer this check. It belongs in the fast pre-merge tier, running before any geometry assertion, because a CRS failure invalidates every geometric result that follows and reporting it first turns a long red build into a short accurate one.

There is a second placement worth adding: the same assertions, run as a start-up guard rather than as tests. A guard that executes before collection and aborts the session on failure produces one clear message instead of a suite full of geometry failures that are all consequences of the same environment problem. The tests remain useful as a record of what is asserted; the guard is what makes a mis-provisioned runner comprehensible.

For the fuller sweep — dozens of points across the domain rather than five — a nightly schedule is the right home. Its purpose is different: not to gate a change, but to watch the residual distribution over time, so a slow degradation caused by a drifting data package becomes visible as a trend rather than as a sudden threshold breach months later.

Failure modes and edge cases

  1. Comparing residuals in degrees. A tolerance of 0.01 is one centimetre in a metric projection and roughly a kilometre in EPSG:4326. Project before measuring, and carry the unit in the constant’s name.
  2. A monument near the diagonal. If the reference point’s easting and northing are numerically similar, an axis-order swap produces almost no drift and passes. Choose a point whose coordinates differ substantially.
  3. Reusing one transformer for both directions. Constructing a single transformer and inverting its arguments can select a different operation than constructing the inverse explicitly. Build both, and assert both.
  4. Silent infinity. Without errcheck=True, an out-of-domain coordinate returns infinity, the residual computes as NaN, and a comparison against a threshold is False — so the assertion passes. Always construct with the check enabled.
  5. Fixture points inside a single zone for multi-zone data. A dataset spanning two UTM zones needs a fixture set per zone; one set measured in the first zone says nothing about the second.
  6. Treating a stable non-zero residual as acceptable. A residual that is small but consistently non-zero at the origin indicates a real transformation difference, not floating-point noise. Noise is random in direction; a constant is not.

Conclusion

A round-trip assertion is worth writing only if it can fail for the right reasons. Choosing fixture points by role rather than by sampling, anchoring at least one comparison against an externally-published coordinate, constructing transformers with explicit axis order and error checking, and recording the residual distribution rather than a maximum turns a check that always passes into one that names its own defect — the fidelity half of coordinate reference system testing.