Vertical Datum & Geoid Model Handling Jump to heading#

Horizontal coordinates arrive with a CRS attached often enough that a pipeline can be built around the expectation. Heights arrive as a column called elev. That asymmetry is the whole problem: an elevation of 214.7 is meaningless until you know whether it is measured from an ellipsoid or from a tidal datum, which realization of that datum, in what units, and via which geoid model — and the difference between those answers is routinely tens of metres. A flood model built on mismatched vertical references does not fail; it produces plausible, wrong answers. This stage makes height a first-class datum-bearing quantity inside CRS Normalization & Sync.

The scope boundary is narrow and deliberate. Datum transformation fallback chains owns the horizontal transformation graph and what to do when a grid is missing; this page owns the vertical component — declaring it, converting it, and stating its accuracy. Unit conversion and tolerance thresholds owns the foot-versus-metre question generally; here it recurs with a vertical twist, because a US survey foot elevation mislabelled as international feet drifts by about 0.6 mm per kilometre of height — negligible for a building, material for an aircraft procedure.

What a Height Actually Needs to Be Interpretable Jump to heading#

Four facts, and a dataset that carries fewer than four is carrying an opinion rather than a measurement.

The four facts a height column must carry are: the vertical datum (EVRF2019, NAVD88, a local tide gauge), the height type (ellipsoidal or orthometric — or, less commonly, normal or dynamic), the unit, and, where a conversion has occurred, the geoid model that produced it. Only the third is usually recorded. The convention this stage enforces is that a height column has a compound CRS or it has nothing: EPSG:25832 + EPSG:5730 says ETRS89 / UTM 32N horizontally and EVRF2007 vertically, and it is a statement a machine can act on.

Declarative Configuration Manifest Jump to heading#

# vertical_policy.yaml — pyproj >=3.6, PROJ >=9.3
policy_version: "1.3.0"            # MANDATORY
target_compound_crs: "EPSG:25832+EPSG:9389"   # MANDATORY: ETRS89/UTM32N + EVRF2019
geoid_models:                      # MANDATORY when converting ellipsoidal to orthometric
  - id: "eur_nkg_nkg2015.tif"      # the PROJ grid filename, pinned
    covers: "EPSG:9389"            # the vertical CRS it realizes
    accuracy_m: 0.02               # MANDATORY: published model accuracy
sources:
  - name: "lidar_dtm_2024"
    height_column: "z"             # MANDATORY
    declared_vertical_crs: "EPSG:4937"   # ETRS89 ellipsoidal heights
    height_type: ellipsoidal       # MANDATORY: ellipsoidal | orthometric
    unit: metre                     # MANDATORY
    source_accuracy_m: 0.15
  - name: "municipal_benchmarks"
    height_column: "elev_m"
    declared_vertical_crs: "EPSG:5730"   # EVRF2007 — an older realization
    height_type: orthometric
    unit: metre
    source_accuracy_m: 0.01
guards:
  reject_undeclared: true          # MANDATORY: a height column with no vertical CRS is refused
  plausibility_range_m: [-450, 9000]   # OPTIONAL: rejects sentinels and unit errors
  max_conversion_shift_m: 80.0     # OPTIONAL: a shift beyond this means the wrong model

reject_undeclared: true is the setting that does the work, and it is the one teams disable first when a delivery is late. The alternative — assuming the height is orthometric because it usually is — produces a dataset that is right most of the time and silently 40 m out in the places where it matters.

Preprocessing Requirements Jump to heading#

Horizontal position is settled first. A geoid separation is looked up at a location, so the horizontal coordinates must already be in a known CRS before any vertical conversion runs. Reprojecting afterwards is fine; converting height first is not.

The geoid grid is installed and pinned at build time. Geoid models are PROJ grids exactly like horizontal shift grids, and the same reproducibility argument applies as in handling missing NTv2 grid files in production: fetch at build time into a pinned PROJ_DATA, never at run time. A missing geoid grid does not raise — PROJ falls back to a coarser path or to no conversion at all, and the pipeline produces heights that are wrong by the separation.

Sentinels are removed before plausibility testing. -9999, -32768 and 0.0 are all common “no data” markers in elevation columns. Convert them to nulls at the reader, and count them: a source whose no-data rate jumps between deliveries has changed something upstream.

Execution Engine & Precision Guards Jump to heading#

# crs/vertical.py — pyproj >=3.6, PROJ >=9.3 — Python 3.10+
import logging
from dataclasses import dataclass

from pyproj import CRS, Transformer
from pyproj.transformer import TransformerGroup

logger = logging.getLogger("crs.vertical")


@dataclass(frozen=True)
class VerticalResult:
    height: float | None
    separation: float | None      # N applied, metres; None when no conversion was needed
    accuracy_m: float
    code: str                     # "" on success


