CoordTransform Cadastral & High-Precision Python

Applying Geoid Undulation for Orthometric Heights

Reducing a GNSS ellipsoidal height to an orthometric height is a single subtraction — H=h−NH = h - N — but it is the subtraction that carries the entire vertical datum of a cadastral deliverable, and getting its sign or its uncertainty wrong quietly corrupts every height in the file. This routine is the height-reduction step of Geoid and Vertical Datum Transformations in Python, itself part of the broader Algorithmic Math & Geodetic Workflows pipeline. Here we isolate one self-contained function: it takes an ellipsoidal height hh and a geoid undulation NN — either supplied directly or bilinearly interpolated from a small grid — and returns the orthometric height HH together with a propagated uncertainty, so the result is not just a number but a defensible number.

Reducing ellipsoidal height by a negative undulation A vertical axis with three marks. The ellipsoid reference is near the bottom, the geoid mark lies below it because N is negative in the conterminous United States, and the terrain mark is at the top. The span from ellipsoid to terrain is h; the span from geoid to terrain is H, which is larger than h because a negative N is subtracted. The relation H equals h minus N is shown. terrain ellipsoid geoid (N < 0) h H H = h − N

Figure — in the conterminous United States NN is negative, so subtracting it raises HH above hh.

The Relation and Its Uncertainty

The orthometric height is the ellipsoidal height minus the geoid undulation:

Height error produced by four common mistakes Bar chart of vertical error in metres from four mistakes: using the wrong sign for the undulation, which doubles it, at 86 metres; using a global model where a national one is published, at 0.35 metres; interpolating nearest-neighbour rather than bilinear on a 1 by 1 minute grid, at 0.08 metres; and ignoring the model epoch, at 0.02 metres. The first is catastrophic and obvious, the rest are quiet. 0.01 0.1 1 10 100 m 86 sign flipped 0.35 global vs national 0.08 nearest-neighbour 0.02 epoch ignored

Figure — four ways to get an orthometric height wrong, and how far each one moves it.

H=h−NH = h - N

Because hh and NN are independent measurements, their variances add in quadrature. The standard uncertainty of the reduced height is

σH=σh2+σN2\sigma_H = \sqrt{\sigma_h^2 + \sigma_N^2}

where σh\sigma_h is the GNSS vertical uncertainty and σN\sigma_N combines the model’s stated accuracy with the additional error introduced by interpolating between nodes. Ignoring σN\sigma_N — treating the geoid as if it were exact — understates the vertical error budget and is a common reason a height passes an internal check but fails an agency’s independent verification. The sign of NN is the other hazard: in the conterminous United States the geoid lies below the ellipsoid, so NN runs roughly −8 m-8\,\text{m} to −33 m-33\,\text{m}, and subtracting that negative value correctly makes HH larger than hh. This is the opposite visual sense of the positive-NN North Atlantic case shown in the parent guide, yet the same equation governs both.

Complete Runnable Implementation

The function below is self-contained: it bilinearly interpolates NN from a small in-memory grid (or accepts a directly supplied NN), applies H=h−NH = h - N with deterministic rounding, and propagates the combined uncertainty. It depends only on numpy, dataclasses, and the standard library.

from __future__ import annotations

import math
from dataclasses import dataclass

import numpy as np


@dataclass(frozen=True)
class OrthometricResult:
    H_orthometric_m: float
    undulation_m: float
    sigma_H_m: float


def bilinear_undulation(
    grid: np.ndarray,
    lat0: float,
    lon0: float,
    dlat: float,
    dlon: float,
    lat: float,
    lon: float,
) -> float:
    """Bilinearly interpolate geoid undulation N (metres) from a regular grid.

    grid has shape (nrows, ncols), row 0 at the lower-left origin (lat0, lon0),
    rows running south-to-north, columns west-to-east. Raises if out of extent.
    """
    nrows, ncols = grid.shape
    lat_max = lat0 + (nrows - 1) * dlat
    lon_max = lon0 + (ncols - 1) * dlon
    if not (lat0 <= lat <= lat_max and lon0 <= lon <= lon_max):
        raise ValueError("point outside geoid grid extent")

    fi = (lat - lat0) / dlat
    fj = (lon - lon0) / dlon
    r0 = min(int(math.floor(fi)), nrows - 2)   # clamp so r0+1 stays in range
    c0 = min(int(math.floor(fj)), ncols - 2)
    ty = fi - r0                                # in-cell offset, north positive
    tx = fj - c0                                # in-cell offset, east positive

    n00 = float(grid[r0, c0]);     n10 = float(grid[r0, c0 + 1])
    n01 = float(grid[r0 + 1, c0]); n11 = float(grid[r0 + 1, c0 + 1])
    return (n00 * (1.0 - tx) * (1.0 - ty) + n10 * tx * (1.0 - ty)
            + n01 * (1.0 - tx) * ty + n11 * tx * ty)


def apply_undulation(
    h_ellipsoidal_m: float,
    undulation_m: float,
    sigma_h_m: float = 0.020,
    sigma_N_m: float = 0.015,
    ndigits: int = 4,
) -> OrthometricResult:
    """Reduce an ellipsoidal height to orthometric via H = h - N and propagate
    the combined vertical uncertainty in quadrature.

    Parameters
    ----------
    h_ellipsoidal_m : GNSS ellipsoidal height (m).
    undulation_m    : signed geoid undulation N (m); negative in CONUS.
    sigma_h_m       : GNSS vertical standard uncertainty (m).
    sigma_N_m       : geoid model + interpolation standard uncertainty (m).
    ndigits         : rounding precision (round-half-to-even).
    """
    big_h = round(h_ellipsoidal_m - undulation_m, ndigits)     # H = h - N
    sigma_H = round(math.hypot(sigma_h_m, sigma_N_m), ndigits)  # sqrt(sh^2 + sN^2)
    return OrthometricResult(
        H_orthometric_m=big_h,
        undulation_m=round(undulation_m, ndigits),
        sigma_H_m=sigma_H,
    )

