How to Reproject Between Datums Correctly (and Why the Grid File Matters) #

Problem statement #

The same reprojection, run twice, gives different answers:

from pyproj import Transformer

t = Transformer.from_crs(4326, 27700, always_xy=True)
print(t.transform(-2.2426, 53.4808))

# March, on a laptop with PROJ 9.2 and no network
(383618.507, 398050.393)

# August, in a container with PROJ 9.4 and PROJ_NETWORK=ON
(383620.114, 398048.881)

1.6 metres apart. Same input, same EPSG codes, same library — different transformation pipeline, because one run had the OSTN15 grid file available and the other did not.

For a map at 1:50,000 that is invisible. For a cadastral boundary, a utility asset, or anything with legal consequence, it is the difference between correct and wrong. And nothing warns you: both runs return a plausible float.

Getting datum transformations right means knowing that there is usually more than one way to do it, and choosing deliberately.

Quick answer #

from pyproj import Transformer, TransformerGroup

# see every available pipeline, ranked by accuracy
group = TransformerGroup(4326, 27700, always_xy=True)
for t in group.transformers[:3]:
    print(f"{t.description}\n   accuracy {t.accuracy} m")

Inverse of OSGB36 to WGS 84 (9) + British National Grid
   accuracy 1.0
Inverse of OSGB36 to WGS 84 (5) + British National Grid
   accuracy 3.0
Ballpark geographic offset from WGS 84 to OSGB36 + British National Grid
   accuracy -1.0

import os
os.environ["PROJ_NETWORK"] = "OFF"      # before importing pyproj

conda install -c conda-forge proj-data  # the full grid set, ~600 MB

Always log Transformer.description with any coordinates you publish. It names the pipeline that produced them.

Step-by-step solution #

1. Find out what your transformation actually does #

from pyproj import Transformer

t = Transformer.from_crs("EPSG:4326", "EPSG:27700", always_xy=True)
print(t.description)
print(f"accuracy   {t.accuracy} m")
print(f"operations {t.operations if hasattr(t, 'operations') else '—'}")

Inverse of OSGB36 to WGS 84 (9) + British National Grid
accuracy   1.0

OSGB36 to WGS 84 (9) is the OSTN15 grid-based transformation, accurate to about a metre. Compare with what you get without the grid:

Ballpark geographic offset from WGS 84 to OSGB36 + British National Grid
accuracy   -1.0

"Ballpark" is PROJ telling you it gave up. It means no proper transformation was available, so it treated the two datums as if they coincided — an error of up to about 120 m in Britain. An accuracy of -1.0 means unknown, which in practice means unquantified and large.

Never ship coordinates produced by a ballpark transformation without saying so.

2. Understand what a datum shift is #

A geodetic datum is an ellipsoid plus its position and orientation relative to the earth. Different datums use different ellipsoids, fitted to different regions in different eras:

So converting between datums moves the point — it is a correction, not a rounding. A transformation that ignores it is wrong by the offset.

from pyproj import Transformer, Geod

geod = Geod(ellps="WGS84")
lon, lat = -2.2426, 53.4808

proper = Transformer.from_crs(4326, 27700, always_xy=True)
x1, y1 = proper.transform(lon, lat)
back = Transformer.from_crs(27700, 4326, always_xy=True).transform(x1, y1)
_, _, err = geod.inv(lon, lat, back[0], back[1])
print(f"round trip error: {err * 100:.2f} cm")

round trip error: 0.03 cm

A clean round trip is evidence the same pipeline ran both ways. A round trip that loses metres means the two directions chose different transformations.

3. Make sure the grid files are available #

Grid-based transformations model the datum difference as a field that varies with position, stored as a raster. They are far more accurate than a global 7-parameter formula because the difference genuinely varies — Britain's crust does not shift uniformly.

import pyproj

print(f"data dir  {pyproj.datadir.get_data_dir()}")
print(f"network   {pyproj.network.is_network_enabled()}")

from pathlib import Path
grids = sorted(Path(pyproj.datadir.get_data_dir()).glob("*.tif"))
print(f"grids     {len(grids)} files")
for g in grids[:6]:
    print(f"          {g.name}")

data dir  /opt/conda/share/proj
network   False
grids     412 files
          au_ga_AGQG_20201120.tif
          ca_nrc_ntv2_0.tif
          uk_os_OSGB15_GB.tif
          us_noaa_conus.tif
          …

uk_os_OSGB15_GB.tif is OSTN15. Without it, WGS 84 to OSGB36 falls back to ballpark.

Install the full set:

