Vertical Datums and Height Systems

A vertical datum is the surface a Z value is measured from, and reconciling the several that arrive in one project is the part of the Coordinate Transformation & Spatial Alignment pipeline that most often survives review unnoticed and fails on site.

Horizontal misalignment is visible: an asset drawn in the wrong field is obvious on a map. A vertical error is not. A pipe invert 2.4 m below where the model says it is looks entirely reasonable in plan, in section, and in every automated check that does not know what surface the height was measured from. The datasets that meet in an infrastructure project routinely carry three different height systems at once — GNSS survey in ellipsoidal height, levelled data on a national orthometric datum, and a BIM model on a project datum whose origin was chosen for drafting convenience — and none of the three formats involved compels anyone to say which is which.

Prerequisites

• Python 3.9+ with type hints and pathlib.
• pyproj>=3.5.0, built against PROJ 9.x. Earlier PROJ releases handle compound coordinate reference systems less consistently, and the grid-fetching behaviour differs.
• Network access or a pinned grid package for geoid models, if survey-grade results are required.
• Working knowledge of your project’s horizontal CRS and, critically, of which vertical datum each incoming dataset was observed on — this is a question for the surveyor, not something to infer.

# pyproj>=3.5.0 (PROJ 9.x)
pip install "pyproj>=3.5.0"

# Fetch the geoid grids for a bounding box, once, at image build time
projsync --bbox -8,49,2,61 --file uk_os_OSGM15_GB.tif

Verify the installation resolves a vertical transformation before writing pipeline code:

# pyproj>=3.5.0
from pyproj import CRS, Transformer

compound = CRS.from_epsg(7405)      # OSGB36 / British National Grid + ODN height
print(compound.is_compound)          # True
print([sub.name for sub in compound.sub_crs_list])

If is_compound is False, the code you write against it will move X and Y and leave Z alone.

Architectural Overview

Three surfaces are involved, and every height in a project is measured from one of them.

The reference ellipsoid is a smooth mathematical figure fitted to the Earth. It is what a horizontal CRS is defined against, and a GNSS receiver reports position relative to it. Heights measured from it are ellipsoidal heights, conventionally written h.

The geoid is the equipotential surface of the Earth’s gravity field that best approximates mean sea level. It is not smooth — it undulates with the distribution of mass beneath it — and it is what water responds to. Heights measured from it are orthometric heights, written H, and every drainage gradient, flood level and building-height limit is defined against a national realisation of it.

The geoid separation N is the difference between the two at a given horizontal position, so h = H + N. N is not a constant. Over a country it varies by tens of metres and over a large site it can vary by tens of centimetres, which is why applying a single average separation across a project is a defensible approximation for planning and an indefensible one for setting out.

On top of these sits the project datum used inside CAD and BIM: an arbitrary origin, often chosen so that a floor level is a round number. Its relationship to a national datum is a single constant offset, and that constant lives in the project setup rather than in the geometry.

The mechanism that matters in code is the compound CRS: a coordinate reference system that pairs a horizontal CRS with an explicit vertical CRS. PROJ selects a vertical transformation only when both sides of a transformation declare one. Given two plain 2D definitions it has nothing to work with, so it moves X and Y and returns Z unchanged — correctly, quietly, and usually wrongly.

Step-by-Step Implementation

1. Establish what each source is actually measured from

Before any code, record for every input: horizontal CRS, vertical datum, and whether the Z values have already had a geoid correction applied. A GNSS deliverable that has been “corrected to ODN” by the surveyor’s software is already orthometric, and applying the correction again shifts it by twice the separation.

# Python 3.9+
from dataclasses import dataclass

@dataclass(frozen=True)
class HeightSource:
    name: str
    horizontal_epsg: int
    vertical: str          # "ellipsoidal" | "orthometric" | "project"
    vertical_epsg: int | None = None   # required when vertical == "orthometric"
    project_offset_m: float | None = None  # required when vertical == "project"

    def __post_init__(self):
        if self.vertical == "orthometric" and self.vertical_epsg is None:
            raise ValueError(f"{self.name}: orthometric height needs a named vertical datum")
        if self.vertical == "project" and self.project_offset_m is None:
            raise ValueError(f"{self.name}: project height needs a documented offset")

The class refuses to describe a source ambiguously. That refusal is the point: an unlabelled height is the failure this whole page exists to prevent, and it is much cheaper to be stopped here than to be corrected on site.

2. Build compound coordinate reference systems

# pyproj>=3.5.0
from pyproj import CRS

def compound(horizontal_epsg: int, vertical_epsg: int) -> CRS:
    """A horizontal CRS paired with an explicit vertical CRS."""
    crs = CRS.from_string(f"EPSG:{horizontal_epsg}+{vertical_epsg}")
    assert crs.is_compound, "expected a compound CRS"
    return crs

# ETRS89 geographic 3D (ellipsoidal height) -> British National Grid + ODN
SRC = CRS.from_epsg(4937)                 # ETRS89, ellipsoidal height
DST = compound(27700, 5701)               # OSGB36 / BNG + ODN orthometric height

3. Transform in three dimensions

# pyproj>=3.5.0, numpy>=1.24
import numpy as np
from pyproj import Transformer

transformer = Transformer.from_crs(SRC, DST, always_xy=True)

def to_orthometric(lon, lat, h_ellipsoidal):
    """Returns easting, northing and ORTHOMETRIC height."""
    e, n, H = transformer.transform(lon, lat, h_ellipsoidal)
    return np.asarray(e), np.asarray(n), np.asarray(H)

Two details decide whether this is correct. The z argument must actually be passed — omit it and the transformer has no height to operate on. And the separation applied is the one PROJ found for that horizontal position, which means the transformation is only as good as the grid installed; see the validation section for how to confirm which grid was used.

