CRS Roundtrip Testing with pyproj

A CRS roundtrip test reprojects a set of known control points from the source coordinate reference system to the target and back to the source, then asserts that the maximum coordinate drift stays below a tolerance — under 1 cm for grid-based datum transforms, sub-millimetre for a plain projection change. That single assertion catches the three failures that silently corrupt geospatial training labels: axis-order bugs, missing datum-shift grids, and lossy datum operations. Build a forward and a back pyproj Transformer with always_xy=True, run your control points through both, measure the residual, and fail the build when it exceeds the tolerance for the transform class you are exercising. This guide shows the runnable pattern and the tolerance table to calibrate it against.

# Why Roundtrip Drift Predicts Corrupted Labels

Reprojection sits on the critical path of almost every annotation pipeline: imagery arrives in one coordinate reference system, annotators draw in another, and the training exporter writes a third. Every hop is a Transformer call, and a wrong one does not raise — it returns plausible numbers that are quietly off by metres. A polygon shifted two metres north still looks like a valid building footprint in a viewer, but it now disagrees with the pixels it is supposed to label, and the model learns the offset as signal.

The roundtrip is powerful because it is self-checking. A correct forward transform composed with a correct inverse is the identity, so the returned coordinate must equal the original to within numerical noise. Any measurable residual means one of the two directions is wrong, and the magnitude tells you which failure class you hit. A residual of tens of degrees is an axis swap. A residual of one to ten metres is a missing grid falling back to an approximate shift. A residual of a few millimetres on a supposedly exact transform points at a lossy intermediate step. You do not need ground truth to run the check — the source coordinate is its own reference — which is what makes it cheap enough to run on every commit.

# Step-by-Step Implementation

Install the pinned toolchain once:

pip install pyproj==3.6.1 numpy==1.26.4 pytest==8.2.2

# Step 1 — Define Control Points and CRS Pairs

