Coordinate Reference System Testing

15 min read

Coordinate reference system testing is the family of assertions that prove a pipeline’s coordinates mean what its metadata claims they mean. It sits beneath Geospatial QA Fundamentals & Architecture as the transformation-and-projection layer of the core pipeline, and it deserves separate treatment from geometry checking for one reason: a CRS defect is systematic. An invalid ring corrupts one feature and a validity predicate finds it; a datum applied without its grid shift file moves every feature by the same one-to-two metres, and no per-feature check will ever notice. This page covers the three properties that must be asserted independently — identity, fidelity and authority — the tooling that makes each testable with pyproj 3.6+ and GeoPandas 0.14+, and the residual signatures that turn a failed assertion into a diagnosis rather than a mystery.

The reason this layer is so often folded into “geometry checks” is that its failures do not look like failures. The geometry stays valid, the schema stays correct, every predicate returns what it always returned. What changes is where the data is, and that is only visible against an external reference — which is precisely what most spatial suites lack.

The Three Properties, and Why They Fail Independently

A CRS assertion set that tests one property and assumes the others is the normal starting point, and it leaves two of the three failure classes entirely uncovered. The properties are genuinely independent: any one can be correct while the other two are wrong.

Identity is the claim that the CRS the data declares is the CRS the data uses. It is a contract property, checked by comparing the declared authority code against what the pipeline manifest expects. A GeoPackage whose srs_id says EPSG:27700 while its coordinates are plainly in degrees fails identity and passes everything else — no numerical tolerance detects it, because nothing about the numbers is wrong in isolation.

Fidelity is the claim that a forward transform followed by its inverse returns the original coordinate within budget. It is a numerical property, and it catches missing grid files, wrong operation selection where several paths exist between two datums, and axis-order swaps that cancel in one direction but not the other.

Authority is the claim that the transformation PROJ actually selected is the one the organisation intends. PROJ frequently offers multiple operations between a datum pair, with accuracies spanning from centimetres to several metres, and it chooses by internal ranking. Nothing about that choice appears in the output coordinates, so an authority failure is invisible until an upgrade changes the ranking and every coordinate moves.

Identity, fidelity and authority — three independent CRS properties Three cards side by side. Identity asserts that the declared coordinate reference system is the one the pipeline contract expects; its failure mode is a driver silently substituting a default definition, and it is detected by comparing the resolved authority code rather than by checking that any coordinate reference system is present. Fidelity asserts that transforming forward and back returns the original coordinate within the round-trip budget; its failure modes are a missing datum grid file and a swapped axis order, and it is detected by measuring the residual across a fixture set that spans the whole valid domain. Authority asserts that the transformation operation PROJ chose is the one intended; its failure mode is a re-ranking after a data package upgrade, and it is detected by naming the operation explicitly and asserting its reported accuracy. A footer records that each property can hold while the other two fail. Identity ASSERTS declared code = contract code FAILS WHEN a driver defaults a definition coordinates are unchanged DETECTED BY comparing to_epsg(), not presence Fidelity ASSERTS round trip returns the point FAILS WHEN a datum grid is missing an axis order is swapped DETECTED BY residuals across the whole domain Authority ASSERTS the chosen operation is intended FAILS WHEN PROJ re-ranks after an upgrade nothing in the output changes shape DETECTED BY pinning the operation + its accuracy Any one of the three can hold while the other two fail. A suite that asserts only fidelity reports identity and authority defects as passing. Assert them separately, and in this order — the cheap one first.

The ordering matters as much as the coverage. Identity is a string comparison that fails in microseconds; fidelity requires transforming a fixture set; authority requires querying the operation PROJ selected. Running them in cost order means a mislabelled dataset fails immediately rather than after a full numerical pass, and the failure names the actual cause rather than its consequences.

Assertion Reference

The table below maps each property to the concrete call that tests it, the tolerance form it uses, and the failure it is designed to catch. All of it targets pyproj 3.6+ against PROJ 9.x.

