Geoid and Vertical Datum Transformations in Python

Every GNSS receiver reports a height above the reference ellipsoid, yet no cadastral plan, flood model, or levelling network is graded in ellipsoidal metres — they are graded in orthometric height above the geoid, the physical surface that water actually follows. Reconciling those two vertical worlds is the height dimension of Algorithmic Math & Geodetic Workflows, and it turns on a single quantity: the geoid undulation $N$ , the separation between the ellipsoid and the geoid at each point. This page narrows the workflow to one concrete operation — converting an ellipsoidal height $h$ from a GNSS observation into an orthometric height $H$ by reading a gridded geoid model, bilinearly interpolating $N$ at the point, and applying the height relation with the correct sign. The horizontal analogue of this grid interpolation — reading a fixed-layout node array and blending the four surrounding values — is the same discipline used when understanding NTv2 grid shift files in Python; here the surface being interpolated is a height separation rather than a horizontal shift, but the interpolation algebra is identical.

Figure — the ellipsoid, geoid, and terrain at one point: $h$ is measured to the ellipsoid, $H$ to the geoid, and $N$ is their separation.

The Height Relation and Its Sign

The three heights are collinear along the ellipsoidal normal, so they combine by simple subtraction. The orthometric height is the ellipsoidal height minus the geoid undulation:

H = h - N

where $h$ is the GNSS-derived ellipsoidal height (metres above the reference ellipsoid, e.g. GRS80/WGS84), $N$ is the geoid undulation at the point, and $H$ is the orthometric height above the geoid. The undulation carries a sign: $N > 0$ where the geoid lies above the ellipsoid (much of the North Atlantic and Europe, where $N$ reaches roughly $+50\,\text{m}$ ) and $N < 0$ where it lies below (the conterminous United States, where GEOID18 undulations run roughly $-8\,\text{m}$ to $-33\,\text{m}$ ). Because the sign of $N$ flips between regions, the single most damaging error in vertical transformation is adding $N$ where you should subtract — a mistake that displaces every height by $2N$ , tens of metres, while still returning a plausible-looking number. The self-contained function that isolates exactly this sign discipline is built in applying geoid undulation for orthometric heights.

Strictly, $N$ is the geoid undulation used with orthometric height, while the closely related height anomaly $\zeta$ pairs with normal height in a Molodensky-type height system. Regular gridded models such as the U.S. GEOID series distribute $N$ directly on a latitude/longitude mesh; quasi-geoid models distribute $\zeta$ . The interpolation machinery below is agnostic to which quantity the grid stores — but the vertical datum you claim in your deliverable must match the quantity the grid was built to deliver.

The Geoid Model Grid

A geoid model is a regular grid of undulation values sampled at a fixed angular spacing in latitude and longitude. The GTX format (used by NOAA’s geoid models) and the newer GGXF container both reduce to the same logical structure: a lower-left origin $(\phi_0, \lambda_0)$ , a node spacing $(\Delta\phi, \Delta\lambda)$ , a row and column count, and a dense array of undulation values in metres stored row by row. Reading a value at an arbitrary point means locating the grid cell that encloses it and blending its four corner nodes.

Given a query point $(\phi, \lambda)$ , the fractional cell coordinates are

i = \frac{\phi - \phi_0}{\Delta\phi}, \qquad j = \frac{\lambda - \lambda_0}{\Delta\lambda}

with integer parts $(i_0, j_0)$ giving the lower-left node of the enclosing cell and fractional parts $(t_y, t_x) = (i - i_0,\; j - j_0)$ giving the position inside it. Writing the four surrounding undulations as $N_{00}$ (lower-left), $N_{10}$ (lower-right), $N_{01}$ (upper-left), and $N_{11}$ (upper-right), the bilinear interpolation is

N(\phi,\lambda) = N_{00}(1-t_x)(1-t_y) + N_{10}\,t_x(1-t_y) + N_{01}(1-t_x)t_y + N_{11}\,t_x t_y

