SpatiaLite & GeoPackage Automation with Python

Verifying a Reprojection with Known Control Points

Transform a handful of points whose coordinates you hold in both reference systems, measure the geodesic distance between the result and the known value,…

Transform a handful of points whose coordinates you hold in both reference systems, measure the geodesic distance between the result and the known value, and assert that every residual is below a threshold you chose in advance. A residual of a few centimetres means the accurate route ran; one of a metre or two means it did not, and no other check will tell you.

This page belongs to the Datum Transformations & Projection Accuracy guide.

Why This Matters

Every other check in a spatial pipeline verifies structure. The registry rows are present, the reference system is declared, the geometry decodes, the index matches the data — and all of those pass on a container whose coordinates were transformed by the wrong route and sit a metre and a half from where they should.

Control points are the only mechanism that measures the output rather than the configuration. They convert “PROJ says a grid-based route is available” into “the coordinates this pipeline produced are within five centimetres of the published values”, which is a claim about the data rather than about the machine.

Prerequisites

  • At least three control points spread across the working extent, with published coordinates in both systems
  • pyproj 3.4+ with a working proj.db
  • A threshold decided from the capture accuracy, not from what the first run happened to produce
  • The route pinning described in the parent guide

Primary Method

python
# Residual check against published control points
from dataclasses import dataclass
from pyproj import Transformer, Geod


@dataclass(frozen=True)
class ControlPoint:
    name: str
    source: tuple[float, float]     # coordinates in the source system
    expected: tuple[float, float]   # published coordinates in the target system


def residuals(points: list[ControlPoint], src_crs: str, dst_crs: str):
    """Geodesic distance, in metres, between each transformed point and its
    published position. Works for any target system, because the comparison
    is done geodesically rather than in the target's own units."""
    t = Transformer.from_crs(src_crs, dst_crs, always_xy=True)
    geod = Geod(ellps="WGS84")
    to_wgs = Transformer.from_crs(dst_crs, "EPSG:4326", always_xy=True)

    out = []
    for p in points:
        got = t.transform(*p.source)
        got_ll = to_wgs.transform(*got)
        exp_ll = to_wgs.transform(*p.expected)
        _, _, metres = geod.inv(*got_ll, *exp_ll)
        out.append((p.name, metres))
    return out

Measuring geodesically rather than in the target system’s own units is what makes the routine reusable. A residual of “0.043” means nothing until you know whether the units are metres or degrees; a geodesic distance is always metres.

What each layer of checking actually provesA registry check proves the container declares a reference system. A route check proves PROJ can run an accurate transformation on this machine. Only a control-point residual proves that the coordinates in the container are where they should be. The first two are statements about configuration; the third is a statement about the data.registry check — the container declares a systemsays nothing about whether the coordinates match that systemroute check — an accurate transformation is available heresays nothing about which route produced the stored coordinatescontrol-point residual — the coordinates are where they should bethe only check that measures the output rather than the environment
The first two can both pass on a container whose coordinates are a metre and a half out.

Step-by-Step Walkthrough

1. Choose control points that cover the extent

Three points is a working minimum and they should be spread rather than clustered, because the datum shift varies across the area. Points at the corners and the middle of the working extent catch a route that is accurate in one region and not another — which is exactly what a grid with a mismatched area of use produces.

python
CONTROLS = [
    ControlPoint("SW corner",  (170_000.0,  30_000.0), (-5.2011, 50.0930)),
    ControlPoint("centre",     (400_000.0, 300_000.0), (-1.9932, 52.1859)),
    ControlPoint("NE corner",  (620_000.0, 480_000.0), ( 1.0342, 53.7842)),
]

Published values are what make these control points rather than fixtures. Coordinates you generated by running the transformation you are trying to verify prove only that the code is deterministic.

2. Set the threshold from the capture accuracy

The threshold is a decision about acceptable error, and it should be tighter than the accuracy the data claims. A rule that has held up well: set it at a quarter of the capture accuracy, so the transformation contributes a small fraction of the total error budget.

python
CAPTURE_ACCURACY_M = 0.20      # what the survey equipment delivers
THRESHOLD_M = CAPTURE_ACCURACY_M / 4

Setting it from the first run’s observed residual is the anti-pattern. That bakes in whatever route happened to be installed that day, including a wrong one.

3. Run the check and report every point

python
failures = []
for name, metres in residuals(CONTROLS, "EPSG:27700", "EPSG:4326"):
    status = "ok" if metres <= THRESHOLD_M else "FAIL"
    print(f"{status:>4}  {name:<12} residual {metres:7.4f} m")
    if metres > THRESHOLD_M:
        failures.append((name, metres))

if failures:
    raise AssertionError(
        f"{len(failures)} control point(s) beyond {THRESHOLD_M} m: {failures}"
    )

Reporting all of them, rather than raising on the first, is what makes the output diagnostic. One point failing and two passing is a different problem from all three failing by a similar amount.

Reading the residual patternThree diagnostic patterns. All points failing by a similar amount in a similar direction indicates the wrong transformation route ran — usually a missing grid. Points failing progressively worse toward one edge indicates a grid whose area of use does not cover the extent. A single point failing while the others pass usually means that control point's published value is wrong or was transcribed incorrectly.all fail, similar amountand similar directionthe wrong route ranusually a missing gridfix the packagingworse toward one edgecentre fine, corner badgrid area of usedoes not cover the extentfix the route choiceone fails, others passno spatial patternthe control point is wrongtranscription, or a stale valuefix the fixturethe pattern is more informative than the magnitude, which is why every point is reporteda check that raises on the first failure discards the pattern
Three points is the minimum that makes these patterns distinguishable at all.