def build_transformer(source_crs: str, target_crs: str) -> Transformer:
    """Pick the transformation explicitly and prove a geoid path is available."""
    group = TransformerGroup(CRS(source_crs), CRS(target_crs))
    if not group.available_operations:
        raise RuntimeError(f"no available vertical path {source_crs} -> {target_crs}; "
                           f"unavailable: {[o.name for o in group.unavailable_operations]}")
    chosen = group.available_operations[0]
    logger.info("vertical path %s -> %s via %s (accuracy %.3f m)",
                source_crs, target_crs, chosen.name, chosen.accuracy or -1.0)
    return Transformer.from_pipeline(chosen.to_proj4())


def convert_height(x: float, y: float, z: float, transformer: Transformer,
                   spec: dict, model_accuracy: float) -> VerticalResult:
    low, high = spec["plausibility_range_m"]
    if z is None:
        return VerticalResult(None, None, 0.0, "VERT_NULL_HEIGHT")
    if not (low <= z <= high):
        return VerticalResult(None, None, 0.0, "VERT_IMPLAUSIBLE_INPUT")

    _x, _y, z_out = transformer.transform(x, y, z)
    if z_out is None or z_out != z_out:                     # NaN check
        return VerticalResult(None, None, 0.0, "VERT_TRANSFORM_FAILED")

    separation = z - z_out
    if abs(separation) > spec["max_conversion_shift_m"]:
        return VerticalResult(None, separation, 0.0, "VERT_SHIFT_IMPLAUSIBLE")

    total_accuracy = (spec["source_accuracy_m"] ** 2 + model_accuracy ** 2) ** 0.5
    return VerticalResult(z_out, separation, total_accuracy, "")

Three guards deserve naming. The transformer is built from an explicitly chosen operation, not from Transformer.from_crs, because the convenience constructor silently selects a path and a missing geoid grid changes that selection without an error. The shift is checked against a plausibility ceiling: geoid separations are large but bounded — roughly −105 m to +85 m worldwide — so a computed shift of 300 m means the wrong grid or a horizontal position outside the model’s coverage. Accuracy is combined in quadrature, so the published vertical accuracy reflects the geoid model as well as the survey, which is the number a downstream flood model should be using.

Failure Modes & Fallback Routing Jump to heading#

The last row is the failure that survives review most often, because the two realizations look like the same datum in every user-facing label. Recording the exact EPSG code rather than a name is what prevents it, and it is why the manifest names EPSG:9389 rather than “EVRF”.

Merging Sources That Sit on Different Vertical References Jump to heading#

The single-source case is arithmetic. The multi-source case is where vertical work actually goes wrong, because the mismatch does not announce itself: two layers overlay perfectly in plan and disagree by a metre in height, and every derived product quietly inherits the disagreement.

The mismatch is invisible in every horizontal check. A lidar surface on ellipsoidal heights and a benchmark network on a national datum will pass every CRS gate you have, because the horizontal components agree. Nothing in a bounding box, a projection test or a topology rule looks at the third dimension. The only detector is a declared vertical CRS per source and an assertion that they match before the merge — which is exactly what reject_undeclared buys.

Ordering matters when both components change. A dataset moving from an old horizontal datum and an old vertical datum must have the horizontal shift applied first, because the geoid separation is looked up at a position, and the position is about to move. On a national scale that is a metre or two of horizontal movement, which changes the separation by a millimetre or so — negligible for most work and not for geodetic control. Where both shifts are needed, use one PROJ pipeline that carries both steps rather than two passes, so the intermediate state never exists as data. The chaining of vertical and horizontal datum shifts covers the pipeline construction in detail.

Accuracy does not survive a merge unchanged. When a 2 cm benchmark network and a 15 cm lidar surface are combined, the result is not a 2 cm product. Carry accuracy per feature rather than per dataset, combine in quadrature at every conversion, and publish the distribution rather than a single figure — a merged layer whose vertical accuracy ranges from 2 cm to 15 cm is honestly described by that range and dishonestly described by either end of it.

Local datums are more common than anyone expects. Utility asset records, mine surveys, airport datasets and long-lived industrial sites frequently carry heights on a datum established by a private survey decades ago, related to the national network by a single published offset — or by none at all. These cannot be converted; they can only be labelled. Give a local datum its own identifier in the manifest, refuse to merge it into a national-datum layer, and record the relationship as an offset with an accuracy if one exists. The failure mode being avoided is a dataset in which most heights are national and a few hundred are local, with nothing distinguishing them.

A short pre-merge assertion catches all four cases at once:

# crs/vertical_merge.py — Python 3.10+
def assert_vertically_compatible(sources: list[dict]) -> None:
    """Refuse a merge whose inputs do not share one declared vertical reference."""
    declared = {s["name"]: s.get("declared_vertical_crs") for s in sources}
    missing = [name for name, crs in declared.items() if not crs]
    if missing:
        raise ValueError(f"sources with no declared vertical CRS: {missing}")

    distinct = set(declared.values())
    if len(distinct) > 1:
        raise ValueError(
            f"sources sit on different vertical references {sorted(distinct)}; "
            "convert every source to the target compound CRS before merging"
        )

    types = {s["height_type"] for s in sources}
    if len(types) > 1:
        raise ValueError(f"mixed height types {sorted(types)} — convert before merging")