Property Assertion Tolerance form Catches
Identity — declared gdf.crs.to_epsg() == expected Exact Missing or substituted definition
Identity — axis order Transformer.from_crs(..., always_xy=True) Exact, by construction Latitude/longitude inversion
Fidelity — round trip max⁡p∥T−1(T(p))−p∥≤τ\max_p \lVert T^{-1}(T(p)) - p \rVert \le \tau Absolute, CRS units Missing grid, wrong zone
Fidelity — external Compare against a published monument Absolute, metres Compensating errors in both directions
Authority — operation Assert the selected operation’s name Exact Silent re-ranking after an upgrade
Authority — accuracy Assert the operation’s reported accuracy Upper bound, metres A coarse fallback substituted for a grid
Availability — grids Assert the required grid files are present Exact Environment-dependent drift

The round-trip residual for a coordinate pp under a forward transformation TT is

r(p)=∥T−1(T(p))−p∥r(p) = \lVert T^{-1}(T(p)) - p \rVert

and the gate asserts max⁡p∈Pr(p)≤τ\max_{p \in P} r(p) \le \tau over a fixture set PP chosen to span the CRS domain rather than to sit near its origin. That last qualification is where most round-trip suites are weak — a residual measured at the projection’s centre is close to zero for almost any transformation, correct or not.

Choosing the Fixture Set

A round-trip test is only as informative as the points it round-trips, and the default of “a centroid from the data” is close to useless. Projection error is not uniform: it is smallest at the projection origin and grows towards the edges of the valid area, so a check anchored at the middle of a zone passes while data at the zone boundary is metres out.

A fixture set that reliably catches CRS defects has five members, each exercising something the others cannot:

  1. The projection origin. Residual here is near zero for any correct transform, so a failure means something structurally wrong — a swapped axis, a completely incorrect CRS — rather than a precision issue.
  2. A point near the edge of the valid area. Scale distortion is greatest here, and this is where an incorrect zone or an approximate transformation first becomes visible.
  3. A point just outside the valid area. The transform should fail loudly or return infinity; returning a plausible-looking coordinate is itself a defect worth catching, because it means invalid input will pass silently.
  4. A point on a zone or datum boundary. PROJ may select a different operation for two neighbouring features, producing a discontinuity in the middle of a dataset that no single-point test can see.
  5. A published monument with known coordinates in both systems. The only member that validates against external truth rather than against the transformation’s own inverse.
Compensating errors give a zero round-trip residual A three-state cycle. The source coordinate is transformed forward by an operation that omits its datum grid file, arriving one metre from the true projected position; this intermediate value is what the pipeline actually writes to its output. The inverse transformation, omitting the same grid, moves the value back by an identical vector, returning exactly the source coordinate, so the measured round-trip residual is zero and a round-trip assertion passes. Alongside, a published monument coordinate is compared directly against the intermediate value and reveals the one-metre offset, because that comparison does not use the transformation's own inverse as its reference. A footer generalises the point to any assertion that validates a function against its own inverse. source p geographic, known T(p) — written to the file 1 m from truth T⁻¹(T(p)) = p residual = 0 grid omitted → +1 m grid omitted → −1 m the round-trip assertion passes — the two errors cancelled exactly published monument coordinate compares T(p) against external truth → the 1 m is found Any assertion that validates a function against its own inverse is blind to errors the inverse shares — which is most of the interesting ones. One externally-published coordinate is worth more than a thousand generated round trips.

The fifth member is the one teams skip and the one that matters most. A round-trip test asks whether T−1(T(p))T^{-1}(T(p)) returns pp, and the answer is yes for any pair of mutually consistent transformations — including a pair that are both wrong in the same way. When the pipeline selects an operation that omits a datum grid, the forward and inverse both omit it, the residual is essentially zero, and the suite is green while every coordinate sits a metre from where it belongs.

Reading the Residual