conda install -c conda-forge proj-data     # ~600 MB, everything
pip install pyproj                          # bundles only a minimal set

Or fetch what you need:

projsync --list-files
projsync --area-of-use "United Kingdom"
projsync --file uk_os_OSGB15_GB.tif

4. Turn network fetching off for reproducibility #

import os
os.environ["PROJ_NETWORK"] = "OFF"      # must be set before pyproj is imported
import pyproj
print(pyproj.network.is_network_enabled())    # False

With PROJ_NETWORK=ON, PROJ downloads grids from a CDN when a transformation needs one, caches them, and uses them. That is convenient and it means the same code gives different answers depending on network availability and cache state.

import pyproj
from pyproj import Transformer

for enabled in (True, False):
    pyproj.network.set_network_enabled(enabled)
    t = Transformer.from_crs(4326, 27700, always_xy=True)
    x, y = t.transform(-2.2426, 53.4808)
    print(f"network={str(enabled):<5} {x:>11.3f} {y:>11.3f}   {t.description[:44]}")

network=True   383618.507  398050.393   Inverse of OSGB36 to WGS 84 (9) + British Na…
network=False  383620.114  398048.881   Ballpark geographic offset from WGS 84 to OS…

On a machine with no local grid, network access decides whether you get a metre-accurate answer or a ballpark one. For anything reproducible: install the grids and turn the network off, so a missing grid is a visible degradation rather than a silent dependency on connectivity.

5. Choose a transformation explicitly when it matters #

from_crs picks the highest-accuracy available pipeline. When the choice must be recorded and stable, name it:

from pyproj import TransformerGroup, Transformer

group = TransformerGroup("EPSG:4326", "EPSG:27700", always_xy=True)
print(f"{len(group.transformers)} available, "
      f"{len(group.unavailable_operations)} unavailable")
for i, t in enumerate(group.transformers):
    print(f"  [{i}] accuracy {str(t.accuracy):>6} m  {t.description[:56]}")
for op in group.unavailable_operations[:3]:
    print(f"  ✗ {op.name[:50]} — needs {[g.name for g in op.grids]}")

3 available, 1 unavailable
  [0] accuracy    1.0 m  Inverse of OSGB36 to WGS 84 (9) + British National…
  [1] accuracy    3.0 m  Inverse of OSGB36 to WGS 84 (5) + British National…
  [2] accuracy   -1.0 m  Ballpark geographic offset from WGS 84 to OSGB36 …
  ✗ OSGB36 to ETRS89 (2) — needs ['uk_os_OSGB15_GB.tif']

unavailable_operations lists what a missing grid is costing you, and names the file to install.

Pin the pipeline by its identifier:

t = Transformer.from_pipeline(
    "+proj=pipeline "
    "+step +proj=axisswap +order=2,1 "
    "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
    "+step +inv +proj=hgridshift +grids=uk_os_OSGB15_GB.tif "
    "+step +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 "
    "+x_0=400000 +y_0=-100000 +ellps=airy")

That is verbose and completely unambiguous — it will do the same thing on every machine that has the grid, and fail loudly on one that does not.

The gentler option is an accuracy floor:

t = Transformer.from_crs(4326, 27700, always_xy=True, accuracy=5.0)

which raises rather than falling back to ballpark when nothing meets the requirement.

6. Record the pipeline alongside the results #

import json
from datetime import datetime, timezone
from pathlib import Path
import pyproj
from pyproj import Transformer

def transform_with_provenance(gdf, target_crs, out_json=None):
    src = gdf.crs
    t = Transformer.from_crs(src, target_crs, always_xy=True)
    out = gdf.to_crs(target_crs)

    record = {
        "source_crs": src.to_string(),
        "target_crs": pyproj.CRS.from_user_input(target_crs).to_string(),
        "pipeline": t.description,
        "accuracy_m": t.accuracy,
        "proj_version": pyproj.proj_version_str,
        "pyproj_version": pyproj.__version__,
        "proj_data_dir": str(pyproj.datadir.get_data_dir()),
        "proj_network": pyproj.network.is_network_enabled(),
        "when": datetime.now(timezone.utc).isoformat(),
    }
    if "Ballpark" in t.description:
        print(f"  ⚠ ballpark transformation — error may exceed 100 m")
    if out_json:
        Path(out_json).write_text(json.dumps(record, indent=2))
    return out, record

projected, record = transform_with_provenance(
    gdf, 27700, out_json="out/transformation.json")
print(record["pipeline"], "—", record["accuracy_m"], "m")

