Should I add or subtract the geoid undulation?

Always subtract: H = h - N. Adding a negative CONUS undulation lowers H below h and biases every height by twice N while still returning a plausible number. Assert the expected sign relationship for the region as a guard.

Why must I combine the geoid uncertainty with the GNSS uncertainty?

H depends on two independent measurements, h and N, so their variances add in quadrature. Reporting only the GNSS vertical uncertainty understates the budget and can fail an agency's independent verification.

What goes wrong when I mix geoid models or tide systems?

Different geoid releases and different tide systems give undulations that differ by centimetres to decimetres for the same point. Mixing them injects a systematic bias, so record the model and tide system with every reduced height.

Applying Geoid Undulation for Orthometric Heights

Reducing a GNSS ellipsoidal height to an orthometric height is a single subtraction — $H = 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 $h$ and a geoid undulation $N$ — either supplied directly or bilinearly interpolated from a small grid — and returns the orthometric height $H$ together with a propagated uncertainty, so the result is not just a number but a defensible number.

Figure — in the conterminous United States $N$ is negative, so subtracting it raises $H$ above $h$ .

The Relation and Its Uncertainty

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

H = h - N

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

\sigma_H = \sqrt{\sigma_h^2 + \sigma_N^2}

where $\sigma_h$ is the GNSS vertical uncertainty and $\sigma_N$ combines the model’s stated accuracy with the additional error introduced by interpolating between nodes. Ignoring $\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 $N$ is the other hazard: in the conterminous United States the geoid lies below the ellipsoid, so $N$ runs roughly $-8\,\text{m}$ to $-33\,\text{m}$ , and subtracting that negative value correctly makes $H$ larger than $h$ . This is the opposite visual sense of the positive- $N$ 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 $N$ from a small in-memory grid (or accepts a directly supplied $N$ ), applies $H = 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 - 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\,\text{m}$ at latitude $38.5000°$ N, longitude $-98.0000°$ (positive-east). We interpolate $N$ 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\,\text{m}$ — over $28\,\text{m}$ higher than the raw ellipsoidal figure, precisely because the negative undulation is subtracted. The combined uncertainty of $0.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.

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.

Geoid and vertical datum transformations in Python — the parent guide on reading geoid grids and the full height relationship.
Understanding NTv2 grid shift files in Python — the horizontal counterpart of bilinear grid interpolation.
Geodetic conversion math: ellipsoid to Cartesian — where the reduced coordinates feed a geocentric mapping.
Algorithmic Math & Geodetic Workflows — the parent reference for the deterministic geodetic pipeline.

Applying Geoid Undulation for Orthometric Heights

The Relation and Its Uncertainty #

Complete Runnable Implementation #

Parameter Reference #

Worked Example #

Validation Check #

Common Mistakes #

Related #