Datum Shifts and Transformation Grids in PyProj

A coordinate transformation between two datums is not arithmetic — it is a measurement, and PROJ will quietly hand you a cheaper, less accurate measurement rather than fail when the data it needs is missing. This guide is for anyone whose output has to be trusted to the metre or better: surveyors, utility and cadastral teams, and anyone joining historic records to modern ones. It sits under Coordinate Systems with PyProj in Mastering Core Geospatial Python Libraries, and picks up where Fixing PyProj CRS Transformation Errors leaves off — that guide covers transforms that are wrong, this one covers transforms that are merely worse than you think.

Why This Approach / What Goes Wrong

Changing projection while keeping the datum — geographic lon/lat to a UTM zone on the same ellipsoid, say — is a conversion. It is a closed-form formula, exact to floating-point precision, and the only thing that can go wrong is picking the wrong zone (handled in choosing a UTM zone automatically in Python). Changing the datum is a different animal. NAD27, NAD83, ETRS89 and the ITRF realizations each define the shape and the origin of the reference surface differently, so the numeric coordinates of a physical benchmark differ between them even though the brass disk in the ground never moved. There is no formula for that difference — only models fitted to observations.

Two examples set the scale. NAD27 uses the Clarke 1866 ellipsoid with an origin fitted to a single station in Kansas; NAD83 uses GRS80 with a geocentric origin. At a benchmark near 40°N, 105°W the published difference is about 45.7 metres. ETRS89 versus ITRF2014 is a subtler case: ETRS89 is pinned to the stable part of the Eurasian plate, while ITRF tracks the Earth as a whole, so the two frames drift apart by roughly 2.5 cm per year. At Berlin that is currently about 0.63 metres — invisible on a web map, fatal in a cadastral join.

One benchmark expressed in two datums, with the ground offset annotated Two map panels showing the same physical benchmark plotted under two datums. On the left, a benchmark in the Colorado Front Range sits 45.7 metres apart when expressed in NAD27 versus NAD83, because NAD27 uses the Clarke 1866 ellipsoid and a regional origin while NAD83 uses GRS80 and a geocentric origin; a three-parameter Helmert models that shift to only plus or minus 11 metres. On the right, a benchmark in Berlin sits 0.63 metres apart between ETRS89 and ITRF2014, drawn magnified because the offset is roughly one seventieth of the left panel's; that separation grows about 2.5 centimetres a year because ETRS89 is pinned to the Eurasian plate and ITRF is not. The same benchmark, two datums — how far does it move? Real pyproj output on PROJ 9.7 — the ground never moved, the reference frame did. NAD27 → NAD83 · Colorado Front Range NAD27 NAD83 45.7 m apart Clarke 1866 ellipsoid → GRS80, and a new origin fitted to the whole continent. A 3-parameter Helmert models this to only ±11 m. ETRS89 → ITRF2014 · Berlin ETRS89 ITRF2014 0.63 m apart drawn magnified — 1/70th of the left panel's shift Same plate-fixed datum, different global frame: ETRS89 is pinned to Eurasia, ITRF is not. Separation grows ≈2.5 cm a year since 1989.
Both panels show one physical point under two datums; the left offset is metres of historical survey error, the right is centimetres a year of plate motion.

Datum shifts are modelled two ways. A Helmert transformation converts the coordinate to geocentric XYZ, applies three translations (and optionally three rotations plus a scale factor), and converts back. It is a single rigid body move for an entire continent, so it cannot represent the local distortion baked into a triangulation network surveyed with theodolites over eighty years. A transformation grid — NTv2 in Canada, Australia and much of Europe, NADCON and NADCON5 in the United States, all now distributed as GeoTIFF by PROJ — stores an observed shift per grid cell and bilinearly interpolates between them. The grid captures the distortion; the Helmert cannot. For NAD27 to NAD83 the difference in stated accuracy is 0.15 m for the grid against 11 m for the Helmert.

The trap is what PROJ does when the grid file is not on the machine. It does not raise. It walks down a ranked list of candidate operations and runs the best one whose inputs are present — usually a Helmert routed through the WGS 84 pivot, and in the worst case a "ballpark geographic offset" that copies the coordinates through unchanged. TransformerGroup is the object that makes that list visible before you commit to it.