This is the standard four-point bilinear kernel: linear in each axis, exact at the nodes, and continuous across cell boundaries. It is adequate for models whose spacing is one arc-minute or finer, where the geoid is smooth relative to the node interval; coarser grids or steep gravimetric features may warrant a bicubic kernel, but bilinear is the interoperable baseline that agencies expect and that matches how PROJ itself samples vertical grids.

Step-by-Step Implementation

The implementation uses only numpy, struct, dataclasses, logging, and the standard library so that every precision and sign decision stays explicit and auditable. Build it in four steps; each is a runnable, type-hinted unit.

Step 1 — Model the grid and parse a header

We represent the model as an immutable header plus a float64 undulation array. The reader below decodes a GTX-style binary header (lower-left latitude and longitude, latitude and longitude spacing, row and column counts) and then loads the big-endian float32 payload, promoting it to float64 before any interpolation touches it.

from __future__ import annotations

import logging
import struct
from dataclasses import dataclass
from pathlib import Path

import numpy as np

logger = logging.getLogger("geoid")

GTX_HEADER_STRUCT = ">dddii"   # lat0, lon0, dlat, dlon (deg), nrows, ncols
GTX_HEADER_SIZE = struct.calcsize(GTX_HEADER_STRUCT)   # 40 bytes


@dataclass(frozen=True)
class GeoidGrid:
    """A regular geoid-undulation grid, values in metres, promoted to float64.

    Nodes are stored south-to-north by row and west-to-east by column, with
    row 0 at the lower-left origin (lat0, lon0), matching the GTX convention.
    """

    lat0: float          # latitude of the lower-left node (deg)
    lon0: float          # longitude of the lower-left node (deg, positive east)
    dlat: float          # node spacing in latitude (deg, > 0)
    dlon: float          # node spacing in longitude (deg, > 0)
    values: np.ndarray   # shape (nrows, ncols), float64, undulation N in metres

    @property
    def nrows(self) -> int:
        return int(self.values.shape[0])

    @property
    def ncols(self) -> int:
        return int(self.values.shape[1])

    @property
    def lat_max(self) -> float:
        return self.lat0 + (self.nrows - 1) * self.dlat

    @property
    def lon_max(self) -> float:
        return self.lon0 + (self.ncols - 1) * self.dlon


def load_gtx(path: str | Path) -> GeoidGrid:
    """Read a GTX-style geoid grid into a GeoidGrid with float64 values."""
    raw = Path(path).read_bytes()
    lat0, lon0, dlat, dlon, nrows, ncols = struct.unpack_from(
        GTX_HEADER_STRUCT, raw, 0
    )
    count = nrows * ncols
    # GTX payload is big-endian float32; promote to float64 so the on-disk
    # single precision never leaks into the height budget during interpolation.
    payload = np.frombuffer(
        raw, dtype=np.dtype(">f4"), count=count, offset=GTX_HEADER_SIZE
    ).astype(np.float64).reshape(nrows, ncols)
    return GeoidGrid(lat0=lat0, lon0=lon0, dlat=dlat, dlon=dlon, values=payload)

Step 2 — Locate the enclosing cell and guard the extent

Before interpolating, the query must be inside the grid. A point outside the model’s footprint has no defensible undulation, and extrapolating one silently fabricates height. This step returns the lower-left node indices and the in-cell fractional offsets, raising if the point falls outside.

