Production geodesy: PROJ pipelines, NTv2 grids, time-dependent frames

You saw the geometric foundation in Week 3: ellipsoid vs geoid, why GPS reports one and maps want the other, the EGM2008 model. This week is the engineering layer that makes that foundation usable in production — where you're ingesting data from civil-survey teams (often NAD83 or OSGB36), aviation feeds (WGS84), tide-gauge data (orthometric to a national vertical datum like NAVD88), and satellite-derived geometry (ITRF) all into the same database, and you need them to agree to centimeters.

Quick refresher (Week 3 said this)

For any point you carry three heights: h (ellipsoidal, what GPS reports), N (geoid undulation from a model like EGM2008), and H = h − N (orthometric, what maps and humans use). The WGS84 ellipsoid has equatorial radius 6,378,137 m and flattening 1/298.257223563. The geoid undulates by ±100 m relative to that ellipsoid. None of that is new. What is new this week:

PROJ pipelines: multi-step transforms in one composition

Modern PROJ (8+) exposes transformation pipelines — a string-grammar that composes any chain of operations: ellipsoid swap, datum shift, projection, vertical-datum transform, plate-motion correction. Each step is a verb; PROJ assembles them into one optimized transform:

# 3D pipeline: NAD27 lat/lon/orthometric → WGS84 lat/lon/ellipsoidal
pipeline=+proj=pipeline
  +step +proj=axisswap +order=2,1                  # PROJ wants lon,lat,h
  +step +inv +proj=hgridshift +grids=us_noaa_conus.tif  # NAD27 → NAD83 via NADCON
  +step +proj=vgridshift +grids=us_noaa_g2018u0.tif      # NAVD88 → ellipsoidal
  +step +proj=axisswap +order=2,1                  # back to lat,lon,h

Each grid file is a small binary correction surface (NTv2 for horizontal, VRG for vertical). PROJ ships these via the PROJ data network — invoke once, and PROJ downloads + caches the grids you need. Tells you to do this with --download on first use.

NTv2 grids: the datum-shift workhorse