4. Fold in the project datum offset

A CAD or BIM project height is related to a national datum by one constant, established at project setup. Apply it as a named step rather than folding it into the CRS, because it is a project fact and not a geodetic one:

# Python 3.9+
def project_to_national(z_project: "np.ndarray", offset_m: float) -> "np.ndarray":
    """Project elevation -> national orthometric datum.

    offset_m is the height of the project's zero level above the national datum,
    taken from the project setup document — never inferred from a single element.
    """
    return z_project + offset_m

5. Verify against a levelled benchmark

The transformation is only verified by a point whose height is independently known. Transform the benchmark’s coordinates and compare.

# pyproj>=3.5.0
def verify_benchmark(lon, lat, h_ellipsoidal, published_H, tol=0.02):
    _, _, H = transformer.transform(lon, lat, h_ellipsoidal)
    residual = abs(H - published_H)
    if residual > tol:
        raise AssertionError(
            f"vertical residual {residual:.3f} m exceeds {tol} m — "
            "check the geoid grid and the source height system"
        )
    return residual

Edge Cases and Gotchas

A 2D transformer silently passes height through. This is the dominant failure and it produces output that is entirely plausible: X and Y are correct, Z is unchanged, and no exception is raised. Assert crs.is_compound on both sides of every transformation that carries height, and treat a transformer built from two 2D definitions as a bug when z is in play.

Double-correcting an already-corrected survey. GNSS post-processing software frequently outputs orthometric height directly. Applying the geoid separation to that output shifts it by 2N, which in the UK is roughly 100 m and therefore obvious, and in a region where N is near zero is roughly nothing and therefore invisible. Record per source whether the correction has been applied, as in step 1.

Sign conventions on the separation. h = H + N is the standard relation, but published separation values are occasionally tabulated as H - h. A sign error is exactly 2N, which looks like the double-correction above. Confirm against a benchmark rather than against the documentation.

Missing grids fall back silently. Where a national grid is not installed, PROJ may fall back to a coarser global model and still return a result. The result is a metre-level approximation presented with the same confidence as a centimetre-level one. Check which operation was actually used:

# pyproj>=3.5.0
print(transformer.description)
for grid in transformer.get_last_used_operation().grids:
    print(grid.name, "available:", grid.available)

Tunnel and bridge deliverables mixing systems within one model. Long linear assets can span enough of a country for N to vary measurably along their length. A single applied separation is then wrong at both ends in opposite directions. Transform per coordinate rather than applying a project-wide constant.

IFC OrthogonalHeight is an offset, not a height. It records where the model origin sits relative to the projected coordinate reference system, so it is added to model Z rather than replacing it. Reading it as an absolute height puts the whole model at that elevation.

Validation and Testing

Vertical correctness needs two tests: one that the transformation is three-dimensional at all, and one that it lands where an independent measurement says it should.

# pytest, pyproj>=3.5.0
import pytest
from pyproj import CRS, Transformer

def test_transformer_is_three_dimensional():
    src = CRS.from_epsg(4937)                       # ETRS89 3D
    dst = CRS.from_string("EPSG:27700+5701")        # BNG + ODN
    assert src.is_geographic and dst.is_compound
    t = Transformer.from_crs(src, dst, always_xy=True)
    _, _, z_out = t.transform(-1.5, 52.0, 100.0)
    # A 2D pipeline would return the input height unchanged.
    assert abs(z_out - 100.0) > 1.0, "height passed through — the transform is horizontal"

@pytest.mark.parametrize("lon,lat,h,published_H", [
    (-1.54785, 53.80139, 96.412, 47.31),   # a levelled benchmark, ellipsoidal h -> ODN H
])
def test_benchmark_within_tolerance(lon, lat, h, published_H):
    t = Transformer.from_crs(CRS.from_epsg(4937),
                             CRS.from_string("EPSG:27700+5701"), always_xy=True)
    _, _, H = t.transform(lon, lat, h)
    assert abs(H - published_H) < 0.05

The first test is the more valuable of the two, because it fails loudly on the mistake that otherwise never announces itself. Substitute benchmarks and datums for your own region; the shape of both tests is unchanged.

Beyond unit tests, assert plausibility on every run. Ground level in a project has a known range; a transformed dataset whose heights sit 45 m from where the site is known to be has usually had a separation applied twice or not at all, and a single range assertion catches both.

Performance and Scale

Vertical transformation costs essentially nothing beyond the horizontal one — the geoid lookup is an interpolation in a raster the transformer has already opened — so the performance considerations are about the grids rather than the arithmetic.

Fetch grids at image build time, never at run time. A container that calls out to the PROJ CDN on first use has a cold-start dependency on a network service, and a worker that cannot reach it falls back to a coarser model rather than failing, which converts an outage into a silent accuracy regression. Pin the grid package, install it in the image, and set PROJ_NETWORK=OFF so a missing grid is an error rather than a download.

Build the transformer once and reuse it, as with any PROJ pipeline; the reasoning is identical to that in Reprojecting CAD Coordinates with pyproj Transformer and the saving is larger here because the grid is opened during construction. For large arrays, pass whole numpy arrays to transform() rather than iterating; the interpolation is vectorised internally and per-point calls spend their time in Python.

FAQ

Related Pages

• Coordinate Transformation & Spatial Alignment — the parent pipeline this height work is a stage of
• CRS Normalization Workflows — the horizontal counterpart — detection, validation and one cached transformation
• Unit Conversion Pipelines — height is a length, and inherits every unit trap on that page
• Reading IFC Georeferencing with ifcopenshell — where a BIM model records its orthogonal height offset
• Scale and Rotation Synchronization — the horizontal alignment that a vertical shift is often mistaken for