4. Verify the container, not only the transformer

The check above verifies the transformer object. To verify what the pipeline actually wrote, put a control point into the container as a feature and read it back.

python
# A control point stored as a feature, checked after the reprojection pass
import sqlite3
from shapely import from_wkb

conn = sqlite3.connect("file:survey.gpkg?mode=ro", uri=True)
row = conn.execute(
    "SELECT geom FROM control_points WHERE name = 'centre'"
).fetchone()

pt = from_wkb(gpb_to_wkb(row[0]))
geod = Geod(ellps="WGS84")
_, _, metres = geod.inv(pt.x, pt.y, -1.9932, 52.1859)

assert metres <= THRESHOLD_M, f"stored control point is {metres:.4f} m out"

This is the version worth putting in the publication gate. It survives every step of the pipeline, so it catches a reprojection that ran correctly and a later stage that re-wrote the geometry wrongly.

5. Record the result with the artefact

A residual measured and then discarded proves the run was correct and leaves no evidence. Write the figures into the manifest the pipeline already emits, alongside the PROJ version and the operation used.

python
manifest["crs_verification"] = {
    "source": "EPSG:27700",
    "target": "EPSG:4326",
    "operation": best.description,
    "proj_version": pyproj.proj_version_str,
    "threshold_m": THRESHOLD_M,
    "residuals_m": {name: round(m, 4) for name, m in results},
}

Verification

The check needs its own check: confirm it fails when it should, by running it against a deliberately degraded route.

python
# Force the approximate route and confirm the gate rejects it
from pyproj import Transformer

approximate = Transformer.from_pipeline(
    "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 "
    "+k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy "
    "+step +proj=helmert +x=446.448 +y=-125.157 +z=542.06 "
    "+step +proj=unitconvert +xy_in=rad +xy_out=deg"
)

got = approximate.transform(400_000.0, 300_000.0)
_, _, metres = Geod(ellps="WGS84").inv(*got, -1.9932, 52.1859)
assert metres > THRESHOLD_M, (
    "the gate would not have caught a three-parameter fallback"
)
print(f"degraded route residual {metres:.3f} m — correctly rejected")

That assertion is the one that keeps the gate honest. A threshold set too loosely passes both routes and gives false confidence, and this test is what reveals it.

Choosing a threshold that separates the two routesA scale of residual magnitude. The grid-based route produces residuals of a few centimetres. The seven-parameter fallback produces about a metre. The three-parameter fallback produces several metres. A threshold set between the first and second bands rejects every fallback; one set above the second band passes the approximate routes and gives false confidence.Residual magnitude by routegrid routeunder 0.10 mthreshold0.05 m — set here7-parameter fallbackaround 1 m3-parameter fallbackseveral metresset the threshold from the capture accuracynot from the first run's observed residuala threshold of 2 m passes every route and proves nothing
The gap between the first and third bands is wide, so a threshold anywhere sensible separates them — the failure is choosing one far too loose.

Alternative Approaches or Edge Cases

Where no published control points exist. Use a reference dataset instead: a national mapping agency layer whose positions are authoritative, and compare a sample of features rather than points. The residual is noisier because the features themselves have capture error, but the pattern still distinguishes a route problem from a fixture problem.

Verifying a round trip instead. Transforming out and back and checking closure is cheaper and much weaker — a consistently wrong route round-trips perfectly. It is a useful check for asymmetric grid availability and no substitute for absolute verification.

Extents crossing a grid boundary. Where the working area spans two grids, place control points in each and expect the residuals to differ. A single threshold still applies; what changes is that a failure in one region and not the other is informative rather than confusing.

Troubleshooting

Residuals are around 100 metres

Cause: A datum was ignored entirely — the coordinates were relabelled rather than transformed, or the source system was misidentified. Fix: Check the source reference system against the data’s provenance; a residual of this magnitude is a datum, not a route, as described in Datum Shifts and Why Coordinates Move by Metres.

Residuals are enormous and vary wildly

Cause: Axis order. Longitude and latitude swapped produces residuals that are large and have no spatial pattern. Fix: Confirm always_xy=True is set consistently on every transformer in the path, including the ones inside helper functions.

The check passes in CI and fails in production

Cause: The two environments have different grid availability — the exact failure the check exists to catch, appearing where it should. Fix: Bundle the grids into the runtime image as described in How to Bundle PROJ Grid Files with an Application, and run the check against the runtime image rather than a separate test image.

Frequently Asked Questions

How many control points are enough?

Three to cover an extent and distinguish the failure patterns; five or six if the extent is large or crosses a grid boundary. Beyond that the marginal value falls quickly, because additional points in the same region tell you what the existing ones already did. What matters far more than count is spread — five clustered points are weaker than three spread ones.

Where do published control-point coordinates come from?

National mapping agencies publish coordinates for their geodetic control networks, usually in both the national system and a modern global one, precisely so transformations can be validated. Those are the values to use. Coordinates from a web map, a handheld receiver, or a previous run of your own pipeline are not control points, however precise they look.

Should this run on every build or only on release?

Every build, because it costs milliseconds and the thing it catches — a base-image change altering grid availability — happens without anyone touching the code. It is one of the cheapest checks in a spatial pipeline and the only one that would notice a silent metre-scale shift, which makes the case for running it constantly fairly easy.