Resolving a Vertical Datum Mismatch in GCP Heights

The horizontal residuals are two centimetres and the vertical residuals are all about −43 metres. Or all about −2.1 metres. In either case every control point agrees with every other, the plan position is excellent, and the surface sits at the wrong height by a constant.

That constant is the signature. A vertical datum mismatch is the only error in a survey that is genuinely uniform, and recognising it from that uniformity is faster than any amount of inspection of individual points.

Three vertical references, two ways to mix them

Ellipsoidal height is what a GNSS receiver computes directly: distance from the reference ellipsoid along its normal. It is purely geometric and has no relationship to where water flows.

Orthometric height is measured from the geoid — the equipotential surface mean sea level approximates — and is what every survey deliverable is expressed in. The difference between the two is the geoid undulation, and it ranges from about −105 m to +85 m worldwide.

A local vertical datum is a national realisation tied to specific tide gauges and benchmarks. It differs from a global geoid model by a smaller amount, typically a few decimetres to a couple of metres, and that difference is the second mixing failure.

Those two produce distinguishable magnitudes. A tens-of-metres offset is ellipsoidal against orthometric. A decimetre-to-two-metre offset is one orthometric system against another.

Three vertical references and the two offsets between them Three stacked reference surfaces in cross-section. At the bottom, the reference ellipsoid, drawn as a smooth line. Above it, a global geoid model, undulating and separated from the ellipsoid by tens of metres. Just above or below that, a national vertical datum, separated from the global geoid by a few decimetres to a couple of metres. A point on the terrain is shown with its height measured from each surface, and the two offsets are annotated with their characteristic magnitudes so a measured offset can be attributed to one or the other. global geoid national datum ellipsoid terrain geoid undulation N — tens of metres datum difference — decimetres to metres The magnitude of the offset says which pair of surfaces was mixed. Nothing else in a survey produces an error that is the same at every point.

Figure 1 — Two possible mixings with two distinguishable magnitudes. The measured offset is a diagnosis, not just a correction.

Minimal reproducible solution

Confirm uniformity first — a non-uniform vertical error is not a datum problem and applying a datum fix to it makes things worse.

import numpy as np


def classify_vertical_error(dz: np.ndarray, planimetric_extent_m: float) -> str:
    """Uniform, tilted, or neither, from the residuals alone."""
    dz = np.asarray(dz, dtype=float)
    if dz.size < 4:
        return "too few points to classify"

    mean, spread = float(dz.mean()), float(dz.std())
    if spread > abs(mean) * 0.25 and spread > 0.05:
        return (f"not uniform: mean {mean:.3f} m, spread {spread:.3f} m — "
                "this is not a datum offset")

    magnitude = abs(mean)
    if magnitude > 8.0:
        return (f"uniform {mean:+.2f} m — ellipsoidal height used where "
                "orthometric was expected, or vice versa")
    if magnitude > 0.08:
        return (f"uniform {mean:+.3f} m — two orthometric systems, or a "
                "national datum against a global geoid model")
    return f"uniform {mean:+.3f} m — within normal survey noise"

The spread test comes before the magnitude test on purpose. A block with a systematic tilt has a non-zero mean and a large spread, and treating it as an offset removes the average while leaving the tilt — producing a result that is worse at both ends and correct in the middle. That case belongs in diagnosing systematic tilt in GCP residuals.

With uniformity established, the fix is a transformation, not a subtraction. The undulation varies across a survey — by centimetres over a small site and by decimetres over a corridor — so a single constant is right at one point and wrong everywhere else.

import numpy as np
from pyproj import CRS, Transformer


def to_orthometric(lon: np.ndarray, lat: np.ndarray, h_ellipsoidal: np.ndarray,
                   vertical_crs: str = "EPSG:5773") -> np.ndarray:
    """Ellipsoidal to orthometric heights via a geoid grid.

    Compound source (horizontal + ellipsoidal) to compound target (horizontal +
    the named vertical datum). PROJ applies the grid; without it installed the
    transform silently returns the heights unchanged, so the result is checked.
    """
    src = CRS.from_user_input("EPSG:4979")            # WGS84 3D, ellipsoidal
    tgt = CRS.from_user_input(f"EPSG:4326+{vertical_crs.split(':')[1]}")
    tr = Transformer.from_crs(src, tgt, always_xy=True)

    _, _, h_ortho = tr.transform(lon, lat, h_ellipsoidal)
    h_ortho = np.asarray(h_ortho, dtype=float)

    if np.allclose(h_ortho, h_ellipsoidal, atol=1e-6):
        raise RuntimeError(
            f"heights unchanged — the geoid grid for {vertical_crs} is not "
            "installed. Set PROJ_NETWORK=ON or ship the grid; PROJ falls back "
            "to a horizontal-only transform without raising.")
    if not np.all(np.isfinite(h_ortho)):
        raise RuntimeError("some points fell outside the grid's area of use")
    return h_ortho

