Converting Ellipsoidal to Orthometric Heights in Python

To convert an ellipsoidal height to an orthometric one in Python, build both sides of the transformation as compound coordinate reference systems and pass z to Transformer.transform. PROJ then selects a vertical transformation, looks up the geoid separation for that horizontal position, and returns the orthometric height. Omit either the compound definition or the z argument and the height comes back exactly as it went in, with no error raised. This page is part of the Vertical Datums and Height Systems reference.

How pyproj Decides to Transform a Height

A Transformer is built from two coordinate reference systems, and it can only perform operations those systems describe. Two horizontal definitions describe no vertical relationship, so a transformer built from them has no vertical operation to apply. It is not ignoring the height — it was never given a reason to change it.

A compound CRS supplies that reason. It pairs a horizontal CRS with a vertical one, and PROJ’s operation search then includes vertical transformations: geoid grid interpolations, datum offsets, or a combination. The selected operation is what applies the separation.

The other half of the mechanism is the grid itself. A national geoid model is a raster of separation values, distributed separately from PROJ. When it is present the interpolation is centimetre-level; when it is absent PROJ may fall back to a coarse global model and still return an answer. This is the single most important thing to verify, because the fallback is silent and the difference is metres.

Production-Ready Script

# pyproj>=3.5.0 (PROJ 9.x), numpy>=1.24, Python 3.9+
from __future__ import annotations

import numpy as np
from pyproj import CRS, Transformer


class VerticalTransformError(RuntimeError):
    """Raised when the transformation is not actually three-dimensional."""


def build_vertical_transformer(
    src_epsg: int,
    dst_horizontal_epsg: int,
    dst_vertical_epsg: int,
    *,
    require_grid: bool = True,
) -> Transformer:
    """Transformer from a 3D geographic CRS to a compound projected + vertical CRS.

    src_epsg must be a 3D CRS (ellipsoidal height), e.g. 4937 for ETRS89.
    """
    src = CRS.from_epsg(src_epsg)
    dst = CRS.from_string(f"EPSG:{dst_horizontal_epsg}+{dst_vertical_epsg}")

    if len(src.axis_info) < 3:
        raise VerticalTransformError(
            f"EPSG:{src_epsg} is 2D — use the 3D realisation, or z will pass through"
        )
    if not dst.is_compound:
        raise VerticalTransformError("destination is not a compound CRS")

    transformer = Transformer.from_crs(src, dst, always_xy=True)

    if require_grid:
        op = transformer.get_last_used_operation() if hasattr(
            transformer, "get_last_used_operation") else None
        # The operation is only resolved after a first transform on some PROJ builds,
        # so probe with a representative coordinate inside the grid's extent.
        transformer.transform(0.0, 51.5, 0.0)
        _assert_grids_available(transformer)

    return transformer


def _assert_grids_available(transformer: Transformer) -> None:
    op = transformer.get_last_used_operation()
    missing = [g.short_name for g in op.grids if not g.available]
    if missing:
        raise VerticalTransformError(
            f"geoid grid(s) not installed: {', '.join(missing)} — "
            "PROJ would fall back to a coarse global model"
        )


def to_orthometric(
    transformer: Transformer,
    lon: np.ndarray,
    lat: np.ndarray,
    h_ellipsoidal: np.ndarray,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Ellipsoidal height -> orthometric height, vectorised over whole arrays."""
    east, north, H = transformer.transform(
        np.asarray(lon, dtype=float),
        np.asarray(lat, dtype=float),
        np.asarray(h_ellipsoidal, dtype=float),
    )
    separation = np.asarray(h_ellipsoidal) - np.asarray(H)
    if np.allclose(separation, 0.0):
        raise VerticalTransformError(
            "separation is zero everywhere — the transform is horizontal"
        )
    return np.asarray(east), np.asarray(north), np.asarray(H)


if __name__ == "__main__":
    t = build_vertical_transformer(4937, 27700, 5701)   # ETRS89 3D -> BNG + ODN
    lon = np.array([-1.54785, -1.54600])
    lat = np.array([53.80139, 53.80200])
    h = np.array([96.412, 97.001])
    e, n, H = to_orthometric(t, lon, lat, h)
    for i in range(len(H)):
        print(f"h={h[i]:.3f}  ->  H={H[i]:.3f}  (N={h[i] - H[i]:.3f} m)")

Key implementation notes:

• build_vertical_transformer refuses a 2D source rather than transforming it. That refusal is the whole value of the function: a 2D source is the failure this page exists to prevent.
• The zero-separation assertion in to_orthometric is a second, independent guard. It catches the case where the definitions look right but the operation PROJ chose does nothing vertical.
• always_xy=True is set for the same reason as in any horizontal transformation — it fixes the argument order to longitude, latitude regardless of the authority axis order.
• Whole arrays are passed to transform. The grid interpolation is vectorised internally; per-point calls spend their time in the Python loop.

Compatibility Matrix

Fallback Strategies

1. The source is already orthometric. Applying this conversion produces an error of exactly 2N. There is no way to detect it from the numbers alone at a single point, but a dataset whose heights differ from local mapped ground levels by roughly the known regional separation is the signature. Record the height system per source, as described on the parent page, rather than inferring it.

2. No national grid installed. _assert_grids_available raises. Install the grid package at image build time and set PROJ_NETWORK=OFF so a missing grid is an error rather than a runtime download that may or may not succeed.

3. Points outside the grid extent. A national grid covers a bounded area. Points outside it transform through whatever fallback PROJ finds, which for a site near a border may silently switch model partway through a dataset. Assert the input bounding box against the grid extent before transforming.

4. Mixed vertical datums in one dataset. A merged deliverable can contain sections on different vertical datums. There is no way for one transformer to handle this correctly. Split the dataset by source, transform each with its own transformer, and merge afterwards.

5. Heights that are already project-datum values. A CAD or BIM elevation is not an ellipsoidal height and this conversion does not apply to it. Resolve the project offset first — see Reconciling BIM Project Elevation with a National Datum — and only then treat the values as national-datum heights.

FAQ

Related Pages

• Vertical Datums and Height Systems — parent reference on the three surfaces a height can be measured from
• Applying a Geoid Model with pyproj — sibling guide on installing, pinning and verifying the grid this conversion depends on
• Reprojecting CAD Coordinates with pyproj Transformer — the horizontal transformer whose caching rules apply here too