Geospatial Schema Standardization & Mapping Automation

Chaining Vertical and Horizontal Datum Shifts

Elevation data is where datum handling quietly breaks. A LiDAR tile in NAD27 horizontal / NGVD29 vertical must land in NAD83 / NAVD88 for a modern flood model, and those are two independent shifts — a horizontal NADCON grid and a vertical GEOID-derived grid — that must compose into one ordered 3D transform. Apply them in the wrong order, or apply the vertical shift with the pre-shifted horizontal position, and heights drift by decimetres exactly where a floodplain boundary is decided. This procedure composes the compound transform explicitly and budgets its error. It extends the routing contract in Datum Transformation Fallback Chains, the parent stage inside CRS Normalization & Sync that treats horizontal and vertical passes as separate, independently gated operations.

The steps follow configure (Step 1, compound CRS), execute (Steps 2-3, pipeline and transform), validate (Step 4, accuracy accumulation), and log (Step 5, lineage). Unlike a purely horizontal shift, a 3D chain has two grids that can go missing and two accuracy figures that add — both realities drive the design below.

Ordered Compound 3D Datum Transform Pipeline A left-to-right pipeline. The input carries NAD27 horizontal coordinates and NGVD29 heights. Stage one applies a NADCON5 horizontal grid shift to x and y. Stage two applies a GEOID18 vertical shift to z, sampled at the horizontal position produced by stage one. The output is NAD83 horizontal with NAVD88 orthometric heights. A note below shows the horizontal and vertical component accuracies summing into a combined 3D error budget. Source 3D NAD27 / NGVD29 1 · Horizontal NADCON5 grid → x, y ±0.10 m 2 · Vertical GEOID18 → z at (x′,y′) ±0.05 m Output 3D NAD83 / NAVD88 combined 3D budget = √(horizontal² + vertical²) ≈ 0.11 m → gate against tolerance

Prerequisites checklist

Step 1: Declare the compound CRS pair

A 3D transform needs 3D CRSes on both ends. Build compound CRSes that pair a horizontal component with a vertical one so PROJ sees all three axes and can route both shifts. Guessing the vertical component is the most common upstream defect; resolve it explicitly the way Projection Normalization Workflows resolves ambiguous horizontal codes.

python
# pyproj >= 3.6 — Python 3.10+
from pyproj import CRS

# Compound source: NAD27 horizontal (4267) + NGVD29 height (5702)
source = CRS.from_user_input("EPSG:4267+5702")
# Compound target: NAD83(2011) horizontal (6318) + NAVD88 height (5703)
target = CRS.from_user_input("EPSG:6318+5703")

assert source.is_compound and target.is_compound, "both ends must carry a vertical axis"
assert len(source.axis_info) == 3, "expected lat, lon, height on the source"

Declaration rules:

  • Use the horizontal+vertical compound form ("EPSG:4267+5702"); a bare horizontal code silently drops the height axis.
  • Verify is_compound and a 3-element axis_info before proceeding — a 2D CRS will pass through z untouched and corrupt it invisibly.
  • Never assume NGVD29 or NAVD88 from the horizontal datum; they are independent and must come from source metadata.

Step 2: Compose the ordered pipeline

Order matters. The vertical GEOID shift is sampled at a horizontal position, so the horizontal shift must apply first; the GEOID grid is then read at the already-corrected (x', y'). Let PROJ build the compound operation rather than hand-writing two transforms, and inspect the resulting pipeline to confirm both steps are present.

python
# pyproj >= 3.6 — Python 3.10+
from pyproj.transformer import TransformerGroup

group = TransformerGroup(source, target, always_xy=True)
if not group.transformers:
    raise RuntimeError("no compound 3D path; check GEOID/NADCON grid availability")

pipeline = group.transformers[0]
proj_string = pipeline.to_proj4() or pipeline.definition

# A valid compound pipeline references BOTH a horizontal and a vertical grid step.
assert "hgridshift" in pipeline.definition or "gridshift" in pipeline.definition, \
    "horizontal grid step missing"
assert "vgridshift" in pipeline.definition, "vertical GEOID step missing"

Composition rules:

  • Trust PROJ’s compound operation ordering; it applies +proj=hgridshift before +proj=vgridshift so the geoid is read at the shifted horizontal location.
  • If the pipeline lacks a vgridshift step, the GEOID grid is missing — apply the missing-grid recovery from Handling Missing NTv2 Grid Files in Production, which works identically for vertical grids.
  • Do not hand-concatenate two independent Transformer objects; a manual chain re-samples the geoid at the pre-shift position and introduces a systematic height bias.

Step 3: Transform full 3D coordinates

Reproject x, y, and z in one call so the vertical result reflects the corrected horizontal position. The 2D GeoDataFrame.to_crs ignores z; drive the transformer directly on all three ordinates.

python
# pyproj >= 3.6, geopandas >= 0.14 — Python 3.10+
import geopandas as gpd
from shapely import get_coordinates, set_coordinates