Six months later, when the numbers differ from a colleague's, the JSON says which pipeline, which PROJ version and whether the network was on. Without it, the question is unanswerable. This is the coordinate half of recording run metadata and data lineage.

Code examples #

Example 1: a pre-flight check for transformation accuracy #

import pyproj
from pyproj import CRS, TransformerGroup

def check_transformation(src, dst, *, required_accuracy_m=5.0, verbose=True):
    """Verify a datum transformation is available and accurate enough."""
    src_crs = CRS.from_user_input(src)
    dst_crs = CRS.from_user_input(dst)
    group = TransformerGroup(src_crs, dst_crs, always_xy=True)

    if verbose:
        print(f"{src_crs.name} → {dst_crs.name}")
        print(f"  datum: {src_crs.datum.name if src_crs.datum else '—'}")
        print(f"      →  {dst_crs.datum.name if dst_crs.datum else '—'}")

    usable = [t for t in group.transformers
              if t.accuracy is not None and 0 <= t.accuracy <= required_accuracy_m]
    ballpark = [t for t in group.transformers if "Ballpark" in t.description]

    if verbose:
        for t in group.transformers[:4]:
            acc = "unknown" if t.accuracy in (None, -1.0) else f"{t.accuracy} m"
            mark = "✓" if t in usable else "✗"
            print(f"  {mark} {acc:>9}  {t.description[:56]}")
        for op in group.unavailable_operations[:3]:
            names = [g.name for g in op.grids]
            print(f"  ✗ unavailable: {op.name[:44]}")
            print(f"      needs {names} — projsync --file {names[0] if names else '?'}")

    if not usable:
        raise RuntimeError(
            f"no transformation from {src_crs.name} to {dst_crs.name} meets "
            f"{required_accuracy_m} m. Best available: "
            f"{group.transformers[0].description}. "
            f"Install proj-data or run projsync for the missing grids.")

    if ballpark and group.transformers[0] in ballpark:
        raise RuntimeError(
            f"only a ballpark transformation is available — error may exceed 100 m")

    best = usable[0]
    if verbose:
        print(f"  → using: {best.description} ({best.accuracy} m)")
    return best

check_transformation(4326, 27700, required_accuracy_m=5.0)
check_transformation(4267, 4326, required_accuracy_m=5.0)     # NAD27 → WGS 84

WGS 84 → OSGB36 / British National Grid
  datum: World Geodetic System 1984 ensemble
      →  Ordnance Survey of Great Britain 1936
  ✓       1.0 m  Inverse of OSGB36 to WGS 84 (9) + British National Grid
  ✓       3.0 m  Inverse of OSGB36 to WGS 84 (5) + British National Grid
  ✗   unknown  Ballpark geographic offset from WGS 84 to OSGB36 + Brit…
  → using: Inverse of OSGB36 to WGS 84 (9) (1.0 m)

Running this at the start of a pipeline turns a silent 100 m error into a startup failure that names the grid file to install. The projsync hint is the part that makes it actionable rather than merely alarming.

The accuracy of -1.0 for ballpark is a PROJ convention meaning unknown — which is why the filter checks 0 <= accuracy rather than just an upper bound.

Example 2: quantifying what a wrong transformation costs #

import numpy as np
import geopandas as gpd
from pyproj import Transformer, Geod, CRS

def compare_transformations(points_4326, target_crs, *, ellps="WGS84"):
    """How far apart do the available pipelines put the same points?"""
    geod = Geod(ellps=ellps)
    from pyproj import TransformerGroup
    group = TransformerGroup(4326, target_crs, always_xy=True)

    lons = points_4326.geometry.x.to_numpy()
    lats = points_4326.geometry.y.to_numpy()

    results = {}
    for t in group.transformers:
        xs, ys = t.transform(lons, lats)
        results[t.description[:46]] = (np.asarray(xs), np.asarray(ys), t.accuracy)

    names = list(results)
    best_x, best_y, _ = results[names[0]]
    back = Transformer.from_crs(target_crs, 4326, always_xy=True)

    print(f"{len(names)} pipelines, {len(lons)} points\n")
    for name in names:
        xs, ys, acc = results[name]
        dx, dy = xs - best_x, ys - best_y
        offset = np.hypot(dx, dy)
        blon, blat = back.transform(xs, ys)
        _, _, round_trip = geod.inv(lons, lats, blon, blat)
        acc_s = "unknown" if acc in (None, -1.0) else f"{acc:.1f} m"
        print(f"{name}")
        print(f"  stated accuracy   {acc_s}")
        print(f"  vs best pipeline  mean {offset.mean():>8.3f} m  "
              f"max {offset.max():>8.3f} m")
        print(f"  round-trip error  mean {round_trip.mean():>8.3f} m")
    return results