Control points are ordinary coordinates in the source CRS, chosen to span the extent of your annotation tiles. Use the corners and the centre rather than a single point, because datum residuals vary across an area. Coordinates here are (lon, lat) in [EPSG:4326](https://www.geospatialannotation.com/geospatial-annotation-fundamentals-architecture/coordinate-reference-systems-in-annotation-pipelines/):

from dataclasses import dataclass

@dataclass(frozen=True)
class RoundtripCase:
    source_crs: str
    target_crs: str
    tolerance_m: float          # tolerance expressed in metres on the ground
    control_points: list[tuple[float, float]]   # (x, y) in source CRS order

CENTRAL_EUROPE = RoundtripCase(
    source_crs="EPSG:4326",
    target_crs="EPSG:32633",    # UTM zone 33N
    tolerance_m=1e-3,           # equal-datum projection change: sub-mm expected
    control_points=[
        (13.30, 52.45),         # SW corner
        (13.60, 52.45),         # SE corner
        (13.30, 52.60),         # NW corner
        (13.60, 52.60),         # NE corner
        (13.45, 52.52),         # centre
    ],
)

# Step 2 — Build Forward and Back Transformers

Both transformers must carry always_xy=True so axis order is explicit and identical in each direction. The only difference between them is which CRS is source and which is target:

from pyproj import Transformer

def build_transformers(
    source_crs: str,
    target_crs: str,
) -> tuple[Transformer, Transformer]:
    """Return (forward, back) transformers that share xy axis order."""
    forward = Transformer.from_crs(source_crs, target_crs, always_xy=True)
    back = Transformer.from_crs(target_crs, source_crs, always_xy=True)
    return forward, back

# Step 3 — Run the Roundtrip and Measure Drift

Project each point forward then back, and measure the residual. Because the source CRS here is geographic (degrees), convert the degree residual to an approximate ground distance in metres so the comparison against a metric tolerance is dimensionally correct. A local scale factor is sufficient for a tolerance check:

import math
import numpy as np

def degrees_to_metres(
    dlon: float,
    dlat: float,
    at_lat_deg: float,
) -> float:
    """Approximate ground distance for a small lon/lat delta near a latitude."""
    m_per_deg_lat = 111_132.0
    m_per_deg_lon = 111_320.0 * math.cos(math.radians(at_lat_deg))
    return math.hypot(dlon * m_per_deg_lon, dlat * m_per_deg_lat)

def max_roundtrip_drift(case: RoundtripCase) -> float:
    """Return the largest source→target→source drift, in metres, over all points."""
    forward, back = build_transformers(case.source_crs, case.target_crs)
    drifts: list[float] = []
    for x, y in case.control_points:
        tx, ty = forward.transform(x, y)
        rx, ry = back.transform(tx, ty)
        if case.source_crs.upper().endswith("4326"):
            drift = degrees_to_metres(rx - x, ry - y, at_lat_deg=y)
        else:
            drift = math.hypot(rx - x, ry - y)   # already metric
        drifts.append(drift)
    return float(np.max(drifts))

# Step 4 — Assert the Tolerance

The assertion is one line, but the message matters — when the gate fails months later, the reported drift and tolerance are what let an engineer classify the failure immediately:

def assert_roundtrip(case: RoundtripCase) -> float:
    """Raise AssertionError if max drift exceeds the case tolerance."""
    drift = max_roundtrip_drift(case)
    assert drift <= case.tolerance_m, (
        f"{case.source_crs}->{case.target_crs} roundtrip drift "
        f"{drift:.4g} m exceeds tolerance {case.tolerance_m:.4g} m — "
        f"check axis order and datum-shift grids"
    )
    return drift

# Step 5 — Wrap It as a pytest Gate

Parametrize over every CRS pair your pipeline touches so a single test file guards the whole reprojection surface. Add the grid-based cases you rely on — the datum transformation grids applied with pyproj determine whether those cases hit their tight tolerance or fall back and fail:

import pytest

CASES = [
    CENTRAL_EUROPE,
    RoundtripCase(
        source_crs="EPSG:4267",     # NAD27
        target_crs="EPSG:4326",     # WGS84 (needs an NTv2 grid)
        tolerance_m=0.01,           # grid-based datum transform: sub-cm
        control_points=[(-122.30, 47.60), (-122.30, 47.62), (-122.28, 47.61)],
    ),
]

@pytest.mark.parametrize("case", CASES, ids=lambda c: f"{c.source_crs}->{c.target_crs}")
def test_crs_roundtrip_within_tolerance(case: RoundtripCase) -> None:
    drift = assert_roundtrip(case)
    assert drift >= 0.0

Run it with pytest -q. A green result certifies that every registered transform is invertible to within its tolerance; a red result names the offending CRS pair and the measured drift.

# Tolerance by Transform Type

Drift is not a single number — the achievable residual depends on what the transform actually does. Assert against the class of transform you are exercising, not a universal constant:

The last row is the whole point of the exercise. When the correct grid is present the transform is invertible to a centimetre; when it is missing, pyproj substitutes an approximate operation and the roundtrip no longer cancels. A tolerance of 1 cm on a grid-based case turns that silent substitution into a red build.

# Common Errors and Fixes

Roundtrip drift is tens of degrees, not centimetres Root cause: axis-order mismatch — one transformer was built with always_xy=True and the other without, so latitude and longitude are swapped on the return leg. Fix: set always_xy=True identically on both the forward and back Transformer, and never mix a transformer built one way with one built the other.

Drift is a stable 1–10 m across all points Root cause: the datum-shift grid for the transform is not installed, so pyproj fell back to a null or Helmert approximation and the two directions do not cancel. Fix: install the grids with projsync --all or the pyproj-data package, then confirm with Transformer.from_crs(...).description that the operation names the expected grid rather than a ballpark transform.

The test passes at 1e-6 but should have caught a real shift Root cause: the tolerance was compared in the wrong units — a drift measured in degrees was tested against a tolerance meant for metres, so a real 2 m error (about 2e-5 degrees) slipped under a 1e-6 bar that looked strict. Fix: convert the residual to ground metres before comparison, as degrees_to_metres does, and always state the tolerance in the same units you measure the drift in.

# Related

• Coordinate Reference Systems in Annotation Pipelines — the parent guide covering CRS contracts, datum management, and reprojection patterns that a roundtrip test defends
• Applying Datum Transformation Grids with pyproj — install and verify the NTv2/GTX grids whose absence a roundtrip test is designed to expose
• CI/CD gates for annotation datasets — wire the pytest roundtrip check into a pipeline that blocks a bad reprojection from reaching training
• Calculating IoU Thresholds for Geospatial Object Detection — a companion reprojection-sensitive calculation that assumes the transforms a roundtrip test validates

This guide covers one focused check within Coordinate Reference Systems in Annotation Pipelines, which is itself part of Geospatial Annotation Fundamentals & Architecture.