Applying Datum Transformation Grids with pyproj

A datum shift such as NAD27 to NAD83, or NAD27 to WGS84, is not a simple rotation and scale — it corrects for how two reference frames disagree region by region, and that regional distortion is stored in a transformation grid. In pyproj, the accurate path uses an NTv2 grid (horizontal shifts) or a GTX grid (vertical shifts). When that grid is not installed, PROJ does not fail loudly; it silently falls back to a balloted Helmert transform or a null transform and returns coordinates that are off by one to several metres. The fix is to install the grids with projsync or by setting PROJ_NETWORK=ON, build a pyproj.transformer.TransformerGroup, inspect it for the grid-based operation, select that operation explicitly with always_xy=True, and verify the result against a published control point before trusting the pipeline.

# Why a Missing Grid Silently Corrupts Labels

Annotation pipelines routinely ingest legacy vector data — county parcels, historical survey lines, older aerial mosaics — that carries a pre-1983 datum, then reproject it into the coordinate reference system used for training. If the datum transform is wrong by two metres, every polygon vertex lands two metres from its true position. At a ground sample distance of 20 cm/pixel that is a ten-pixel offset — enough to shift a building footprint off its roof and turn correct annotations into systematically mislabelled training data. The failure is insidious because nothing errors out: the code runs, coordinates come back, and only a control-point check reveals the drift.

The root cause is how PROJ ranks operations. For a given datum pair there may be several candidate pipelines, each with a published accuracy. The most accurate uses a grid; coarser ones use a seven-parameter Helmert transform or assume the datums coincide. PROJ prefers the most accurate operation whose grid is present. Remove the grid and it quietly drops to the next candidate. Two machines with different grid sets therefore compute different coordinates for identical input, which also breaks reproducibility across a team.

# Step-by-Step Implementation

Install the pinned toolchain once:

pip install pyproj==3.6.1

pyproj 3.6.1 ships with PROJ 9.3, which distributes all shift grids as GeoTIFF and can stream them from the PROJ content delivery network on demand.

# Step 1 — Enumerate Candidate Operations

TransformerGroup returns every operation PROJ knows for a datum pair, ranked by accuracy, plus a list of operations whose grids are missing. Inspect it before choosing anything:

from pyproj.transformer import TransformerGroup

def describe_operations(source_crs: str, target_crs: str) -> None:
    """Print each candidate operation, its accuracy, and grid availability."""
    tg: TransformerGroup = TransformerGroup(source_crs, target_crs, always_xy=True)
    print(f"available operations: {len(tg.transformers)}")
    for op in tg.transformers:
        acc: float | None = op.operations[0].grids[0].available if op.operations and op.operations[0].grids else None
        print(f"  {op.description!r}  accuracy={op.accuracy} m")
    print(f"unavailable (grid missing): {len(tg.unavailable_operations)}")
    for op in tg.unavailable_operations:
        print(f"  MISSING -> {op.name}")

describe_operations("EPSG:4267", "EPSG:4326")  # NAD27 -> WGS84

If unavailable_operations is non-empty, the most accurate path cannot run yet and PROJ will substitute a coarser one. That is the silent-fallback condition you must fix.

# Step 2 — Enable the Network or Pre-Stage Grids

The fastest fix is to let PROJ fetch grids on demand from its CDN. Enable it in code or by environment variable:

import pyproj

pyproj.network.set_network_enabled(active=True)
print("network enabled:", pyproj.network.is_network_enabled())

Equivalently, export PROJ_NETWORK=ON before the process starts. For reproducible, offline builds — the right choice for CI and containers — pre-stage the exact grids instead of streaming them at runtime with the projsync command-line tool:

# Download only the grids you need (fast, deterministic)
projsync --file us_noaa_conus.tif
projsync --file ca_nrc_ntv2_0.tif

# Or fetch every grid intersecting a bounding box
projsync --bbox -125,24,-66,50 --target-dir /opt/proj_data

Pin the grid set in your image so every machine resolves the same operation. Point PROJ at a fixed directory with PROJ_DATA=/opt/proj_data if you stage grids outside the default location.

# Step 3 — Select the Grid-Based Operation Explicitly

Do not trust the default Transformer.from_crs() result blindly, because it degrades silently. Filter the group for an available, centimetre-accurate operation and construct the transform from it:

from pyproj import Transformer
from pyproj.transformer import TransformerGroup

def best_grid_transformer(source_crs: str, target_crs: str) -> Transformer:
    """Return the most accurate transformer whose grids are all installed."""
    tg: TransformerGroup = TransformerGroup(source_crs, target_crs, always_xy=True)
    if tg.unavailable_operations:
        names = ", ".join(op.name for op in tg.unavailable_operations)
        raise RuntimeError(f"grid-based operation unavailable, install grids: {names}")
    # tg.transformers is sorted best-first; the top entry is grid-based when present
    return tg.transformers[0]