A CI job that reports only the largest round-trip residual throws away the signal that identifies the defect. Bucket the residuals by position instead, and three unmistakable signatures appear — each with a completely different remedy, and only one of them ever legitimately addressed by adjusting a threshold.

A systematic offset displaces every point by the same vector. Almost always a missing datum grid: the transformation fell back to a coarser operation, and the fallback differs from the intended one by a constant. The fix is to install the grid package, never to widen the budget.

Magnitude-scaled drift grows with distance from the projection origin. This is either genuine precision loss or the wrong zone, and it is the only signature where a tolerance discussion is legitimate — though the honest fix is usually to split the data by zone rather than to accept a larger budget.

A mirrored pattern across the diagonal means the axis order was swapped somewhere. Constructing every transformer with explicit axis ordering removes the whole class, and picking reference points whose two coordinates are not numerically close makes the swap detectable when it does occur.

What to record from a round-trip check A record layout for a coordinate reference system round-trip result. It carries the mean residual, the ninety-fifth percentile residual, the maximum residual, the coordinate at which the maximum occurred, the number of fixture points evaluated, and the PROJ library and data package versions in force during the run. Beneath, two interpretations are given: when the mean sits close to the maximum the displacement is systematic and points to a missing grid file, and when the mean is far below the maximum with the worst point at the edge of the domain the drift is magnitude-scaled and points to precision or zone selection. A closing note states that a single maximum value cannot distinguish the two. the record to emit mean_residual_m p95_residual_m max_residual_m max_at_coordinate points_evaluated proj_version / grid_version mean ≈ max every point moved the same distance → systematic: install the grid package mean ≪ max, worst at the edge error grows with distance from origin → precision or zone: split, then re-measure A single maximum cannot distinguish these two, and they have opposite remedies — one is an environment fix, the other a data-partitioning decision. Recording the distribution costs nothing and turns a red build into a diagnosis.

Making the Check Reproducible

A CRS assertion that gives different answers on two machines is worse than none, because the disagreement reads as flakiness and gets retried away. Three things must be pinned, and only the first is obvious.

The PROJ data package version, which supplies the EPSG database and the grid shift files. A transform resolving to a high-accuracy operation on a machine with the full grid package resolves to a coarser fallback on one without it, and the two differ by exactly the amount the grid was correcting for. Install the package explicitly in the runtime image rather than relying on what the base layer happened to ship, and record its version in every log line. The mechanics of that pinning belong to containerized GIS test runtimes.

The operation itself, not just the endpoints. Requesting a transform between two CRSs lets PROJ choose among candidates by its own ranking, which changes between releases. Where accuracy matters, select the operation explicitly and assert on its reported accuracy so an upgrade that would silently substitute a different path fails instead of succeeding differently.

Axis-order handling at construction. Whether a transformer yields easting-northing or northing-easting depends on the authority definition unless stated otherwise, and the two are indistinguishable when a test point sits near the diagonal. Construct with explicit ordering, and choose reference points that are unambiguous.

A start-up guard that asserts the grid files the pipeline depends on are actually present turns the most common cause of environment-dependent drift into a clear message at the beginning of a run, rather than a puzzling numeric difference at the end of it.

Where These Checks Belong in the Pipeline

CRS assertions are cheap and catch systematic defects, which makes them ideal for the fast pre-merge tier described under CI/CD spatial quality gates. Identity checks in particular should run first in the whole suite: they complete in milliseconds, and a CRS mismatch invalidates every geometric result that follows, so failing on it early converts a long red build into a short accurate one.

Fidelity checks belong in the same tier when the fixture set is small — five points transform in negligible time — with a fuller sweep across the CRS domain scheduled nightly. Authority checks are naturally scheduled rather than pre-merge, because what they detect is an environment change rather than a code change, and the environment changes on the image’s cadence rather than on the pull request’s.

