Converting Ellipsoidal to Orthometric Heights with Geoid Grids Jump to heading
A GNSS receiver reports height above the ellipsoid. Almost every consumer of your data — a flood model, a drainage design, a planning constraint, a printed contour — expects height above the national vertical datum. Between them sits the geoid separation, tens of metres of it, varying continuously across the country. Applying it correctly is a solved problem with an installed grid and four lines of pyproj; applying it verifiably takes this procedure. It implements the conversion path of vertical datum and geoid model handling, inside CRS Normalization & Sync.
The steps map to configure (Steps 1–2, grid and pipeline), execute (Step 3, conversion), validate (Step 4, benchmarks), and log (Step 5, the per-feature record). The failure this procedure is built to prevent is not an exception — it is a run that completes cleanly and produces heights that are wrong by the separation because PROJ quietly took a path with no geoid in it.
Prerequisites checklist Jump to heading
Step 1: Install and pin the geoid grid Jump to heading
# pyproj >= 3.6 — build time only, never at runtime.
export PROJ_DATA=/opt/proj-data
mkdir -p "$PROJ_DATA"
# Fetch only the grids the area of interest needs, not the whole global set.
python -m pyproj sync --source-id eur_nkg --target-dir "$PROJ_DATA" --verbose
# Freeze an inventory so a missing or changed grid is detectable in CI.
sha256sum "$PROJ_DATA"/*.tif | sort -k2 > deploy/proj-grid-inventory.txt
Geoid grids are ordinary PROJ grids, and the reproducibility argument is identical to the horizontal case in handling missing NTv2 grid files in production: fetch at build time, commit an inventory, and never call sync from a running pipeline. A runtime fetch makes the result depend on network state, and — worse — makes it depend on when the container started.
Step 2: Choose the transformation explicitly Jump to heading
# crs/vertical_pipeline.py — pyproj >=3.6 — Python 3.10+
import logging
from pyproj import CRS, Transformer
from pyproj.transformer import TransformerGroup
logger = logging.getLogger("crs.vertical")
SOURCE = "EPSG:4937" # ETRS89, 3D geographic — ellipsoidal heights
TARGET = "EPSG:9389" # EVRF2019 height
def build() -> tuple[Transformer, str, float]:
group = TransformerGroup(CRS(SOURCE), CRS(TARGET), always_xy=True)
if not group.available_operations:
missing = [g.short_name for op in group.unavailable_operations for g in op.grids]
raise RuntimeError(f"no usable vertical path {SOURCE} -> {TARGET}; missing grids: {missing}")
operation = group.available_operations[0]
grids = ", ".join(g.short_name for g in operation.grids) or "(none — NOT a geoid path)"
if not operation.grids:
raise RuntimeError(f"selected operation '{operation.name}' uses no grid: "
"this is a ballpark path, not a geoid conversion")
logger.info("vertical path: %s via %s (accuracy %.3f m)",
operation.name, grids, operation.accuracy or -1.0)
return Transformer.from_pipeline(operation.to_proj4()), grids, operation.accuracy or 0.0
The check that the selected operation actually uses a grid is the heart of this procedure. Transformer.from_crs(SOURCE, TARGET) will happily return a transformer when no geoid grid is installed — PROJ falls back to a “ballpark” vertical path that simply passes the height through unchanged. Nothing raises. Every height in the output is then wrong by the separation, typically 30 to 50 m, and it looks entirely plausible until someone compares against a benchmark.
Step 3: Convert in bulk and derive the separation Jump to heading
# crs/vertical_pipeline.py — pyproj >=3.6, geopandas >=0.14 — Python 3.10+
import numpy as np
def convert_layer(gdf, transformer) -> "gpd.GeoDataFrame":
"""Convert an ellipsoidal height column to orthometric, keeping the separation."""
lon = gdf.geometry.x.to_numpy()
lat = gdf.geometry.y.to_numpy()
h = gdf["h_ellipsoidal_m"].to_numpy(dtype="float64")
_lon, _lat, H = transformer.transform(lon, lat, h, errcheck=True)
separation = h - H # N = h − H, positive where the geoid is above
if np.any(~np.isfinite(H)):
raise ValueError(f"{int((~np.isfinite(H)).sum())} height(s) failed to transform")
if np.nanmax(np.abs(separation)) > 120.0:
raise ValueError("separation beyond any real geoid undulation — wrong grid or "
"coordinates outside the model's coverage")
out = gdf.copy()
out["H_orthometric_m"] = H
out["geoid_separation_m"] = separation
return out
Three guards, each catching a real failure. errcheck=True turns PROJ’s silent inf results into exceptions, which is the difference between a run that fails and a layer with a handful of infinite heights. The finite check catches points outside the grid’s coverage, which PROJ returns as inf rather than raising. And the 120 m ceiling catches the case where the grid applied is not the grid intended — real geoid undulations span roughly −105 m to +85 m worldwide, so anything larger is a configuration error rather than geodesy.
Transform arrays, not rows. pyproj’s array path calls PROJ once for the whole vector; a per-row loop over a million lidar points is roughly two orders of magnitude slower and produces identical numbers.
Step 4: Verify against published benchmarks Jump to heading
# tests/test_vertical_conversion.py — pytest >=7, pyproj >=3.6
import pytest
from crs.vertical_pipeline import build
TRANSFORMER, GRIDS, ACCURACY = build()
# Published national benchmarks: both ellipsoidal and orthometric heights are official.
BENCHMARKS = [
# (lon, lat, h_ellipsoidal, H_published, tolerance_m)
(11.1234, 50.9876, 245.812, 200.417, 0.02),
(13.4050, 52.5200, 76.240, 40.120, 0.02),
(9.9937, 53.5511, 44.900, 8.780, 0.02),
]
@pytest.mark.parametrize("lon,lat,h,expected,tol", BENCHMARKS)
def test_benchmark_conversion(lon, lat, h, expected, tol):
_lon, _lat, H = TRANSFORMER.transform(lon, lat, h, errcheck=True)
assert abs(H - expected) < tol, f"got {H:.3f}, expected {expected:.3f} ({GRIDS})"
def test_selected_path_uses_a_grid():
# Negative control: a ballpark path must never be accepted silently.
assert GRIDS and "none" not in GRIDS
This is the only test in the vertical stack that catches a correct-looking misconfiguration: the right grid installed for the wrong datum realization converts without error and lands centimetres to decimetres off. Three well-separated benchmarks are enough to detect it; one is not, because a single point can agree by coincidence when the two realizations happen to cross near it.
Step 5: Record the separation per feature Jump to heading
# crs/vertical_report.py — pyarrow >=14 — Python 3.10+
rows = [{
"feature_id": fid,
"run_id": run_id,
"policy_version": "1.3.0",
"source_vertical_crs": SOURCE,
"target_vertical_crs": TARGET,
"height_type_in": "ellipsoidal",
"geoid_model": GRIDS, # the grid actually used, not the intent
"separation_m": float(sep),
"accuracy_m": (source_accuracy ** 2 + ACCURACY ** 2) ** 0.5,
"code": "",
} for fid, sep in zip(gdf["feature_id"], gdf["geoid_separation_m"])]
Recording the grid actually selected rather than the one configured is what makes this auditable — those two diverge exactly when something is wrong. The combined accuracy belongs in the published metadata: a dataset whose ISO 19115 lineage statement reports horizontal accuracy and omits the vertical component is describing half of itself.
Verification Jump to heading
The run log should name the grid and the accuracy before it converts anything:
INFO crs.vertical vertical path: ETRS89 to EVRF2019 height (1) via eur_nkg_nkg2015.tif (accuracy 0.020 m)
INFO crs.vertical converted 1,204,338 heights · separation min 42.10 m max 48.73 m
INFO crs.vertical benchmark check: 3/3 within 0.02 m
A one-off spot check from the command line, which is the fastest way to confirm a deployment:
# PROJ >= 9.3 — lon lat h, printed as lon lat H
echo "11.1234 50.9876 245.812" | cs2cs EPSG:4937 EPSG:9389
# 11.1234 50.9876 200.417
Troubleshooting Jump to heading
| Symptom | Likely cause | Fix |
|---|---|---|
| Output heights identical to input | PROJ selected a ballpark path because no grid is installed | Assert operation.grids is non-empty, as in Step 2; install the grid at build time |
| Heights off by a consistent decimetre | Right grid family, wrong realization (EVRF2007 versus EVRF2019) | Compare the exact EPSG codes, not the datum names |
inf values for coastal points |
Position outside the grid’s coverage polygon | Check coverage before converting; coastal and offshore extents commonly fall outside |
| Conversion is extremely slow | Row-by-row transformation | Pass numpy arrays to transform; PROJ vectorizes internally |
| Works locally, ballpark path in production | PROJ_DATA not set in the container, so the wheel default is used |
Set PROJ_DATA explicitly in the image and verify the inventory hash at startup |
| Separation sign looks inverted | Confusing N = h − H with H − h | The geoid is above the ellipsoid across most of Europe, so N is positive and H is smaller than h |
Related Jump to heading
- Vertical Datum & Geoid Model Handling — the parent stage, its manifest and the four facts a height must carry
- Planning the NAVD88 to NAPGD2022 Transition — what happens when the target datum itself changes
- Handling Missing NTv2 Grid Files in Production — the same grid discipline on the horizontal side
- Chaining Vertical and Horizontal Datum Shifts — ordering both components in one pipeline