def locate_cell(grid: GeoidGrid, lat: float, lon: float) -> tuple[int, int, float, float]:
    """Return (row0, col0, ty, tx): the lower-left node of the enclosing cell
    and the fractional offsets within it. Raises if the point is out of extent."""
    if not (grid.lat0 <= lat <= grid.lat_max):
        raise ValueError(f"latitude {lat} outside grid [{grid.lat0}, {grid.lat_max}]")
    if not (grid.lon0 <= lon <= grid.lon_max):
        raise ValueError(f"longitude {lon} outside grid [{grid.lon0}, {grid.lon_max}]")

    fi = (lat - grid.lat0) / grid.dlat       # fractional row index
    fj = (lon - grid.lon0) / grid.dlon       # fractional column index
    row0 = min(int(np.floor(fi)), grid.nrows - 2)   # clamp so row0+1 is valid
    col0 = min(int(np.floor(fj)), grid.ncols - 2)   # clamp so col0+1 is valid
    ty = fi - row0                            # 0..1 within the cell, north+
    tx = fj - col0                            # 0..1 within the cell, east+
    return row0, col0, ty, tx

Step 3 — Bilinearly interpolate the undulation

With the cell located, blend its four corner undulations using the kernel derived above. The function checks for the model’s null sentinel so a query that lands on an unmodelled node is rejected rather than averaged into a fictitious value.

NULL_SENTINEL = -88.8888   # common GTX flag for an undefined node


def interpolate_undulation(grid: GeoidGrid, lat: float, lon: float) -> float:
    """Bilinearly interpolate the geoid undulation N (metres) at (lat, lon)."""
    row0, col0, ty, tx = locate_cell(grid, lat, lon)
    n00 = grid.values[row0, col0]         # lower-left
    n10 = grid.values[row0, col0 + 1]     # lower-right
    n01 = grid.values[row0 + 1, col0]     # upper-left
    n11 = grid.values[row0 + 1, col0 + 1] # upper-right

    corners = np.array([n00, n10, n01, n11], dtype=np.float64)
    if np.any(np.isclose(corners, NULL_SENTINEL, atol=1e-3)):
        raise ValueError("query cell contains an undefined geoid node")

    # Four-point bilinear kernel: exact at nodes, continuous across cell edges.
    return float(
        n00 * (1.0 - tx) * (1.0 - ty)
        + n10 * tx * (1.0 - ty)
        + n01 * (1.0 - tx) * ty
        + n11 * tx * ty
    )

Step 4 — Apply the height relation and log the conversion

The final step applies $H = h - N$ , rounds to a declared precision, and emits a structured audit record naming the model. Every certified height conversion should be reconstructable from its log line alone: the input ellipsoidal height, the interpolated undulation, the model identity, and the resulting orthometric height.

import json
from dataclasses import asdict


@dataclass(frozen=True)
class HeightConversion:
    lat: float
    lon: float
    h_ellipsoidal_m: float
    undulation_m: float
    H_orthometric_m: float
    geoid_model: str


def to_orthometric(
    grid: GeoidGrid,
    lat: float,
    lon: float,
    h_ellipsoidal_m: float,
    geoid_model: str,
    ndigits: int = 4,
) -> HeightConversion:
    """Convert an ellipsoidal height to orthometric via H = h - N and log it."""
    n = interpolate_undulation(grid, lat, lon)
    big_h = round(h_ellipsoidal_m - n, ndigits)   # H = h - N (round-half-to-even)
    result = HeightConversion(
        lat=lat, lon=lon,
        h_ellipsoidal_m=h_ellipsoidal_m,
        undulation_m=round(n, ndigits),
        H_orthometric_m=big_h,
        geoid_model=geoid_model,
    )
    logger.info("geoid_convert %s", json.dumps(asdict(result)))
    return result

Parameter and Return-Value Reference

Name	Type	Units	Valid range	Cadastral significance
`lat` / `lon`	`float`	degrees	inside grid extent	positive-east longitude; a point outside the footprint has no defensible $N$
`dlat` / `dlon`	`float`	degrees	> 0	node spacing; drives interpolation error, must match the model header exactly
`h_ellipsoidal_m`	`float`	metres	any finite	GNSS height above the ellipsoid — never a levelled height
`undulation_m` (`N`)	`float`	metres	≈ −110 to +90 globally	signed geoid–ellipsoid separation; its sign convention governs the whole result
`H_orthometric_m`	`float`	metres	any finite	the certified height above the geoid, $H = h - N$
`geoid_model`	`str`	—	e.g. `GEOID18`	model identity must be logged; mixing models corrupts the vertical datum
`ty` / `tx`	`float`	—	0 ≤ t ≤ 1	in-cell fractional position; outside [0,1] means a cell-location bug
`interpolate_undulation()` →	`float`	metres	finite or raises	raises on out-of-extent or null-node queries, never extrapolates

