Choosing a Transformation Pipeline When Datums Differ

Between two datums there are usually several published transformations, differing in accuracy by a factor of thirty — and the one your tool picks depends on which grid files happen to be installed. This page is a focused companion to CRS conflict resolution across branches.

Concept & Context

A coordinate reference system conversion between two systems on the same datum is arithmetic and exact. A conversion between datums is a physical model of the relationship between two realisations of the Earth’s shape, and there are always several — a coarse global three-parameter shift, a seven-parameter transformation fitted regionally, and a grid file interpolating measured differences at high density.

PROJ chooses among them by picking the best available operation for the extent, given what is installed. That is a sensible default and a poor guarantee. The same command on two machines with different grid packages produces coordinates up to a metre apart, and neither reports anything unusual. In a versioned repository this shows up as a full-layer diff that nobody made, which is why detecting silent reprojection errors exists as a gate.

The remedy is short: enumerate the candidates, choose one deliberately with its accuracy in view, pin it by name, and record it with the output.

The candidates between two national systems, and what each costs Horizontal bars of stated accuracy for three published transformations between the same pair of systems: a grid-based national operation at a decimetre, a regional seven-parameter transformation at three metres, and a global three-parameter shift at ten. STATED ACCURACY OF EACH CANDIDATE National grid shift 0.1 m only chosen if the grid is installed 7-parameter regional 3 m 3-parameter global 10 m the fallback nobody chose A hundredfold spread, and which one you get depends on which grid package the machine happens to have.

Core Algorithmic Pipeline

  1. Enumerate operations between the source and target CRS for the data’s bounding box.
  2. Filter to those whose area of use covers the whole extent, and whose grids are actually available.
  3. Rank by stated accuracy, then by whether the operation is published by the relevant national authority.
  4. Pin the chosen operation by its authority code in the pipeline parameters.
  5. Record it alongside the output, together with the accuracy and the grid versions used.
Selecting an operation you can defend A chain of four stages: enumerate every operation between the two systems over the data extent, filter to those whose grids are available and whose area of use covers it, rank by stated accuracy, and pin the choice by name. Enumerate over the real extent candidates Filter grids present, area covers usable Rank by stated accuracy best Pin by name, in params Passing the extent matters: without it, an operation valid for one country is offered everywhere.

Working Implementation

"""Enumerate, choose and pin a datum transformation."""
from __future__ import annotations

import json
from dataclasses import dataclass, asdict

import geopandas as gpd
from pyproj import CRS, Transformer
from pyproj.transformer import TransformerGroup


@dataclass
class Candidate:
    name: str
    accuracy_m: float | None
    grids_available: bool
    definition: str

    def usable(self) -> bool:
        return self.grids_available and self.accuracy_m is not None


def candidates(source: str, target: str, bbox: tuple[float, float, float, float]
               ) -> list[Candidate]:
    """Every operation PROJ knows between two systems, over this extent.

    area_of_interest matters: an operation valid for one country is offered for
    that country and ranked poorly elsewhere, and passing no extent hides that.
    """
    from pyproj.aoi import AreaOfInterest

    group = TransformerGroup(
        CRS.from_user_input(source),
        CRS.from_user_input(target),
        area_of_interest=AreaOfInterest(*bbox),
    )

    out = []
    for t in group.transformers:
        out.append(Candidate(
            name=t.description,
            accuracy_m=None if t.accuracy in (None, -1) else float(t.accuracy),
            grids_available=True,
            definition=t.definition,
        ))
    for missing in group.unavailable_operations:
        out.append(Candidate(
            name=missing.name,
            accuracy_m=None if missing.accuracy in (None, -1) else float(missing.accuracy),
            grids_available=False,
            definition="",
        ))
    return out


def choose(source: str, target: str, bbox, required_accuracy_m: float) -> Candidate:
    """Best usable operation, refusing to proceed if none meets the requirement."""
    options = [c for c in candidates(source, target, bbox) if c.usable()]
    if not options:
        raise RuntimeError(
            f"no usable transformation from {source} to {target} over this extent — "
            "the required grid files are not installed"
        )

    best = min(options, key=lambda c: c.accuracy_m)
    if best.accuracy_m > required_accuracy_m:
        raise RuntimeError(
            f"best available transformation is accurate to {best.accuracy_m} m, "
            f"but this data is maintained to {required_accuracy_m} m. Install the "
            "published grid for this area rather than accepting the shortfall."
        )
    return best


