Applying a Geoid Grid with pyproj Transformers

The conversion from ellipsoidal to orthometric height is one function call, and it is a function call with four ways to be silently wrong: the wrong CRS pair, an axis-order surprise, a grid that was not installed, and a transformer rebuilt inside a loop so slowly that somebody replaces it with a constant.

This page covers doing it correctly and confirming it, as the implementation detail behind geoid models and vertical datum automation.

The four ways it goes wrong

Wrong source CRS. EPSG:4326 is two-dimensional and carries no height; a transformation from it cannot convert a height because there is no height to convert. EPSG:4979 is its three-dimensional counterpart and is what a GNSS position belongs to.

Axis order. EPSG definitions specify an axis order, and for geographic systems that order is latitude-then-longitude. Code written assuming longitude first gets coordinates swapped unless always_xy=True is passed, and the result for a site at 51° north, 1° west is a position in the Indian Ocean.

Missing grid. PROJ falls back to a lower-accuracy operation when the best one’s grid is absent, and completes without complaint. The difference is metres.

Rebuilt transformer. Constructing a transformer is expensive; using it is cheap. A transformer created inside a per-point loop makes a bulk conversion thousands of times slower, which is the usual reason somebody abandons it for a constant undulation.

Four failure modes of a geoid transformation and their symptoms Four rows pairing a mistake with its symptom. Using a two-dimensional source CRS leaves the height unchanged, so the output looks plausible and nothing was converted. Ignoring axis order swaps latitude and longitude, placing the survey thousands of kilometres away, which is at least obvious. A missing grid causes a silent fallback to a lower-accuracy transformation, wrong by metres. Rebuilding the transformer inside a loop is thousands of times slower, which usually leads somebody to replace it with a constant undulation and reintroduce a tilt. mistake symptom 2D source CRS (EPSG:4326) height unchanged, output plausible axis order assumed survey lands in another hemisphere grid not installed silent fallback, wrong by metres transformer rebuilt per point thousands of times slower Two of the four are obvious and two are not. The obvious ones cost an hour; the silent ones cost a delivered survey.

Figure 1 — Four mistakes, and which of them announce themselves.

Minimal reproducible solution

import numpy as np
from pyproj import CRS
from pyproj.transformer import Transformer, TransformerGroup


def build_height_transformer(horizontal_epsg: int, vertical_epsg: int) -> dict:
    """A transformer from 3D WGS84 to a projected compound CRS, checked.

    Built once and reused. The availability check is what distinguishes a
    grid-based conversion from a silent fallback, and it costs nothing after
    construction.
    """
    source = CRS.from_epsg(4979)                       # 3D geographic; 4326 has no height
    target = CRS.from_user_input(f"EPSG:{horizontal_epsg}+{vertical_epsg}")

    group = TransformerGroup(source, target, always_xy=True)
    if not group.best_available:
        missing = [g.short_name for op in group.unavailable_operations
                   for g in op.grids if not g.available]
        raise RuntimeError(
            f"the best transformation needs grids that are not installed: "
            f"{sorted(set(missing))} — results would use a lower-accuracy fallback")

    return {"transformer": Transformer.from_crs(source, target, always_xy=True),
            "description": group.transformers[0].description,
            "accuracy_m": group.transformers[0].target_crs and
                          getattr(group.transformers[0], "accuracy", None),
            "target": target.to_string()}


def convert_heights(transformer, lon: np.ndarray, lat: np.ndarray,
                    ellipsoidal_h: np.ndarray) -> dict:
    """Bulk conversion through a pre-built transformer."""
    x, y, z = transformer.transform(np.asarray(lon), np.asarray(lat),
                                    np.asarray(ellipsoidal_h))
    undulation = np.asarray(ellipsoidal_h) - np.asarray(z)
    return {"x": x, "y": y, "orthometric_h": z, "undulation_m": undulation,
            "undulation_mean_m": float(np.mean(undulation)),
            "undulation_range_m": float(np.ptp(undulation))}

Raising rather than warning on an unavailable grid is the choice that matters. A pipeline that logs a warning and continues produces a survey whose heights are wrong by metres, and the warning is in a log nobody read.

Confirming the grid was actually applied

The undulation range is the cheapest confirmation, and it is definitive. A grid-based conversion over a site of any size produces a range of at least millimetres; a fallback based on a constant or a low-order model produces a range of exactly zero or something implausibly smooth.

import numpy as np


def confirm_grid_applied(undulation_m: np.ndarray, site_extent_m: float) -> dict:
    """Does the undulation vary the way a real geoid model would?

    A geoid surface varies smoothly but measurably: roughly a millimetre per
    hundred metres at the low end. A range of exactly zero means a constant
    was applied; an implausibly large range means the wrong model.
    """
    u = np.asarray(undulation_m)
    u = u[np.isfinite(u)]
    if u.size < 10:
        return {"note": "too few points to judge"}

    spread = float(np.ptp(u))
    expected_min = site_extent_m * 1e-5          # about 1 mm per 100 m
    return {"mean_m": float(np.mean(u)), "range_m": spread,
            "expected_minimum_m": expected_min,
            "grid_applied": spread > expected_min,
            "note": ("the undulation varies across the site, as a grid would"
                     if spread > expected_min else
                     "the undulation is constant — a grid was not applied")}