def transform_3d(gdf: gpd.GeoDataFrame, pipeline) -> gpd.GeoDataFrame:
    """Apply a compound 3D transform to xyz coordinates, preserving z."""
    out = gdf.copy()
    geoms = out.geometry.to_numpy()
    for i, geom in enumerate(geoms):
        coords = get_coordinates(geom, include_z=True)  # (N, 3)
        x, y, z = coords[:, 0], coords[:, 1], coords[:, 2]
        xt, yt, zt = pipeline.transform(x, y, z)         # all three together
        geoms[i] = set_coordinates(geom, list(zip(xt, yt, zt)))
    out["geometry"] = geoms
    out = out.set_crs(target, allow_override=True)
    return out

Transform rules:

  • Pass z into pipeline.transform; omitting it silently leaves heights on the source vertical datum.
  • Use include_z=True when extracting coordinates — a 2D read discards the height you are trying to shift.
  • Reject any feature whose z is NaN after transform; it fell outside the GEOID grid’s coverage and cannot be trusted.

Step 4: Accumulate and gate the error budget

The horizontal and vertical shifts each carry an accuracy; the compound 3D error is their combination, not either alone. Sum them (root-sum-square for independent components) and gate against tolerance. The threshold matrices are standardized in Unit Conversion & Tolerance Thresholds; this step enforces them for the compound case.

python
# pyproj >= 3.6 — Python 3.10+
import math
from dataclasses import dataclass


@dataclass(frozen=True)
class Budget3D:
    horizontal_m: float
    vertical_m: float

    @property
    def combined_m(self) -> float:
        return math.hypot(self.horizontal_m, self.vertical_m)


def gate_budget(budget: Budget3D, tol_h_m: float, tol_v_m: float) -> None:
    if budget.horizontal_m > tol_h_m:
        raise ValueError(f"horizontal {budget.horizontal_m:.3f} m > {tol_h_m} m")
    if budget.vertical_m > tol_v_m:
        raise ValueError(f"vertical {budget.vertical_m:.3f} m > {tol_v_m} m")
    # Combined budget is informational for the audit record, gated per-component above.

Budgeting rules:

  • Gate the horizontal and vertical components separately — a passing combined figure can still hide a vertical failure that a flood model cares about.
  • Treat the components as independent and combine with hypot; do not simply add, which overstates the budget.
  • When either component’s accuracy is -1 (unknown), reject under strict mode exactly as the parent Datum Transformation Fallback Chains does.

Step 5: Log the compound lineage

A compound transform needs a richer lineage row than a 2D one: both operation names and the GEOID model, so an auditor can reproduce the exact heights.

python
# pyarrow >= 14 — Python 3.10+
import datetime as dt
import pyarrow as pa
import pyarrow.parquet as pq


def log_compound_lineage(path: str, pipeline, budget: "Budget3D") -> None:
    row = {
        "operation": pipeline.description,
        "horizontal_grid": "us_noaa_nadcon5",
        "vertical_grid": "us_noaa_g2018",
        "horizontal_accuracy_m": budget.horizontal_m,
        "vertical_accuracy_m": budget.vertical_m,
        "combined_accuracy_m": budget.combined_m,
        "timestamp_utc": dt.datetime.now(dt.UTC).isoformat(),
    }
    pq.write_to_dataset(pa.Table.from_pylist([row]), root_path=path)

Route these rows to the append-only store described in Geospatial Compliance Reporting & Audit Trails; naming the GEOID model is mandatory to defend orthometric heights in an ISO 19115 data-quality review.

Verification

Prove the height axis actually moved and by a plausible amount:

python
out = transform_3d(gdf, pipeline)

# 1. Output is on the compound target.
assert out.crs.is_compound and out.crs.equals(target)

# 2. Heights changed — NGVD29->NAVD88 differs by roughly a decimetre to a metre in CONUS.
src_z = get_coordinates(gdf.geometry.iloc[0], include_z=True)[0, 2]
out_z = get_coordinates(out.geometry.iloc[0], include_z=True)[0, 2]
assert abs(out_z - src_z) > 0.01, "vertical datum shift did not apply — check the pipeline"

# 3. No height fell outside GEOID coverage.
assert not any(
    math.isnan(get_coordinates(g, include_z=True)[:, 2]).any() for g in out.geometry
)

On the CLI, echo "-96 40 100" | cs2cs EPSG:4267+5702 EPSG:6318+5703 should return a third ordinate that differs from 100, confirming the vertical step ran; an unchanged height means the vgridshift step is missing.

Troubleshooting

Symptom Likely cause Fix
Heights unchanged after transform 2D to_crs used, or a bare horizontal CRS dropped the z axis Use the compound "EPSG:h+v" form and drive pipeline.transform(x, y, z) (Steps 1, 3).
Vertical off by a constant ~0.1–1 m Wrong or missing GEOID model; vgridshift step absent from the pipeline Install the correct g2018 grid and re-inspect the pipeline definition (Step 2).
Systematic height bias across the tile Manual two-transform chain sampled the geoid at the pre-shift horizontal position Let TransformerGroup build one ordered compound operation (Step 2).
NaN z after transform Point outside the GEOID grid’s coverage area Clip to the vertical grid extent or quarantine the feature (Step 3).
Combined budget passes but flood model rejects heights Gated only the combined figure, masking a vertical component failure Gate horizontal and vertical tolerances independently (Step 4).