def reproject(path: str, out_path: str, target: str,
              required_accuracy_m: float = 0.1) -> dict:
    gdf = gpd.read_file(path)
    source = gdf.crs.to_string()
    bbox = tuple(gdf.total_bounds[[0, 1, 2, 3]])

    picked = choose(source, target, bbox, required_accuracy_m)

    # Pin the operation explicitly: not to_crs(), which re-picks per call.
    transformer = Transformer.from_pipeline(picked.definition)
    reprojected = gdf.copy()
    reprojected["geometry"] = gdf.geometry.apply(
        lambda g: _apply(g, transformer)
    )
    reprojected.set_crs(target, allow_override=True, inplace=True)
    reprojected.to_file(out_path, driver="GPKG")

    import pyproj
    record = {
        "source_crs": source,
        "target_crs": target,
        "operation": asdict(picked),
        "proj_version": pyproj.proj_version_str,
        "proj_data": pyproj.datadir.get_data_dir().rsplit("/", 1)[-1],
        "extent": list(bbox),
    }
    with open(out_path + ".transform.json", "w", encoding="utf-8") as fh:
        json.dump(record, fh, indent=2)
    return record


def _apply(geom, transformer):
    from shapely.ops import transform as shapely_transform
    return shapely_transform(
        lambda x, y, z=None: transformer.transform(x, y), geom
    )

A representative enumeration between two national systems:

operation                                          accuracy  grids
OSGB36 to ETRS89 (OSTN15)                              0.10   yes
OSGB36 to WGS 84 (9)  [7-parameter, national]          3.00   yes
OSGB36 to WGS 84 (1)  [3-parameter, global]           10.00   yes

Thirty times’ difference between the first and second row, and the default picks the first only if the grid is installed. On a runner that never installed it, the same pipeline silently produces the 3 m answer.

Validation & Output Verification

# The chosen operation must be available and recorded
jq -r '.operation.name, .operation.accuracy_m, .proj_data' \
   data/interim/parcels_3035.gpkg.transform.json

# The same transformation must be selected on the runner and on a workstation
python -c "
from transform import choose
c = choose('EPSG:27700', 'EPSG:3035', (-8.0, 49.8, 2.0, 60.9), 0.1)
print(c.name, c.accuracy_m)
"
# Round-trip error must be well under the maintenance tolerance
import geopandas as gpd, numpy as np
from transform import reproject

reproject("data/parcels.gpkg", "/tmp/fwd.gpkg", "EPSG:3035", required_accuracy_m=0.1)
reproject("/tmp/fwd.gpkg", "/tmp/back.gpkg", "EPSG:27700", required_accuracy_m=0.1)

a = gpd.read_file("data/parcels.gpkg").geometry
b = gpd.read_file("/tmp/back.gpkg").geometry
worst = max(pa.centroid.distance(pb.centroid) for pa, pb in zip(a, b))
print(f"worst round-trip displacement: {worst * 1000:.1f} mm")
assert worst < 0.02, "round trip exceeds the survey tolerance"

The round-trip check does not prove the transformation is right — a consistently wrong operation round-trips perfectly. It proves it is stable, which is a necessary first condition; correctness comes from choosing the authority-published operation and recording it.

Failure Modes

  • Coordinates differ between CI and a workstation — symptom: a diff that appears only in one environment. Root cause: different grid packages installed, so PROJ picked different operations. Fix: pin the operation by name and cache the grid package, as in caching GDAL and GeoPandas environments.

  • A layer is accurate in one region and not another — symptom: displacement varying systematically across the extent. Root cause: one operation applied across more than one area of use. Fix: split by area of use and transform each part with the operation valid there.

  • The transformation is not recorded — symptom: a reprojected layer nobody can reproduce. Root cause: to_crs() called without capturing the chosen operation. Fix: record the operation, accuracy and grid version with the output.

  • A required grid silently missing — symptom: a metre-scale shift with no error. Root cause: PROJ fell back to a coarser operation. Fix: fail when the best available operation does not meet the required accuracy, as choose() does.

One operation across two areas of use Two panes over the same national extent. In the first, a single grid-based operation covers only part of the extent, so features outside it fall back to a coarser transformation without warning. In the second, the layer is split by area of use and each part transformed with the operation valid there. ONE OPERATION silent fallback The grid covers the left half; the right half quietly gets a coarser operation. SPLIT BY AREA OF USE Each part transformed with the operation valid there, and both recorded. The first pane produces a layer that is accurate in one region and wrong in another, with no error anywhere.

Back to CRS Conflict Resolution Across Branches