Three ways a geoid transformation silently does nothing useful Three rows. The grid file not being found leaves the transformer falling back to a null transformation, which returns the input heights unchanged and reports success, so the output is ellipsoidal heights labelled orthometric. The wrong grid for the region returns a separation that is smoothly wrong by metres rather than obviously wrong, because geoid models differ regionally and all of them produce plausible-looking numbers. The transformation applied in the wrong direction subtracts where it should add, doubling the error rather than removing it, and is detectable only against a benchmark. the grid was not found a null transformation returns the input unchanged, reporting success the wrong grid for the region smoothly wrong by metres; every geoid model looks plausible applied in the wrong direction doubles the error instead of removing it All three succeed. Only a benchmark comparison distinguishes them from a correct result.

Figure 3 — Three silent failures with one shared remedy.

Edge-case matrix

Situation Symptom Handling
EPSG:4326 as source Heights unchanged Use EPSG:4979
always_xy omitted Coordinates swapped Always pass it, or match the axis order
Grid missing Silent fallback Check the transformer group and raise
Transformer in a loop Thousands of times slower Build once, reuse
Points outside the grid extent Fallback or NaN Check the site is within coverage
Network fetching enabled Result depends on a download Bake grids in, PROJ_NETWORK=OFF
Compound CRS unsupported downstream Product loses the vertical datum Write it as a tag as well
Dynamic datum, no epoch Default epoch used Pass both epochs explicitly

The grid-extent row is worth a check on any site near a national boundary. National geoid models cover a national area, and a survey that straddles an edge has some points converted with the grid and some with a fallback — producing a step in the middle of the survey that looks like a reconstruction fault.

import numpy as np


def check_grid_coverage(undulation_m: np.ndarray, positions_xy: np.ndarray,
                        *, max_step_m: float = 0.1) -> dict:
    """Look for a discontinuity in the undulation, which means a coverage edge."""
    u = np.asarray(undulation_m)
    order = np.argsort(positions_xy[:, 0])
    steps = np.abs(np.diff(u[order]))
    worst = float(steps.max()) if steps.size else 0.0
    return {"max_step_m": worst,
            "coverage_edge_suspected": worst > max_step_m,
            "note": ("undulation varies smoothly" if worst <= max_step_m else
                     "a step in the undulation suggests part of the site falls "
                     "outside the geoid grid's coverage")}

Verification snippet

import numpy as np


def round_trip_check(transformer, inverse_transformer,
                     lon: np.ndarray, lat: np.ndarray, h: np.ndarray,
                     *, tol_m: float = 0.001) -> dict:
    """Transform and transform back; the result must return to the input.

    A round trip catches an axis-order error, a wrong CRS pair and a
    non-invertible fallback in one test, and it needs no external reference.
    """
    x, y, z = transformer.transform(lon, lat, h)
    lon2, lat2, h2 = inverse_transformer.transform(x, y, z)

    dlon = np.abs(np.asarray(lon2) - np.asarray(lon)).max()
    dlat = np.abs(np.asarray(lat2) - np.asarray(lat)).max()
    dh = float(np.abs(np.asarray(h2) - np.asarray(h)).max())

    return {"max_lon_error_deg": float(dlon), "max_lat_error_deg": float(dlat),
            "max_height_error_m": dh,
            "ok": dh < tol_m and dlon < 1e-9 and dlat < 1e-9,
            "note": ("round trip is clean" if dh < tol_m else
                     "the transformation is not invertible to tolerance — check the "
                     "CRS pair and the grid")}
Geoid undulation across a site, with a grid and with a constant Two profiles of geoid undulation across a four kilometre site. The grid-based conversion produces a smooth curve varying from forty-eight point two metres at one end to forty-eight point nine at the other, a range of seventy centimetres. The constant approximation is a flat line at forty-eight point five metres. A shaded region shows the resulting height error, reaching plus and minus thirty-five centimetres at the two ends of the site — a tilt rather than an offset. 01 km 2 km3 km4 km position across the site undulation geoid grid — 48.2 to 48.9 m constant 48.5 m −35 cm +35 cm A constant undulation produces a tilt, which is harder to spot than an offset.

Figure 2 — What the grid is doing, on a site of ordinary size.

Using it at survey scale

Bulk conversion is where the transformer’s cost structure matters. Constructing one is milliseconds; transforming a million points through a constructed one is also milliseconds, because pyproj passes arrays to PROJ in one call.

import numpy as np


def convert_survey(transformer, coords: np.ndarray, *, chunk: int = 1_000_000) -> np.ndarray:
    """Convert a large array of positions in chunks, reusing one transformer.

    Chunking bounds the peak memory rather than the runtime — PROJ handles a
    million points per call comfortably, and the chunk size exists so a
    hundred-million-point cloud does not need its own copy in memory.
    """
    out = np.empty_like(coords, dtype=np.float64)
    for start in range(0, len(coords), chunk):
        end = min(start + chunk, len(coords))
        x, y, z = transformer.transform(coords[start:end, 0],
                                        coords[start:end, 1],
                                        coords[start:end, 2])
        out[start:end] = np.column_stack([x, y, z])
    return out

When to escalate

  • The required grid is not distributed with PROJ. Some national grids are licensed separately. Obtain them, install them into the PROJ data directory, and bake them into the processing image.
  • The site straddles a grid boundary. Two models over the two halves will differ; choose one and state it, or split the deliverable.
  • The transformation is correct and a downstream tool loses the vertical datum. Write it as a raster tag as well as in the CRS, so the information survives a tool that only understands horizontal systems.
  • A client’s data was converted with an older model version. The difference is legitimate on both sides. Record which version produced which product rather than reconciling to a single number.

← Geoid Models and Vertical Datum Automation