Spatial Epi Public Health GIS

Pinning PROJ and GDAL Versions for Reproducible Results

This guide, part of Testing Geospatial Analysis Code, deals with the dependency that most often changes a surveillance pipeline’s output without anyone changing the pipeline: the coordinate transformation engine and its data.

Problem Context & Constraints

The projection engine underneath almost every Python geospatial stack does more than arithmetic. For a given pair of coordinate reference systems it chooses a transformation pipeline — which datum shift, which grid file, which method — from a database of candidates, based on what is available and what the area of use is. That choice depends on three things that change independently: the library version, the version of its database, and which grid files happen to be installed.

A change in any of them can move coordinates. Installing a new geodetic grid can switch a transformation from a ballistic approximation to a grid-based one and shift points by a metre or more. A database update can change the preferred transformation between two datums. A new library default can alter how an ambiguous request is resolved. None of these raise an error. All of them change which tract some addresses land in.

Three constraints follow. Pin the library and its data together, because pinning one without the other leaves the dominant source of variation free. Record the transformation actually used, not just the systems requested, because the same request can resolve differently on two machines. And treat an upgrade as a change to the analysis, validated by a differential comparison against the pinned version.

One Request, Several Possible Answers A single request to transform coordinates from one datum to another enters a resolver. The resolver consults three independently versioned inputs: the library version, the transformation database version, and the set of installed grid files. Depending on those, it chooses one of several transformation pipelines — a grid-based shift, a fourteen-parameter Helmert, or a ballistic no-op — each producing coordinates that differ by up to a few metres. The note records that pinning the library alone leaves the database and grids free to change the answer. The engine chooses the math; three versioned inputs decide what it chooses Transform request datum A to datum B Pipeline resolver picks the best available library version database grid files grid-based shift parameterised Helmert ballistic approximation answers differ by up to metres Pin all three, and record which pipeline was actually used

Prerequisites

  • python 3.11
  • pyproj 3.6.1 with its bundled or system projection library, rasterio 1.3.10, pandas 2.2.2
  • A lock file or container image mechanism that pins native libraries as well as Python packages
  • A regression fixture of coordinates spanning the jurisdiction, including its edges
  • Access to the transformation grids the jurisdiction requires

Step-by-Step Solution

Capture the full environment at runtime, record the resolved transformation for every pair the pipeline uses, and compare two environments on a fixture before accepting an upgrade.

# Record the projection environment and the transformation actually used; diff upgrades.
# Pinned: pyproj==3.6.1, rasterio==1.3.10, pandas==2.2.2, numpy==1.26.4
import json
import logging

import numpy as np
import pandas as pd
import pyproj
from pyproj import CRS, Transformer
from pyproj.transformer import TransformerGroup

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("env.proj")


def projection_environment() -> dict:
    """Everything that can change a coordinate, captured at runtime."""
    import rasterio
    env = {
        "pyproj": pyproj.__version__,
        "proj": pyproj.proj_version_str,
        "proj_database": pyproj.database.get_database_metadata("DATABASE.LAYOUT.VERSION.MAJOR"),
        "epsg_version": pyproj.database.get_database_metadata("EPSG.VERSION"),
        "gdal": rasterio.__gdal_version__,
        "data_dir": pyproj.datadir.get_data_dir(),
        "network_enabled": pyproj.network.is_network_enabled(),
    }
    log.info("projection environment: %s", json.dumps(env, sort_keys=True))
    return env


def resolved_transformation(src: str, dst: str, area=None) -> dict:
    """Which pipeline the engine will actually use, and what it could have used."""
    group = TransformerGroup(CRS(src), CRS(dst), always_xy=True, area_of_interest=area)
    chosen = group.transformers[0]
    missing = [op.name for op in group.unavailable_operations]
    out = {"src": src, "dst": dst, "chosen": chosen.description,
           "accuracy_m": chosen.accuracy, "pipeline": chosen.to_proj4() if hasattr(
               chosen, "to_proj4") else str(chosen),
           "better_but_unavailable": missing}
    if missing:
        log.warning("%s -> %s: more accurate operations unavailable (grids missing): %s",
                    src, dst, missing[:3])
    log.info("%s -> %s resolved to '%s' (accuracy %s m)",
             src, dst, chosen.description, chosen.accuracy)
    return out


def transform_fixture(points: pd.DataFrame, src: str, dst: str) -> pd.DataFrame:
    t = Transformer.from_crs(src, dst, always_xy=True)
    x, y = t.transform(points["lon"].to_numpy(), points["lat"].to_numpy())
    return pd.DataFrame({"pid": points["pid"], "x": x, "y": y})