The five NAD27 to NAD83 candidate operations ranked by accuracy and grid requirement A five-row matrix listing every candidate operation PyProj enumerates for NAD27 to NAD83 over the Colorado Front Range. NAD27 to NAD83 seven uses a NADCON5 shift grid at 0.15 metre accuracy and needs us_noaa_nadcon5_nad27_nad83_1986_conus.tif. NAD27 to NAD83 one uses a NADCON grid at 0.15 metres and needs us_noaa_conus.tif. NAD27 to NAD83 four uses an NTv2 grid at 1.5 metres and needs ca_nrc_ntv2_0.tif. NAD27 to NAD83 three uses an NTv1 grid at 2.0 metres and needs ca_nrc_ntv1_can.tif. All four are marked missing on a fresh install. The fifth row, the ballpark geographic offset, needs no grid, is always present, has unknown accuracy and simply copies the coordinates. Every NAD27 → NAD83 candidate PyProj knows about Ranked by stated accuracy — and on a fresh install only the bottom row can actually run. CANDIDATE OPERATION METHOD ACCURACY GRID FILE REQUIRED ON DISK? NAD27 to NAD83 (7) NADCON5 shift grid 0.15 m us_noaa_nadcon5_nad27_ nad83_1986_conus.tif ✗ missing NAD27 to NAD83 (1) NADCON shift grid 0.15 m us_noaa_conus.tif ✗ missing NAD27 to NAD83 (4) NTv2 shift grid · Canada 1.5 m ca_nrc_ntv2_0.tif ✗ missing NAD27 to NAD83 (3) NTv1 shift grid · Canada 2.0 m ca_nrc_ntv1_can.tif ✗ missing Ballpark geographic offset copies the coordinates unknown none ✓ always The only offline-usable row moves the point by exactly zero metres — a 45.7 m error, reported as success.
Accuracy and availability are independent axes: the most accurate operations are the ones that need a file you probably do not have.

Prerequisites

conda install -c conda-forge "pyproj=3.7.*" "geopandas=1.0.*"

Install the whole stack from conda-forge rather than layering a pip wheel over a conda GDAL: two proj.db files on one machine means two different operation tables, and the accuracy you assert in CI is not the accuracy you get in production. Confirm with python -c "import pyproj; pyproj.show_versions()".

Step-by-Step Implementation

The worked example moves a survey benchmark on the Colorado Front Range from NAD27 (EPSG:4267) to NAD83 (EPSG:4269) — the transformation that most often silently degrades when historic records are joined to modern ones.

1. Enumerate the candidate operations before you transform anything.

TransformerGroup asks PROJ the same question Transformer.from_crs() asks, but returns the whole ranked list instead of one answer. Pass an AreaOfInterest so continental candidates that do not cover your data are filtered out.

import pyproj
from pyproj.aoi import AreaOfInterest
from pyproj.transformer import TransformerGroup

pyproj.network.set_network_enabled(False)   # inspect what is on this machine only

FRONT_RANGE = AreaOfInterest(
    west_lon_degree=-105.5, south_lat_degree=39.5,
    east_lon_degree=-104.5, north_lat_degree=40.5,
)

group = TransformerGroup(
    "EPSG:4267", "EPSG:4269", always_xy=True, area_of_interest=FRONT_RANGE
)
print("best operation usable:", group.best_available)
print("usable:", len(group.transformers), "| blocked:", len(group.unavailable_operations))
# best operation usable: False
# usable: 1 | blocked: 4

always_xy=True matters here for the same reason it matters everywhere in PyProj — EPSG:4267 and EPSG:4269 both declare latitude-first axis order, so without it your (lon, lat) tuples arrive transposed. Constructing the group also emits a UserWarning naming the first missing grid; treat that warning as an error in CI rather than letting it scroll past.

2. Read the accuracy ladder and the exact filenames you are missing.

Each entry in unavailable_operations is a CoordinateOperation carrying .accuracy in metres and a .grids list of Grid records with .short_name, .url, .available and .open_license.

for tf in group.transformers:
    print(f"USABLE  {tf.accuracy:>6.2f} m  {tf.description}")

for op in group.unavailable_operations:
    grids = ", ".join(g.short_name for g in op.grids)
    print(f"BLOCKED {op.accuracy:>6.2f} m  {op.name:<20} needs {grids}")

# USABLE   -1.00 m  axis order change (2D) + Ballpark geographic offset from NAD27 to NAD83 + axis order change (2D)
# BLOCKED   0.15 m  NAD27 to NAD83 (7)   needs us_noaa_nadcon5_nad27_nad83_1986_conus.tif
# BLOCKED   0.15 m  NAD27 to NAD83 (1)   needs us_noaa_conus.tif
# BLOCKED   1.50 m  NAD27 to NAD83 (4)   needs ca_nrc_ntv2_0.tif
# BLOCKED   2.00 m  NAD27 to NAD83 (3)   needs ca_nrc_ntv1_can.tif

An accuracy of -1.0 means PROJ declines to state one — the signature of a ballpark operation. Any pipeline whose description contains "Ballpark" applies no datum shift whatsoever.

3. Fetch the grids.

There are three ways, and they suit different deployments. The declarative one is PROJ_NETWORK=ON, which lets PROJ read grid tiles straight from cdn.proj.org over HTTP range requests and cache them in cache.db under the user data directory — no full download, but a network dependency at transform time. The reproducible one is the pyproj sync command-line tool, which pulls whole files into the same directory ahead of time:

# See what covers your extent before downloading anything
python -m pyproj sync --bbox -105.5,39.5,-104.5,40.5 --source-id us_noaa --list-files
# filename | source_id | area_of_use
# ----------------------------------
# us_noaa_conus.tif | us_noaa | USA - Conterminous
# us_noaa_nadcon5_nad27_nad83_1986_conus.tif | us_noaa | USA - Conterminous
# ...

# Then fetch just the file you need into ~/.local/share/proj
python -m pyproj sync --file us_noaa_nadcon5_nad27_nad83_1986_conus --verbose

The third is programmatic: TransformerGroup.download_grids() downloads exactly the grids the blocked operations named, and nothing else.

import pyproj

pyproj.network.set_network_enabled(True)     # or export PROJ_NETWORK=ON
group.download_grids(verbose=True)           # open_license=True by default
print("grids land in:", pyproj.datadir.get_user_data_dir())
# grids land in: /home/you/.local/share/proj

In a container, prefer pyproj sync at build time plus PROJ_NETWORK=OFF at runtime: the image is then self-contained and the transform cannot silently change because a CDN fetch failed. Baking the grids in also removes the per-process cache warm-up that otherwise shows up as latency on the first request.

4. Build a transformer that refuses to guess.

Two keyword arguments narrow what PROJ is allowed to pick, and they do different jobs. allow_ballpark=False rejects only the null-offset operation. only_best=True rejects everything except the top-ranked operation — including the Helmert fallback — so a missing grid becomes visible instead of becoming a 3-metre error.

from pyproj import Transformer

strict = Transformer.from_crs(
    "EPSG:4267", "EPSG:4269",
    always_xy=True,
    area_of_interest=FRONT_RANGE,
    only_best=True,          # top-ranked operation or nothing
    allow_ballpark=False,    # never the no-op offset
)
benchmark_lon, benchmark_lat = -105.0, 40.0
nad83_lon, nad83_lat = strict.transform(benchmark_lon, benchmark_lat)
print(f"{nad83_lon:.8f}, {nad83_lat:.8f}")
# -105.00053503, 39.99998439

5. Prove which pipeline actually ran.

Transformer.description is useless for this: for a CRS-to-CRS transformer it reads "unavailable until proj_trans is called", because PROJ defers operation selection to the first coordinate. get_last_used_operation() returns the operation that was genuinely applied, and it is the only honest answer.

op = strict.get_last_used_operation()
print(op.description)
print("stated accuracy:", op.accuracy, "m")
# axis order change (2D) + NAD27 to NAD83 (7) + axis order change (2D)
# stated accuracy: 0.15 m

Run the same benchmark through a default transformer with the grids absent and the network off, and the contrast is stark — NAD27 to WGS 84 (6) + Inverse of NAD83 to WGS 84 (1), an 11-metre three-parameter Helmert routed through the WGS 84 pivot, reported without a warning and landing 3.34 m from the grid-based answer.

The fallback ladder PROJ walks when a transformation grid is missing Four descending rungs show the order PROJ tries operations for NAD27 to NAD83. Rung one is the NADCON5 grid at 0.15 metre accuracy, needing us_noaa_nadcon5_nad27_nad83_1986_conus.tif. Rung two is a lower-ranked NTv2 or NADCON grid between 0.15 and 2 metres. Rung three is a three-parameter Helmert routed through the WGS 84 pivot at 11 metres, baked into proj.db, always available and never warning. Rung four is the ballpark geographic offset, which copies coordinates unchanged and leaves a 45.7 metre error. Three guards sit alongside: only_best equals True stops at rung one and returns infinity rather than raising when the grid is missing; allow_ballpark equals False blocks only rung four, not the Helmert; and get_last_used_operation is the only reliable way to see which rung ran. Enabling PROJ_NETWORK or running pyproj sync makes rungs one and two reachable. The fallback ladder PROJ walks — silently The first rung whose inputs are present wins, and nothing is written to a log. 1 NADCON5 grid — NAD27 to NAD83 (7) stated accuracy 0.15 m needs us_noaa_nadcon5_nad27_nad83_1986_conus.tif 2 Lower-ranked grids — 0.15 m to 2.0 m tried only when rung 1's file is absent us_noaa_conus.tif · ca_nrc_ntv2_0.tif 3 Helmert via the WGS 84 pivot — 11 m NAD27 to WGS 84 (6) + inverse of NAD83 to WGS 84 (1) 3 parameters in proj.db — always there, never warns 4 Ballpark geographic offset — accuracy unknown copies the coordinates unchanged — 45.7 m off here only_best=True stops here. Grid missing? You get inf, not an exception. allow_ballpark=False does NOT block this rung — it only blocks rung 4. get_last_used_operation() the only reliable way to see which rung actually ran. Set PROJ_NETWORK=ON, or run pyproj sync at build time, and rungs 1–2 become reachable — the ladder stops at the top instead of walking all the way down.
Each rung is a real operation PROJ will run without complaint; the guards on the right are the only things that stop the descent.