transformer: Transformer = best_grid_transformer("EPSG:4267", "EPSG:4326")
print(transformer.description)

Passing only_best=True to Transformer.from_crs() is the stricter alternative: PROJ then refuses to substitute a lower-ranked operation and errors instead of degrading.

# Step 4 — Verify Against a Published Control Point

A datum transform is only trustworthy once it reproduces an official answer. Pick a coordinate whose value is published by a national survey and assert the residual:

def verify_control_point(transformer: Transformer,
                         src_lon: float, src_lat: float,
                         expected_lon: float, expected_lat: float,
                         tol_deg: float = 5e-6) -> None:
    """Assert the transform reproduces a known control point (~0.5 m at tol_deg=5e-6)."""
    out_lon, out_lat = transformer.transform(src_lon, src_lat)
    dlon: float = abs(out_lon - expected_lon)
    dlat: float = abs(out_lat - expected_lat)
    assert dlon < tol_deg and dlat < tol_deg, (
        f"control point drift too large: dlon={dlon:.7f}, dlat={dlat:.7f} degrees"
    )
    print(f"verified: ({out_lon:.7f}, {out_lat:.7f})")

Because always_xy=True is set on the group, inputs and outputs are always (longitude, latitude) regardless of each CRS’s declared axis order — the single most common source of a swapped-coordinate bug.

# Step 5 — Compare Grid vs Fallback to Quantify the Error

Make the silent degradation visible by transforming the same point both ways and logging the difference. Run this assertion in CI so a missing grid fails the build instead of shipping metre-scale error:

from pyproj import Transformer, Geod

def grid_vs_fallback_metres(source_crs: str, target_crs: str,
                            lon: float, lat: float) -> float:
    """Return the ground distance in metres between grid and fallback results."""
    grid: Transformer = best_grid_transformer(source_crs, target_crs)
    fallback: Transformer = Transformer.from_crs(
        source_crs, target_crs, always_xy=True, allow_ballpark=True
    )
    gx, gy = grid.transform(lon, lat)
    fx, fy = fallback.transform(lon, lat)
    geod: Geod = Geod(ellps="WGS84")
    _, _, dist = geod.inv(gx, gy, fx, fy)
    return dist

delta: float = grid_vs_fallback_metres("EPSG:4267", "EPSG:4326", -96.0, 41.0)
print(f"grid vs fallback: {delta:.2f} m")

# Common Datum Pairs and Their Grids

The error-without-grid column is the ground distance you inherit if PROJ falls back. First reproject into a metric CRS such as EPSG:32614 (UTM 14N) before measuring any of these residuals as distances.

Horizontal grids (NTv2) and vertical grids (GTX) are independent: a pipeline that reprojects points needs the NTv2 grid, while one that converts orthometric to ellipsoidal heights needs the GTX geoid grid. Getting the horizontal shift right does nothing for elevations.

# Common Errors and Fixes

Transform succeeds but coordinates are off by metres Root cause: the NTv2 grid is not installed, so PROJ used a balloted Helmert or null fallback. Fix: build a TransformerGroup, check that unavailable_operations is empty, and install the named grid with projsync or PROJ_NETWORK=ON.

inf or unchanged coordinates returned Root cause: allow_ballpark=False combined with a missing grid means no operation qualifies, so PROJ returns inf. Fix: install the grid, or intentionally set allow_ballpark=True only when metre-scale accuracy is acceptable and documented.

Answers differ between a laptop and CI Root cause: different grid sets on each machine cause PROJ to pick different operations. Fix: pin grids into the container image with projsync, set PROJ_DATA to that directory, and assert the chosen transformer.description in a test.

Latitude and longitude come back swapped Root cause: the target CRS declares latitude-first axis order and always_xy was not set. Fix: pass always_xy=True to TransformerGroup and Transformer.from_crs() so every call is (lon, lat).

Horizontal fix applied but elevations still wrong Root cause: an NTv2 horizontal grid does not correct vertical datums. Fix: add the matching GTX geoid grid (for example us_noaa_g2018u0.tif) and run a separate vertical transform.

# Related

• CRS Roundtrip Testing with pyproj — assert that a reproject-out-and-back cycle stays under tolerance, the natural companion check to a verified datum grid
• Calculating IoU Thresholds for Geospatial Object Detection — why area and overlap must be computed in a metric CRS after any datum correction
• Coordinate Reference Systems in Annotation Pipelines — the broader topic area covering CRS contracts, axis order, and reprojection across a pipeline
• Preserving Metadata Across Dataset Versions — record the exact datum, grid file, and operation description alongside every versioned export

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