Run it before the merge rather than after, and run it on the declared metadata rather than on the values: a height column whose numbers look plausible is not evidence of anything, and by the time a value looks wrong the merge has already happened.

Compliance Reporting Output Jump to heading#

# crs/vertical_report.py — pyarrow >=14 — Python 3.10+
import pyarrow as pa

VERTICAL_REPORT_SCHEMA = pa.schema([
    ("feature_id",        pa.string()),
    ("run_id",            pa.string()),
    ("policy_version",    pa.string()),
    ("source_vertical_crs", pa.string()),   # exact EPSG code, never a name
    ("target_vertical_crs", pa.string()),
    ("height_type_in",    pa.string()),     # ellipsoidal | orthometric
    ("geoid_model",       pa.string()),     # the pinned grid filename actually used
    ("separation_m",      pa.float64()),    # N applied at this location
    ("accuracy_m",        pa.float64()),    # combined in quadrature
    ("code",              pa.string()),
])

Storing the separation per feature rather than per run looks extravagant until the first time someone asks whether a particular building’s height was converted with the 2015 or the 2020 model. It is one float per feature, and it makes the conversion reproducible without re-running anything. The accuracy_m column is the one that must reach the published metadata: an ISO 19115 lineage statement that reports horizontal accuracy and omits vertical accuracy is describing half the dataset.

CI Integration Jump to heading#

# tests/test_vertical.py — pytest >=7, pyproj >=3.6
def test_geoid_grid_is_installed_and_pinned():
    group = TransformerGroup(CRS("EPSG:4937"), CRS("EPSG:9389"))
    assert group.available_operations, \
        f"missing geoid grids: {[o.grids for o in group.unavailable_operations]}"


def test_known_benchmark_converts_within_a_centimetre():
    # A published benchmark with both h and H recorded by the national agency.
    result = convert_height(x=500_000.0, y=5_600_000.0, z=245.812,
                            transformer=TRANSFORMER, spec=SPEC, model_accuracy=0.02)
    assert abs(result.height - 200.417) < 0.01
    assert abs(result.separation - 45.395) < 0.01


def test_feet_in_a_metre_column_is_rejected():
    # Negative control: 806 "metres" is 245.8 m in feet — plausible-looking, wrong.
    result = convert_height(500_000.0, 5_600_000.0, 806.5, TRANSFORMER, SPEC, 0.02)
    assert result.code in {"VERT_IMPLAUSIBLE_INPUT", "VERT_SHIFT_IMPLAUSIBLE"} or \
           abs(result.height - 761.1) < 1.0     # documented: the value converts, but wrongly

The benchmark test is the one worth the effort of finding real published values for. It is the only assertion in the suite that catches a correct-looking misconfiguration — a geoid model that is installed, applied without error, and simply not the one the target datum requires.

Deeper Implementation Walkthroughs Jump to heading#

Converting ellipsoidal to orthometric heights with geoid grids implements the conversion end to end, including grid installation, coverage checks and the benchmark verification. Planning the NAVD88 to NAPGD2022 transition treats the datum change as a migration project rather than a code change — dual-publication, impact assessment, and how to state which datum a delivery is on.

Frequently Asked Questions Jump to heading#

Why can a GPS receiver report a height that differs from the sign on the mountain? Because they measure from different surfaces. A raw GNSS position is ellipsoidal; the sign quotes an orthometric height above a national vertical datum. The difference is the geoid separation at that location, which ranges roughly from −105 m to +85 m worldwide and changes by tens of metres across a single country. Neither number is wrong; a dataset that mixes them without labelling is.

Is a compound EPSG code always available for our combination? Not always, and where it is not, the correct answer is to declare the horizontal and vertical codes as a pair rather than to pick an approximate compound code. PROJ handles a pipeline built from two explicit codes perfectly well, and the pair records exactly what was used. Picking “the nearest compound code” is how a realization mismatch enters a dataset with full documentation and no error.

Can we skip the geoid and just apply a constant offset? Over a small area with a locally-fitted constant and a stated accuracy, this is a legitimate engineering approximation, and it must be recorded as one — the constant, the area it applies to, and the residual against benchmarks. What it must never be is an undocumented shortcut, because the resulting heights look exactly like properly converted ones and cannot be corrected later without knowing the constant.

How do we handle sources that predate the current realization? Convert them into the target realization at ingest, record both the source and target codes per feature, and keep the original value. Retro-converting a published dataset is far harder than carrying the original, and the per-feature record is what makes a future re-conversion possible when the next realization arrives.

Related Jump to heading#

• Coordinate Reference System Normalization & Sync — the parent section and the horizontal pipeline this stage runs beside
• Converting Ellipsoidal to Orthometric Heights with Geoid Grids — the conversion procedure with benchmark verification
• Planning the NAVD88 to NAPGD2022 Transition — running two vertical datums during a national change
• Chaining Vertical and Horizontal Datum Shifts — ordering the two components in one pipeline
• Handling Missing NTv2 Grid Files in Production — the same grid-pinning discipline applied horizontally
• Detecting US Survey Foot vs International Foot Errors — the unit failure that also appears in height columns