sites = gpd.read_file("sites.gpkg").to_crs(4326)
compare_transformations(sites, 27700)

3 pipelines, 412 points

Inverse of OSGB36 to WGS 84 (9) + British Nati
  stated accuracy   1.0 m
  vs best pipeline  mean    0.000 m  max    0.000 m
  round-trip error  mean    0.000 m

Inverse of OSGB36 to WGS 84 (5) + British Nati
  stated accuracy   3.0 m
  vs best pipeline  mean    0.842 m  max    1.204 m
  round-trip error  mean    0.001 m

Ballpark geographic offset from WGS 84 to OSGB
  stated accuracy   unknown
  vs best pipeline  mean  102.418 m  max  118.882 m
  round-trip error  mean    0.002 m

Three readings. The 7-parameter transformation is within about a metre of the grid-based one, which matches its stated 3 m accuracy — usable for mapping, not for survey. The ballpark result is 102 m out on average and 119 m at worst, which is the full magnitude of the OSGB36 offset.

The round-trip errors are the trap. All three are near zero, because each pipeline is internally consistent — the ballpark transformation inverts itself perfectly while placing every point 100 m from where it belongs. A round-trip test proves consistency, never correctness, and that distinction catches people out.

Example 3: pinning transformations for a project #

# transforms.py — one definition, checked at import
import os
os.environ.setdefault("PROJ_NETWORK", "OFF")     # before pyproj is imported

import pyproj
from pyproj import CRS, Transformer, TransformerGroup

MIN_ACCURACY_M = 5.0

REQUIRED = [
    (4326, 27700, "GPS to British National Grid"),
    (27700, 4326, "British National Grid to GPS"),
    (4326, 3857, "GPS to Web Mercator (display only)"),
]

def _build(src, dst, label):
    group = TransformerGroup(CRS.from_user_input(src), CRS.from_user_input(dst),
                             always_xy=True)
    usable = [t for t in group.transformers
              if t.accuracy is not None and 0 <= t.accuracy <= MIN_ACCURACY_M]
    if not usable:
        missing = [g.name for op in group.unavailable_operations for g in op.grids]
        raise RuntimeError(
            f"{label}: no transformation within {MIN_ACCURACY_M} m. "
            f"Best available: {group.transformers[0].description}. "
            + (f"Missing grids: {sorted(set(missing))}. "
               f"Run: projsync --file {sorted(set(missing))[0]}" if missing else ""))
    return usable[0]

TRANSFORMS = {}
for src, dst, label in REQUIRED:
    t = _build(src, dst, label)
    TRANSFORMS[(src, dst)] = t

def transform_points(lons, lats, src=4326, dst=27700):
    return TRANSFORMS[(src, dst)].transform(lons, lats)

def provenance():
    return {
        "proj": pyproj.proj_version_str,
        "pyproj": pyproj.__version__,
        "data_dir": str(pyproj.datadir.get_data_dir()),
        "network": pyproj.network.is_network_enabled(),
        "pipelines": {f"{s}->{d}": t.description for (s, d), t in TRANSFORMS.items()},
    }

if __name__ == "__main__":
    import json
    print(json.dumps(provenance(), indent=2))

{
  "proj": "9.4.0",
  "pyproj": "3.6.1",
  "data_dir": "/opt/conda/share/proj",
  "network": false,
  "pipelines": {
    "4326->27700": "Inverse of OSGB36 to WGS 84 (9) + British National Grid",
    "27700->4326": "British National Grid + OSGB36 to WGS 84 (9)"
  }
}

Building the transformers at import means a container missing OSTN15 fails on startup, with a message naming the file and the projsync command — not four hours later with coordinates that are quietly 100 m out.

Setting PROJ_NETWORK with setdefault before importing pyproj respects an explicit override while defaulting to off. It must happen before the import: PROJ reads the variable at initialisation, and setting it afterwards has no effect.

The provenance() function is what to write next to any published output, and what to compare when two machines disagree.

Explanation #

A coordinate transformation between two CRS is really two operations, and confusing them explains most of the surprises here.

The first is a projection change — converting between angles on an ellipsoid and coordinates on a plane. It is a deterministic formula with no ambiguity: the same input always gives the same output, to floating-point precision.

The second is a datum change — moving between two different models of the earth's shape and position. That is not a formula but an empirical measurement, and it is where the ambiguity lives. WGS 84 is geocentric and global; OSGB36 uses the Airy 1830 ellipsoid, fitted to Britain in the nineteenth century and positioned to match a network of triangulation stations. The relationship between them varies across the country because the historical survey has its own distortions, so it can be modelled in several ways, at several accuracies.