A constant 7-parameter Helmert transformation is good to a meter or two. For sub-meter, every national survey agency publishes a grid-shift file: NAD27→NAD83 (NADCON in the US; the original NTv2 format was developed by NRCan for Canada), HARN/HPGN refinements, OSGB36→ETRS89 (OSTN15 in the UK), GDA94→GDA2020 (Australia's adjustment for plate-tectonic drift). These are dense correction surfaces that capture local distortion better than any global formula can — and despite being called "NTv2 grids" by convention (PROJ accepts the format universally), the underlying datum-shift methodology differs by country.

The trap: if you skip the grid shift and treat NAD83 as WGS84 (a tempting "they're close" shortcut), you accumulate 1–2 m of horizontal error in CONUS — and that error compounds when you stack data sources. The grid file is small (~500 KB to a few MB per country). Use it.

Time-dependent reference frames: ITRF

The deepest cut: WGS84 itself isn't constant. The most recent realization is ITRF2020 (International Terrestrial Reference Frame, IERS), refined every ~6 years. Tectonic plates move at 1–10 cm/year, so a coordinate measured "in WGS84" in 2010 isn't the same physical point as one measured "in WGS84" in 2026 — even at the same surveyed monument.

For sub-cm-class work — InSAR time series (Week 23), GNSS time series, satellite-conjunction analysis — you have to specify the ITRF realization (ITRF2008, ITRF2014, ITRF2020) and the epoch (e.g. ITRF2014@2020.5 means "ITRF2014, valid at mid-2020"). PROJ lets you encode this explicitly:

from pyproj import Transformer

# Transform a coordinate from ITRF2008@2005 to ITRF2020@2024
t = Transformer.from_pipeline(
    "+proj=pipeline "
    "+step +proj=cart +ellps=GRS80 "
    "+step +proj=helmert +convention=position_vector "
    "  +x=0.0014 +y=0.0009 +z=-0.0014 +s=-0.00042 "  # ITRF2008→ITRF2020 epoch params
    "  +rx=0.000008 +ry=0.000000 +rz=0.000000 "
    "  +t_epoch=2010.0 +drx=... "  # epoch-rate parameters
    "+step +inv +proj=cart +ellps=GRS80"
)

Most application code never needs this level of detail. The exceptions are precise positioning (RTK GPS, GNSS networks), geodynamic monitoring (volcano deformation, plate boundaries), and any cross-decade time-series. For those, the wrong ITRF realization shows up as a 5–20 cm drift in your residuals that defies all other explanation.

Datum conflicts in multi-source pipelines

The class of production bug that destroys a Monday: two feeds arrive in your pipeline, both labelled "WGS84", and they disagree by 30 cm. Why? One feed is actually EPSG:4326 (geographic WGS84) and the other is EPSG:4979 (3D WGS84 with ellipsoidal height). Or one is ITRF2008-epoch-2005 and the other is ITRF2014-epoch-2020. Or one snuck NAVD88 vertical-datum elevations into a field labelled "elevation_m" without saying so.

The hygiene that prevents this:

• Always store the CRS explicitly — EPSG code or full PROJ string, never an English-language name. Postgres has a srid column on every geometry; use it. Avoid SRID=0 ("unknown") in production.
• Reject mixed-srid writes at the database layer with a check constraint.
• Stamp every external ingest with a provenance record: source, ingestion timestamp, original CRS, declared epoch (if applicable). PROJ-pipeline corrections become a documented step.
• Reconciliation tests: pull a known landmark from each upstream feed; compute pairwise distances; alarm if any pair exceeds the source-specified precision.

EGM2008 corrections in pyproj

The vertical-datum transform you sketched in Week 3 — GPS ellipsoidal → orthometric — collapses into a single composed CRS in pyproj. The +5773 at the end of the target CRS is EPSG's code for "EGM2008 height":

from pyproj import Transformer

# Ellipsoid height (from GPS) → orthometric height (above EGM2008 geoid)
to_orthometric = Transformer.from_crs('EPSG:4979', 'EPSG:4326+5773', always_xy=True)
lon, lat, h = -118.6, 36.5, 4421  # Mt. Whitney, GPS-measured
lon, lat, H = to_orthometric.transform(lon, lat, h)
print(f"Orthometric height: {H:.1f} m")  # ≈ 4421 - N(36.5N, -118.6) ≈ 4393 m

For higher accuracy in the US specifically, swap EGM2008 for GEOID18 (NGS's regional model, ~2 cm accuracy in CONUS) by including the us_noaa_g2018u0.tif grid in a PROJ pipeline rather than the global EGM2008 grid. For mainland Europe, swap for the EVRF realisation (nkg_eur_egm2008 or per-country grids like the German de_adv_GCG2016). The point: never assume one geoid model fits everywhere — for sub-meter vertical work, the regional model always beats the global one.

Why this matters for space GIS

For most space-domain applications, ellipsoidal height is fine. The exceptions:

• Comparing satellite-altitude-derived sub-satellite positions to terrain elevation (orthometric).
• Computing line-of-sight blocked by terrain: the terrain DEM is orthometric; the satellite altitude is ellipsoidal. You must convert one to match.
• Reporting "altitude" to non-technical users — they expect mean-sea-level. Use orthometric and label it.
• Precise positioning of ground stations (sub-meter), where the 100 m geoid offset would otherwise be a meaningful error source.

The lab

You'll take a set of GPS-measured points along the Sierra Nevada (provided), compute their ellipsoidal heights, apply EGM2008 geoid correction to get orthometric heights, and compare with USGS-published orthometric heights at the same coordinates. The agreement should be within ±2 m if the GPS measurements were of survey quality. By the end, you'll never confuse the three heights again.