One arrangement worth adopting: run the CRS suite as its own job with its own name, so a failure is attributable at a glance. A failure labelled crs-round-trip sends an engineer to the right place immediately; the same failure buried inside a general geometry job sends them to look at polygons.

Common Failure Modes and Gotchas

  1. Testing only the centroid. A round trip at the projection origin is near-exact for almost any transformation, so a suite anchored there passes while data at the zone edge is metres out. Span the domain deliberately.
  2. Checking that a CRS exists rather than which one it is. A driver that writes a default definition produces an artefact that opens cleanly everywhere and attributes the data to the wrong system. Compare the authority code.
  3. Treating a systematic offset as a tolerance problem. Widening the budget until a constant one-metre displacement passes removes the sensitivity that would have caught the next defect, and leaves the original one in place.
  4. Relying on the transformation’s own inverse. Two compensating errors give a zero residual. At least one externally-published coordinate is required for the check to mean anything.
  5. Unpinned grid packages. Two runtimes reporting identical PROJ library versions can disagree by more than a metre if their data packages differ, and nothing in the usual diagnostic output makes that visible.
  6. Assuming the axis order. Test points near the diagonal make a swap undetectable, and the resulting data is mirrored rather than merely offset — a defect that renders as plausible-looking geometry in the wrong hemisphere.

Frequently Asked Questions

How tight should the round-trip budget be?

Derive it rather than inherit it. The floor is the pipeline’s own numerical noise, measured by transforming an unchanged fixture set twice and taking the 99th percentile of the differences. The ceiling is the smallest positional difference a consumer can detect — the map scale the data is rendered at, or the survey standard it must satisfy. Any value between those is defensible; anything below the floor produces false failures and anything above the ceiling lets real defects pass. If the two bounds have crossed, the pipeline has a precision problem that no threshold can express.

Do we need to test transformations we do not perform ourselves?

Yes, when a driver performs them on your behalf. Writing a projected layer to GeoJSON triggers a reprojection to geographic coordinates whether or not your code asked for one, and that transformation is subject to every failure mode described here. The rule of thumb: if a format has an implied CRS, the write is a transformation and it belongs in the fidelity suite.

What is a reasonable size for the fixture set?

Five to a dozen points per CRS pair, chosen for what each exercises rather than sampled. Adding more points of the same kind does not increase coverage — a hundred coordinates clustered near the projection origin all report the same near-zero residual. What increases coverage is adding a kind of point: another zone boundary, another edge of the valid area, another published monument.

How do we handle datasets that legitimately span multiple zones?

Test per zone rather than per dataset, and make the zone boundary an explicit fixture. A single budget across several zones is either too loose in the middle of each or too tight at the edges, and it hides the discontinuity that appears when PROJ selects different operations for adjacent features. Partitioning the assertion the same way the data is partitioned keeps each measurement interpretable.

Should CRS assertions run against the data or against the transformation?

Both, and they answer different questions. Asserting against the transformation — round-tripping known coordinates — verifies the environment and the operation selection, and it works without any data at all. Asserting against the data verifies that this particular dataset carries what it claims. The first belongs in a start-up guard that runs before the suite; the second belongs with the contract checks that run over every batch.

Is a vertical datum worth testing separately?

If the pipeline carries elevation, yes, and it is routinely forgotten. A horizontal transformation that is correct says nothing about whether heights were converted between an ellipsoidal and a geoid-based reference, and the two can differ by tens of metres. Treat the vertical component as its own identity-and-fidelity pair, with its own fixture points and its own budget.

Conclusion

Coordinate reference system testing is a separate assertion family because its failures are separate: systematic rather than per-feature, invisible to validity predicates, and detectable only against an external reference or a deliberately-spanning fixture set. Asserting identity, fidelity and authority independently — in that order, with the residual distribution recorded rather than a single maximum, against a pinned PROJ data package — turns the most expensive class of spatial defect into one that fails fast and names its own cause. It is the transformation-and-projection layer that the rest of geospatial QA fundamentals and architecture assumes is already in place.