Worked Example: a North Sea GNSS Height

Consider a GNSS observation over the southern North Sea where the geoid lies well above the GRS80 ellipsoid, so $N$ is strongly positive. The receiver reports an ellipsoidal height of $142.500\,\text{m}$ at latitude $54.0000°$ N, longitude $2.0000°$ E. Assume a one-arc-minute model whose interpolated undulation at the point is $N = 48.200\,\text{m}$ .

# Synthesise a small one-arc-minute grid around the point for a runnable demo.
lat0, lon0 = 53.9, 1.9
dlat = dlon = 1.0 / 60.0                     # one arc-minute in degrees
rows, cols = 13, 13
base = np.full((rows, cols), 48.200, dtype=np.float64)  # ~48.2 m undulation
grid = GeoidGrid(lat0=lat0, lon0=lon0, dlat=dlat, dlon=dlon, values=base)

conv = to_orthometric(grid, lat=54.0000, lon=2.0000,
                      h_ellipsoidal_m=142.500, geoid_model="DEMO_NSea")
print(conv.undulation_m, conv.H_orthometric_m)
# 48.2 94.3   ->  H = 142.500 - 48.200 = 94.300 m

The orthometric height, $94.300\,\text{m}$ , is what appears on a survey plan or a flood datum, while the raw GNSS figure of $142.500\,\text{m}$ is meaningless to a levelling crew. The $48.2\,\text{m}$ separation is entirely explained by the geoid sitting above the ellipsoid in this region — a positive $N$ that would be a negative $N$ of similar magnitude if the same computation were run over the central United States. The complementary horizontal computation that carries these same coordinates into a geocentric Cartesian frame is set out in geodetic conversion math: ellipsoid to Cartesian.

Verification and Residual Analysis

A vertical transformation is only trustworthy once it is pinned two ways: a round-trip identity that proves the interpolation and height algebra are self-consistent, and a comparison against a published benchmark whose orthometric height is independently known. The round-trip recovers $h$ from $H$ by re-adding the same $N$ ; the residual against a control benchmark is the signed height difference compared to the vertical tolerance.

def verify_height(
    conv: HeightConversion,
    grid: GeoidGrid,
    control_H_m: float,
    tolerance_m: float = 0.030,   # 30 mm typical vertical control gate
) -> dict[str, object]:
    """Round-trip the conversion and compare H against a control benchmark."""
    # Round-trip: re-add the interpolated N and recover the ellipsoidal height.
    n = interpolate_undulation(grid, conv.lat, conv.lon)
    h_recovered = conv.H_orthometric_m + n            # inverse of H = h - N
    round_trip_err = abs(h_recovered - conv.h_ellipsoidal_m)
    assert round_trip_err < 1e-6, "height algebra is not self-consistent"

    residual = conv.H_orthometric_m - control_H_m     # signed vertical residual
    record = {
        "geoid_model": conv.geoid_model,
        "residual_m": round(residual, 4),
        "tolerance_m": tolerance_m,
        "round_trip_err_m": round(round_trip_err, 9),
        "passed": abs(residual) <= tolerance_m,
    }
    logger.info("geoid_verify %s", json.dumps(record))
    if not record["passed"]:
        raise ValueError(f"vertical residual {residual:.4f} m exceeds {tolerance_m} m")
    return record


rec = verify_height(conv, grid, control_H_m=94.312)
assert rec["passed"]        # residual -0.012 m, inside the 30 mm gate