Parameter Reference

Name Type Units Valid range Notes
h_ellipsoidal_m float metres any finite GNSS height above the ellipsoid, never a levelled height
undulation_m (N) float metres ≈ −110 to +90 signed; negative across the conterminous United States
sigma_h_m float metres ≥ 0 GNSS vertical standard uncertainty (1σ)
sigma_N_m float metres ≥ 0 model accuracy combined with interpolation error
ndigits int — ≥ 0 deterministic rounding precision
H_orthometric_m float metres any finite reduced orthometric height, H=h−NH = h - N
sigma_H_m float metres ≥ 0 combined vertical uncertainty in quadrature

Worked Example

Take a GNSS observation in central Kansas, where GEOID18 undulations are strongly negative. The receiver reports h=152.340 mh = 152.340\,\text{m} at latitude 38.5000°38.5000° N, longitude −98.0000°-98.0000° (positive-east). We interpolate NN from a small grid and reduce the height.

# Small grid near the point: all nodes at -28.740 m for a runnable demo.
lat0, lon0 = 38.4, -98.1
dlat = dlon = 1.0 / 60.0                       # one arc-minute
demo = np.full((13, 13), -28.740, dtype=np.float64)

n = bilinear_undulation(demo, lat0, lon0, dlat, dlon, lat=38.5, lon=-98.0)
res = apply_undulation(h_ellipsoidal_m=152.340, undulation_m=n)
print(res.undulation_m, res.H_orthometric_m, res.sigma_H_m)
# -28.74 181.08 0.025   ->  H = 152.340 - (-28.740) = 181.080 m

The reduced orthometric height is 181.080 m181.080\,\text{m} — over 28 m28\,\text{m} higher than the raw ellipsoidal figure, precisely because the negative undulation is subtracted. The combined uncertainty of 0.025 m0.025\,\text{m} reflects both the GNSS height noise and the geoid contribution, not the GNSS receiver alone.

Validation Check

Gate the result against a vertical control tolerance and confirm the sign behaved as expected for the region before the height is committed to a deliverable.

Looking up an undulation, including the out-of-extent case A position enters the lookup. The grid extent is checked first. Inside the extent, the four surrounding nodes are read and bilinearly interpolated to give the undulation. Outside it, the routine raises: extrapolating a geoid model past its published extent produces a number with no defensible uncertainty. A node carrying the model no-data sentinel is treated the same way as being outside the extent. Position lat, lon Inside the grid extent? Read the 4 surrounding nodes Raise — outside extent no defensible uncertainty Any node is no-data? Bilinear interpolation undulation N in metres yes no no yes

Figure — the lookup path, including the case most implementations forget.

assert math.isfinite(res.H_orthometric_m), "non-finite orthometric height"
# In CONUS N < 0, so H must exceed h; a smaller H means an inverted sign.
assert res.H_orthometric_m > 152.340, "negative-N region: H should exceed h"
assert res.sigma_H_m >= 0.020, "uncertainty cannot be below the GNSS component"
CONTROL_H = 181.075
assert abs(res.H_orthometric_m - CONTROL_H) < 0.030, "residual exceeds 30 mm gate"

Common Mistakes

Adding N instead of subtracting it
The relation is H = h - N, never h + N. Adding a negative CONUS undulation lowers H below h and shifts every height by 2N — well over 50 m in some regions — while still returning a finite, plausible-looking number. Assert the expected sign relationship for the region as a guard, as the validation check above does.
Mixing geoid models or tide systems
A GEOID12B undulation and a GEOID18 undulation for the same point differ by a real, published amount, and models referenced to different tide systems (mean-tide, zero-tide, tide-free) differ by a further few centimetres. Interpolating h from one datum against N from an inconsistent model injects a systematic bias. Record the model and tide system with every reduced height.
Ignoring the interpolation uncertainty in N
Treating the geoid as exact and reporting only the GNSS vertical uncertainty understates the true budget. The undulation carries both the model's stated accuracy and the error of interpolating between nodes; combine them into sigma_N_m and propagate in quadrature, or an internally clean height will fail an agency's independent check.

Frequently Asked Questions

Is the sign of the undulation the same everywhere?

No, which is exactly why the sign convention has to be checked rather than assumed. The geoid sits above the ellipsoid over much of Europe and below it over much of the Pacific, so an implementation that gets the sign backwards will look correct in one region and be out by twice the undulation in another. Test against a published control point in the region you actually work in, and prefer a model interface that names the convention in its documentation.

How accurate is a modern geoid model?

For a recent national model over well-surveyed terrain, one to three centimetres is typical; global models are usually good to a few decimetres. That uncertainty belongs in the height budget, and it is normally the largest single term in a GNSS-derived orthometric height — larger than the GNSS observation itself. Quoting a height to the millimetre while omitting the model uncertainty overstates what was measured.

Can I interpolate the undulation with nearest-neighbour?

Not for survey work. On a one-arc-minute model, nearest-neighbour sampling can be out by several centimetres between nodes, which is the same order as the model’s own uncertainty and entirely avoidable. Bilinear interpolation over the four surrounding nodes costs nothing and removes most of it; bicubic buys a little more where the surface is strongly curved.