Spatial Epi Public Health GIS

Handling Datum Shifts Between NAD83 Realizations

This guide, part of Coordinate Reference Systems for Public Health, addresses a source of silent coordinate error that most surveillance pipelines carry without knowing it: treating every dataset labelled “NAD83” as though it used the same reference frame.

Problem Context & Constraints

NAD83 has been re-realized several times since its first adjustment, as geodetic networks improved and satellite positioning replaced ground surveys. The original 1986 adjustment, the high-accuracy regional networks of the 1990s, the 2007 national readjustment and the 2011 realization all carry the NAD83 name, and the coordinates of the same physical point differ between them — by tens of centimetres in most places, by a metre or more in some regions and older layers.

For a lot of public health work, a metre does not matter. It does matter in three places. Boundary assignment for points near a line: a case point half a metre from a tract boundary can land on either side depending on which realization the boundary file and the geocoder each assumed. Overlays of fine layers: parcel polygons in one realization overlaid on building footprints in another produce slivers and misassignments at exactly the resolution those layers exist for. Change detection: comparing coordinates from two geocoding vintages that used different realizations shows movement that is entirely a datum artefact, as the drift checks in Detecting Coordinate Drift Between Geocoding Vintages would report.

Three constraints shape the fix. Identify the realization of every layer, which is harder than it sounds because many files declare only “NAD83” without specifying which. Transform with the published grid-based method, not a ballistic no-op, and confirm the grid is actually installed. And decide whether the shift matters for the analysis at hand, because transforming everything to the latest realization is correct but not always necessary.

One Point, Several NAD83 Positions A single physical case location plotted with its coordinates as expressed in four NAD83 realizations: the original 1986 adjustment, a regional high-accuracy network, the 2007 readjustment and the 2011 realization. The positions cluster within about a metre of one another. A tract boundary runs through the cluster, so the point falls in one tract under two realizations and in the neighbouring tract under the other two. Every layer labelled NAD83 must have its realization identified before points and boundaries are combined. "NAD83" can place one address on either side of a line tract A tract B 1986 regional network 2007 2011 about one metre spread identify each layer's realization before combining points with boundaries

Prerequisites

  • python 3.11
  • pyproj 3.6.1 with the transformation grids for your region installed, geopandas 0.14.4
  • Metadata or provenance for each layer sufficient to identify its realization
  • The pipeline’s target realization, usually the most recent one used by your authoritative boundary source
  • A handful of surveyed control points, or published coordinates of known monuments, for verification

Step-by-Step Solution

Record the realization of every input, build explicit transformations between realizations with the grid-based method, verify the grid was used, and measure the shift on your own data before deciding whether to apply it.

# Identify and transform between NAD83 realizations with a verified grid method.
# Pinned: pyproj==3.6.1, geopandas==0.14.4, numpy==1.26.4
import logging

import numpy as np
import geopandas as gpd
from pyproj import CRS
from pyproj.transformer import TransformerGroup

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

# Geographic CRSs for the realizations. A file that says only "NAD83" is EPSG:4269,
# which is ambiguous about realization — resolve it from provenance, not the label.
REALIZATIONS = {
    "NAD83(1986)": "EPSG:4269",
    "NAD83(HARN)": "EPSG:4152",
    "NAD83(NSRS2007)": "EPSG:4759",
    "NAD83(2011)": "EPSG:6318",
}


def layer_realization(declared_crs: CRS, provenance_hint: str | None) -> str:
    """Resolve a layer's realization. The declared CRS alone is often insufficient."""
    name = declared_crs.name
    for key in REALIZATIONS:
        if key in name:
            return key
    if provenance_hint in REALIZATIONS:
        log.warning("CRS declares only '%s'; using provenance hint %s", name, provenance_hint)
        return provenance_hint
    raise ValueError(f"cannot determine NAD83 realization for '{name}' — "
                     "record it from the data source before using this layer")


def best_transformer(src: str, dst: str):
    """Choose the most accurate AVAILABLE operation and refuse a silent no-op."""
    group = TransformerGroup(CRS(REALIZATIONS[src]), CRS(REALIZATIONS[dst]), always_xy=True)
    if group.unavailable_operations:
        log.warning("%s -> %s: better operations need grids that are not installed: %s",
                    src, dst, [o.name for o in group.unavailable_operations][:3])
    t = group.transformers[0]
    desc = t.description.lower()
    if "ballpark" in desc or "null" in desc:
        raise RuntimeError(f"{src} -> {dst} resolved to '{t.description}': "
                           "a no-op, install the transformation grid")
    log.info("%s -> %s via '%s' (stated accuracy %s m)", src, dst, t.description, t.accuracy)
    return t