The signed residual matters: a systematic bias of one sign across many benchmarks usually means the wrong model or a tide-system offset (below), whereas scatter of both signs is ordinary interpolation and levelling noise. Emitting the model name, residual, tolerance, and round-trip error is the minimum audit payload a vertical deliverable needs. Propagating the interpolation uncertainty formally into that residual budget is treated in the dedicated routine on applying geoid undulation for orthometric heights.

Troubleshooting and Gotchas

Every orthometric height is off by roughly twice the undulation

The sign of N was inverted — you added it instead of subtracting, or vice versa. The relation is H = h - N. In CONUS N is negative (about −8 to −33 m), so subtracting a negative correctly raises H above h; in the North Atlantic N is positive and H falls below h. A blanket offset of 2N is the fingerprint of a flipped sign.

Heights are consistently biased by a few centimetres to a decimetre

This is almost always the wrong geoid model or a superseded release: a GEOID12B height differs from a GEOID18 height by a real, published amount. Confirm the model identity logged in your audit record matches the datum your deliverable claims, and never interpolate one region's point through another region's grid.

A valid-looking coastal point raises an out-of-extent error

Geoid grids are published on a finite footprint and often carry -88.8888 null nodes over open ocean. If the query cell touches a null node or falls outside the header extent, reject it rather than extrapolate. Check the point against lat0, lon0, lat_max, and lon_max, and confirm longitude sign uses the model's positive-east convention.

GNSS heights disagree with a published model by a constant offset

Different geoid products are referenced to different tide systems — mean-tide, zero-tide, or tide-free — and the permanent-tide correction can reach a few centimetres, latitude-dependent. Mixing a tide-free model with mean-tide ellipsoidal heights injects a systematic bias. Record the tide system alongside the model name and reduce both quantities to the same system before differencing.

Frequently Asked Questions

Why not just let pyproj apply the vertical transformation?

For a supported CRS pair, PROJ does exactly this and you should use it in production. Implementing the grid read and bilinear blend yourself is what lets you audit the interpolation, attach a defensible uncertainty, and handle models PROJ has no operation for — the same reason cadastral pipelines sometimes read grids directly rather than through the transformer API.

Is bilinear interpolation accurate enough for survey-grade heights?

For one-arc-minute or finer geoid models the geoid is smooth relative to the node spacing, so bilinear interpolation holds well within a centimetre and matches how PROJ samples the same grids. Coarse or highly variable gravimetric geoids can justify a bicubic kernel, but bilinear is the interoperable baseline agencies expect.

What is the difference between geoid undulation and height anomaly?

The undulation N pairs with orthometric height on the geoid; the height anomaly ζ pairs with normal height on the quasi-geoid. Regular geoid grids store N, quasi-geoid grids store ζ. The interpolation is identical, but the vertical datum you claim must match the quantity the grid was built to deliver.

Applying geoid undulation for orthometric heights — the self-contained conversion function with sign discipline and uncertainty propagation.
Understanding NTv2 grid shift files in Python — the horizontal analogue of reading a fixed-layout grid and bilinearly blending nodes.
Geodetic conversion math: ellipsoid to Cartesian — the forward mapping that carries the same coordinates into a geocentric frame.
Error distribution modeling in Python — propagating the vertical residual into a full covariance budget.
Algorithmic Math & Geodetic Workflows — the parent reference covering the end-to-end geodetic computation pipeline.

Geoid and Vertical Datum Transformations in Python

The Height Relation and Its Sign #

The Geoid Model Grid #

Step-by-Step Implementation #

Step 1 — Model the grid and parse a header #

Step 2 — Locate the enclosing cell and guard the extent #

Step 3 — Bilinearly interpolate the undulation #

Step 4 — Apply the height relation and log the conversion #

Parameter and Return-Value Reference #

Worked Example: a North Sea GNSS Height #

Verification and Residual Analysis #

Troubleshooting and Gotchas #

Frequently Asked Questions #

Related #