6. Convert the shift to metres so the number means something.

Degrees of difference are unreadable. pyproj.Geod turns the before/after pair into a ground distance on the target ellipsoid, which is the number you can put in a test and in a report.

from pyproj import Geod

geod = Geod(ellps="GRS80")
_, _, shift_m = geod.inv(benchmark_lon, benchmark_lat, nad83_lon, nad83_lat)
print(f"datum shift at this benchmark: {shift_m:.2f} m")
# datum shift at this benchmark: 45.72 m

Verification

The test that belongs in CI asserts three separate things: that a result came back at all, that the operation which ran is the grid one, and that the ground shift is physically plausible for the region.

import math
import pyproj
from pyproj import Geod, Transformer
from pyproj.aoi import AreaOfInterest

pyproj.network.set_network_enabled(True)
FRONT_RANGE = AreaOfInterest(-105.5, 39.5, -104.5, 40.5)
benchmark_lon, benchmark_lat = -105.0, 40.0

tf = Transformer.from_crs("EPSG:4267", "EPSG:4269", always_xy=True,
                          area_of_interest=FRONT_RANGE, only_best=True)
nad83_lon, nad83_lat = tf.transform(benchmark_lon, benchmark_lat)

# 1. only_best returns inf instead of raising — catch that first
assert math.isfinite(nad83_lon) and math.isfinite(nad83_lat), \
    "best operation unusable: grid missing and network disabled"

# 2. The pipeline that ran is grid-based, not a Helmert or a null offset
op = tf.get_last_used_operation()
assert "Ballpark" not in op.description, f"no datum shift applied: {op.description}"
assert "NAD27 to NAD83" in op.description, op.description
assert 0 < op.accuracy <= 0.5, f"stated accuracy {op.accuracy} m is too loose"

# 3. The shift is physically plausible for the conterminous US
shift_m = Geod(ellps="GRS80").inv(benchmark_lon, benchmark_lat, nad83_lon, nad83_lat)[2]
assert 10.0 < shift_m < 120.0, f"implausible datum shift: {shift_m:.2f} m"

print(f"pipeline : {op.description}")
print(f"accuracy : {op.accuracy} m   ground shift: {shift_m:.2f} m")
# pipeline : axis order change (2D) + NAD27 to NAD83 (7) + axis order change (2D)
# accuracy : 0.15 m   ground shift: 45.72 m

The op.accuracy assertion is the load-bearing one. A shift-magnitude check alone passes happily on the 11 m Helmert, because 48.96 m is just as "plausible" as 45.72 m — only the stated accuracy of the chosen operation distinguishes them.

Edge Cases & Debugging

Frequently Asked Questions

Do I need grid files if everything I publish is a web map? Almost certainly not. A slippy map tile is a few metres per pixel at city zoom levels, so an 11 m Helmert or even a 1–2 m NAD83-versus-WGS 84 ballpark is invisible. Grids start to matter the moment a coordinate is compared against a survey, a legal boundary, a utility trace or a GNSS fix — anywhere a metre changes an answer. The rule of thumb: if your pipeline ends in a rendered tile, skip the grids; if it ends in a number someone acts on, fetch them and assert the accuracy.

What is the difference between the old .gsb/.las grids and the .tif ones? Only the container. PROJ 7 migrated the whole grid catalogue to GeoTIFF so a single format could carry NTv2, NADCON, geoid and deformation models with consistent metadata and support HTTP range reads — that last property is what makes PROJ_NETWORK=ON practical, since PROJ fetches only the tiles covering your data rather than the whole file. The shift values are identical; a legacy ca_nrc_ntv2_0.gsb and the current ca_nrc_ntv2_0.tif produce the same answer.

How do I stop a PROJ upgrade from silently changing my numbers? Pin pyproj exactly in your lockfile, bake the grids into the image at a known version rather than pulling them from the CDN at runtime, and record op.description and op.accuracy for each CRS pair as a fixture your test suite compares against. PROJ genuinely does add and re-rank operations between releases — NAD27 to NAD83 (7) overtook (1) as the preferred operation when NADCON5 was added — so an unpinned build can change results without a line of your code changing.

Does GeoPandas .to_crs() use the grids too? Yes. GeoPandas delegates entirely to PyProj, so .to_crs() walks the same ladder and honours PROJ_NETWORK and the same data directories. What it does not give you is the inspection surface: there is no way to pass only_best= through .to_crs(). Validate the CRS pair once with TransformerGroup at start-up and fail fast there, then let .to_crs() run on a machine you have already proved is correctly provisioned.

Related