The unchanged-heights check is the guard that matters most, and it is easy to omit. A missing geoid grid does not raise; PROJ resolves a lower-accuracy path that leaves height alone, so the code runs, the numbers come back, and the correction was never applied. Asserting that the values moved is the only reliable detection.

Edge-case matrix

Situation Signature Correct action
Ellipsoidal used as orthometric Uniform, tens of metres Transform with a geoid grid
National datum vs global geoid Uniform, 0.1–2 m Transform between the two vertical CRS
Geoid grid missing Heights unchanged after transform Install the grid or enable network
Point outside the grid extent Height returns non-finite Different grid, or the point is misplaced
Antenna height not subtracted Uniform, 1.5–2.5 m Not a datum issue — instrument setup
Uniform offset plus a tilt Non-zero mean, large spread Fix the tilt; the offset may be a symptom
Two surveys, two datums Step between two groups Classify per group, not per project
Offset changes between epochs Small, growing Vertical plate motion; needs epochs

The antenna-height row is the one that costs the most time, because it mimics a small datum difference exactly. A rover pole height not entered, or entered and applied twice, produces a uniform offset of one to two and a half metres — the same range as a national datum difference — and the two are distinguished only by checking the field notes.

Verification snippet

After the transform, confirm that the applied correction matches the undulation the model predicts, and that the residual spread did not grow.

import numpy as np


def assert_correction_is_the_undulation(h_before: np.ndarray, h_after: np.ndarray,
                                        expected_n: np.ndarray,
                                        tol_m: float = 0.05) -> None:
    """The height change must equal the geoid undulation, point by point."""
    applied = h_before - h_after
    departure = np.abs(applied - expected_n)
    worst = float(departure.max())
    assert worst <= tol_m, (
        f"applied correction departs from the modelled undulation by up to "
        f"{worst:.3f} m — a different grid was used than the one expected")

    # The correction must not have added scatter: it is a smooth field.
    assert float(applied.std()) < float(expected_n.std()) * 1.5 + 0.01, (
        "the correction is noisier than the undulation field — check that "
        "points were not transformed individually with different operations")

The second assertion catches an unusual but real failure: transforming points one at a time, where PROJ may select a different operation per point if some fall outside a grid’s extent. The result is a correction field with steps in it, which is invisible in a mean and obvious in a standard deviation.

A constant subtraction against a grid transformation A survey corridor eight kilometres long, with the true geoid undulation varying smoothly from forty-two point one metres at one end to forty-two point nine at the other. Subtracting a single constant of forty-two point five metres leaves a residual error that is zero in the middle and four decimetres at each end, with opposite signs. Applying the geoid grid removes the undulation everywhere, leaving a flat residual near zero. A note observes that on a small compact site the two are indistinguishable, which is why the constant survives until it is used on a long job. distance along the corridor → residual height error 0 constant subtraction +0.4 m at the far end −0.4 m here geoid grid applied On a compact site the two agree to millimetres, which is why the constant survives until the first long corridor.

Figure 2 — Why the fix is a transformation. A constant is the correct answer at exactly one point on the site and drifts away from it in both directions.

When to escalate

  • The offset persists after a correct transformation. The declared source datum is wrong rather than the transformation. Check what the receiver was actually configured to output — a base station entered at the wrong height produces a uniform offset that no vertical transformation will remove, because the input was already wrong.
  • Two groups of control disagree by a constant. Two surveys or two instruments with different setups. Classify each group separately; a single mean across both describes neither.
  • The offset is uniform and small and matches the antenna height. Not a datum problem. Confirm against the field notes before transforming anything, because applying a geoid correction on top of an instrument-setup error leaves the survey wrong by their difference.

← Troubleshooting GCP and Coordinate Errors

Reading the cause from the magnitude A logarithmic scale of uniform vertical offset magnitude, from one centimetre to one hundred metres, with four labelled regions. Below eight centimetres is normal survey noise. From about one decimetre to two metres covers both a national datum difference and an unsubtracted antenna height, which overlap and must be separated by checking the field notes. From two to eight metres is unusual and often indicates two errors combined. Above eight metres is an ellipsoidal against orthometric mixing. A note states that the overlap region is the one that costs time. 1 cm 10 cm 1 m 10 m 100 m survey noise national datum or antenna height these overlap — check the field notes unusual ellipsoidal vs orthometric The magnitude narrows the cause to one or two candidates before anything is opened. Only the amber region is ambiguous, and it is resolved by a question rather than by a computation.

Figure 3 — The magnitude as a first-pass diagnosis. Three of the four regions identify the cause outright; the fourth narrows it to two candidates that a field note settles.