def compare_environments(before: pd.DataFrame, after: pd.DataFrame,
                         tolerance_m: float = 0.05) -> dict:
    """Differential check for an upgrade: did any fixture point move materially?"""
    m = before.merge(after, on="pid", suffixes=("_old", "_new"))
    d = np.hypot(m["x_new"] - m["x_old"], m["y_new"] - m["y_old"])
    out = {"n": int(len(m)), "max_shift_m": float(d.max()),
           "p99_shift_m": float(np.percentile(d, 99)),
           "n_over_tolerance": int((d > tolerance_m).sum()),
           "tolerance_m": tolerance_m}
    (log.info if out["n_over_tolerance"] == 0 else log.error)(
        "upgrade diff: max=%.4f m p99=%.4f m over_tol=%d of %d",
        out["max_shift_m"], out["p99_shift_m"], out["n_over_tolerance"], out["n"])
    return out

Recording the resolved transformation, not just the requested systems, is what makes a run reproducible in fact rather than in principle. Two runs that both asked for the same source and target systems may have used different operations; the provenance row that says which operation was chosen, with its stated accuracy, is the only evidence that they did or did not.

The warning about better operations being unavailable is worth treating as a finding. It means the engine fell back to a less accurate transformation because a grid file was not installed — which is common in minimal container images and is exactly how two environments running the same code produce coordinates a metre apart.

What an Upgrade Diff Looks Like A histogram of the distance each fixture point moved between the pinned environment and a candidate upgrade. Most points moved by less than a millimetre, reflecting floating-point differences only. A distinct group in the south-east of the jurisdiction moved by about one point two metres, because the upgrade installed a newer geodetic grid covering that region. The tolerance of five centimetres is marked, and the moved group lies entirely beyond it, so the upgrade must be treated as a change to the analysis rather than accepted silently. Most points did not move; a whole region moved by a metre 5 cm tolerance floating point south-east region new grid installed <1 mm 5 cm ~1.2 m distance each fixture point moved between environments

Validation & Edge Cases

Network-enabled grid fetching. Some configurations download grids on demand, which means the transformation used depends on network availability at run time. Disable it for production and install the required grids explicitly.

Pinning Python packages but not native libraries. A lock file that pins the Python binding but lets the system library float leaves the most important component unpinned. Pin the container image or use wheels that bundle the library, and record the native version at runtime.

Fixtures that do not cover the grid regions. A fixture concentrated in one metropolitan area will miss a shift confined to a region covered by a newly installed grid. Spread fixture points across the whole jurisdiction, including its edges.

Axis order. A change in how an ambiguous system’s axis order is interpreted swaps latitude and longitude for every point, which the differential comparison catches as an enormous shift. Always request an explicit axis order.

Developer laptops that differ from production. An analyst reproducing a result on a laptop with a system-installed library of a different version will get slightly different coordinates and may conclude the pipeline is wrong. Publish the container image used in production and make it the supported way to reproduce any run, rather than documenting the version numbers and hoping.

Raster reprojection and resampling defaults. The raster library’s resampling behaviour and its handling of nodata at edges also change between versions, and an exposure surface reprojected under two versions can differ at every cell boundary. Include a small raster in the regression fixture and compare it cell by cell alongside the point transformations.

Epoch-dependent coordinate systems. Dynamic datums tie coordinates to an observation epoch, and a transformation that ignores the epoch is only approximately right. Where the pipeline uses such a system, record the epoch with the data, and include it in the environment check so an upgrade that starts honouring epochs is detected rather than silently changing results.

Accepting a real improvement. Sometimes the upgrade is more correct — a better grid is genuinely more accurate. Accept it deliberately, rerun the regression fixture, update stored outputs in a reviewed change, and record the date so trend analyses can account for it.

INFO env.proj projection environment: {"epsg_version": "v11.004", "gdal": "3.8.5", "network_enabled": false, "proj": "9.4.0", "pyproj": "3.6.1"}
WARNING env.proj EPSG:4269 -> EPSG:6318: more accurate operations unavailable (grids missing): ['NAD83 to NAD83(2011) (NADCON5)']
ERROR env.proj upgrade diff: max=1.2184 m p99=1.1702 m over_tol=318 of 4000

The recommended practice — one published container image as the supported way to reproduce a run — replaces a set of version numbers that drift independently:

One Image, Every Environment On the left, three environments each assembled independently: production, an analyst's laptop and a reviewer's workstation, each with a slightly different projection library, database and grid set, so the same code yields coordinates that differ between them. On the right, one pinned container image containing the library, its database and the required grids is used by all three, so the coordinates agree everywhere and a reproduction is a true reproduction. Reproduce with the production image, not with a version list Assembled separately production: library 9.4, grids A+B laptop: library 9.2, grid A only reviewer: library 9.5, network grids three sets of coordinates One pinned image library + database + grids one digest production, laptop, reviewer one set of coordinates record the image digest in every run's provenance

Compliance Notes

Write the projection environment and every resolved transformation into the run’s provenance row. For a published product, include the transformation description and its stated accuracy in the metadata, because a recipient reprojecting your data needs to know which operation you used to reproduce your coordinates.

Treat an engine upgrade with the same governance as a method change. It can move case locations across area boundaries, which can move counts across suppression thresholds and change cluster membership. The differential comparison report, the decision to accept or reject, and the reviewer belong in the change record — and where the upgrade is accepted mid-series, the date belongs in the metadata of every trend product that spans it.