That is why more than one transformation exists, and why they disagree:

• A null transformation treats the datums as coincident. In Britain it is out by up to 120 m. PROJ calls this "Ballpark".
• A 7-parameter Helmert transformation applies a global rotation, translation and scale. Better — about 3 m for Britain — but it cannot capture local distortion.
• A grid-based transformation stores the measured offset at thousands of points and interpolates. OSTN15 achieves about 1 m, because it models the actual survey rather than an idealised relationship.

The grid is data, not code. It ships separately, is updated occasionally, and is the thing most likely to be missing on a fresh machine — which is exactly why the same code gives different answers in different environments.

PROJ_NETWORK turns that into a runtime dependency. With it on, PROJ fetches grids from a CDN as needed. This is genuinely helpful for interactive work and unacceptable for anything reproducible, because the result now depends on connectivity and cache state at the moment of execution. Installing proj-data and turning the network off makes a missing grid a visible failure instead.

The round-trip test is the trap worth naming. Transforming there and back returns almost exactly the original coordinates under every pipeline, including the ballpark one, because each is internally consistent. A perfect round trip proves the transformation is invertible, not that it is right. The only real check is Transformer.description and accuracy — which is why logging them alongside results is the practice that matters most.

Finally, note that the accuracy required depends entirely on the use. A 1:50,000 map has a line width worth about 25 m on the ground, so a 3 m transformation is beyond adequate. A cadastral boundary, a buried utility, or anything with legal weight needs the grid. Knowing which situation you are in is the decision; everything above is how to act on it.

Edge cases or notes #

• "Ballpark" means PROJ found no proper transformation. Error can exceed 100 m; never publish without saying so.
• accuracy = -1.0 means unknown, not excellent. Filter with 0 <= accuracy <= limit.
• PROJ_NETWORK is read at import. Set it before import pyproj, or use pyproj.network.set_network_enabled.
• A perfect round trip proves consistency, not correctness. Ballpark round-trips perfectly too.
• always_xy=True forces lon/lat order. Without it, EPSG:4326's declared axis order is lat/lon, and results silently swap.
• TransformerGroup.unavailable_operations names the grid files you are missing.
• projsync --area-of-use "United Kingdom" fetches only the relevant grids instead of 600 MB.
• ETRS89 and WGS 84 are diverging by about 2.5 cm per year as the Eurasian plate moves — roughly 0.9 m since 1989.
• +towgs84= in a PROJ string is a 7-parameter shift, less accurate than a grid and often the reason results differ.
• GeoPandas' to_crs uses the same machinery, so everything here applies to it.

Internal links #

• Coordinate reference systems explained for Python GIS — datums and ellipsoids
• Projected vs geographic CRS — the projection half of a transformation
• How to reproject spatial data in Python (GeoPandas) — the everyday operation
• Reproducible GIS environments explained — pinning PROJ and its data
• How to work with a custom or local CRS in Python — defining a CRS with a datum
• How to record run metadata and data lineage in a GIS pipeline — storing the pipeline with results
• EPSG codes explained: how to choose the right CRS in Python — what a code specifies
• How to measure distance accurately in Python — where a metre matters

FAQ #

Why do I get different coordinates on different machines? #

Different transformation pipelines, because the grid files available differ. Print Transformer.description on both — one is probably using OSTN15 and the other falling back to ballpark.

What does "Ballpark" mean in a transformation description? #

PROJ found no proper datum transformation and treated the two datums as coincident. In Britain that is an error of up to about 120 m. Install the grid file.

How do I install the transformation grids? #

conda install -c conda-forge proj-data for the full set, or projsync --area-of-use "United Kingdom" for just the relevant ones.

Should PROJ_NETWORK be on or off? #

Off for anything reproducible. With it on, grids are fetched at runtime, so results depend on connectivity and cache state. Install proj-data instead.

My round trip is accurate to a millimetre. Does that mean the transformation is right? #

No. Every pipeline is internally consistent and inverts itself cleanly, including the ballpark one. Check description and accuracy, not the round trip.

How accurate does my transformation need to be? #

For a 1:50,000 map, a 3 m Helmert is more than enough. For cadastral, utility or legal work, use the grid-based transformation — about 1 m — and record which one you used.

What is always_xy=True for? #

It forces longitude-then-latitude ordering. EPSG:4326 formally declares latitude first, so without it your coordinates can silently swap.