def shift_between(points: gpd.GeoDataFrame, src: str, dst: str) -> dict:
    """How far does moving between realizations actually shift THIS data?"""
    t = best_transformer(src, dst)
    x0 = points.geometry.x.to_numpy()
    y0 = points.geometry.y.to_numpy()
    x1, y1 = t.transform(x0, y0)
    # Convert the lon/lat difference to metres at each latitude.
    dy = (y1 - y0) * 111_320.0
    dx = (x1 - x0) * 111_320.0 * np.cos(np.radians(y0))
    d = np.hypot(dx, dy)
    out = {"median_m": float(np.median(d)), "max_m": float(d.max()),
           "p95_m": float(np.percentile(d, 95))}
    log.info("%s -> %s shift on this data: median %.3f m, p95 %.3f m, max %.3f m",
             src, dst, out["median_m"], out["p95_m"], out["max_m"])
    return out

The refusal to accept a no-op transformation is the most important line. When the grid for a realization change is not installed, the transformation engine does not fail — it falls back to a ballpark operation that leaves coordinates unchanged and reports success. A pipeline that believes it transformed a layer from one realization to another, when it actually passed the coordinates through untouched, has the datum error it thought it had fixed and no signal that anything is wrong.

Measuring the shift on the agency’s own data turns an abstract geodetic concern into a decision. A median shift of 3 cm in a state whose tract boundaries are generalised to tens of metres is not worth a pipeline change. A median of 90 cm in a state whose parcel-level geocodes feed a tract assignment near boundaries probably is.

The Silent No-Op When a Grid Is Missing Two outcomes of the same transformation request between realizations. With the regional grid installed, the engine applies a grid-based shift and moves each point by its correct local offset of several decimetres. Without the grid, the engine falls back to a ballpark operation that leaves coordinates unchanged, and reports success. The output files are indistinguishable in structure, and only a check on the chosen operation's description or a comparison with control points reveals that nothing was transformed. A missing grid does not fail — it quietly does nothing Grid installed operation: grid-based shift points move 0.3 to 0.9 m correct Grid missing operation: ballpark fallback points move 0.0 m reported as success inspect the chosen operation and refuse a ballpark or null result

Validation & Edge Cases

Layers that declare only “NAD83”. The generic geographic code does not say which realization, and different producers use it to mean different things. Resolve from the source’s documentation or contact the producer; never assume.

Projected systems built on different realizations. A state plane zone exists in several realizations with different EPSG codes. Two files in “the same state plane zone” can differ by a metre if their underlying realizations differ.

Web services that reproject on the fly. A map service that serves data in a web projection has already chosen a transformation, possibly a no-op. Request data in its native system where possible and do the transformation yourself with a verified method.

Control point verification. Published coordinates of geodetic monuments in each realization let you check a transformation end to end. Transform the monument’s older coordinates and compare with its published newer ones; agreement within the stated accuracy confirms the grid is in use.

Time-dependent frames. The newer terrestrial reference frames that will replace NAD83 are tied to an epoch. Coordinates without an epoch are ambiguous at the centimetre level in those frames, which matters for future-proofing a pipeline even where it does not yet affect surveillance results.

WARNING crs.nad83 CRS declares only 'NAD83'; using provenance hint NAD83(HARN)
INFO crs.nad83 NAD83(HARN) -> NAD83(2011) via 'NAD83(HARN) to NAD83(2011) (NADCON5)' (stated accuracy 0.05 m)
INFO crs.nad83 NAD83(HARN) -> NAD83(2011) shift on this data: median 0.412 m, p95 0.688 m, max 0.931 m
Does the Shift Matter Here? A decision based on three measured quantities. If the median realization shift on the agency's own data is much smaller than both the geocode positional accuracy and the generalisation of the boundary file, the shift does not change any assignment and can be documented rather than applied. If the shift is comparable to either, which occurs with parcel-level geocodes and detailed boundaries, the layers must be brought to a common realization before any point-in-polygon step. Compare the shift with the precision it could disturb measured shift median 0.4 m geocode accuracy rooftop about 5 m boundary detail parcel-aligned, sub-metre shift far below both document, need not apply shift near either unify realizations first with parcel-aligned boundaries, a 0.4 m shift is near the boundary detail: unify

Compliance Notes

Record for every layer its NAD83 realization and how it was determined — declared, documented or inferred from provenance — along with the target realization, the transformation operation actually used and its stated accuracy, and the measured shift on your data. Where a layer’s realization could not be determined, record that as a known limitation rather than assigning one silently.

Include the realization in published metadata. A partner combining your case points with their own parcel layer needs to know which realization your coordinates are in, and “NAD83” alone does not tell them. The one line it takes to state the realization prevents a class of sub-metre misalignment that is otherwise invisible until someone overlays two layers at building scale.