Checking Vertical Accuracy Against Benchmarks

Every internal check a photogrammetry pipeline performs is a consistency check: the residuals, the control fit, the overlap agreement. None of them can detect an error that shifted everything by the same amount, and a wrong vertical datum does exactly that.

The only test that catches it compares the survey against heights established independently — published benchmarks, a levelling run, or control measured by a method that does not share the survey’s assumptions. This page covers doing that comparison properly: which points to use, how to read the result, and how to state it in a way that survives a client’s own check. It completes geoid models and vertical datum automation.

What makes a comparison independent

A comparison is only as good as the independence of the reference, and three things can quietly destroy it.

The reference was used as control. A point that constrained the adjustment will agree with the survey by construction, whatever the datum. Reference points must be held out.

The reference came from the same GNSS session. A base station with a wrong height propagates the same error into both the control and the reference, and they agree perfectly while both are wrong.

The reference was converted with the same geoid model. If the published value was itself derived from an ellipsoidal measurement through the model under test, the comparison tests nothing about the model.

A genuinely independent reference is a levelled benchmark: a height established by spirit levelling from a national network, which never touches a satellite or a geoid model.

What each kind of reference can and cannot detect Four reference types against what they can detect. A control point used in the adjustment detects nothing, because it agrees by construction. A held-out checkpoint from the same GNSS session detects reconstruction error but not a datum or base-station error, because it shares them. A checkpoint converted with the same geoid model detects reconstruction and base-station error but not a model error. A levelled benchmark from a national network detects all three, because its height was established without a satellite or a geoid model. reconstruction base station geoid model control used in the adjustment checkpoint, same GNSS session checkpoint, same geoid model levelled benchmark nono no yesno no yesyes no yesyes yes Only the bottom row tests the datum, because only it was established without one. Every row above it shares an assumption with the survey and agrees with it for that reason.

Figure 1 — Independence is a property of how the reference was measured, not of whether it was held out.

Minimal reproducible solution

import numpy as np


def benchmark_comparison(survey_h: np.ndarray, published_h: np.ndarray,
                         *, published_sigma_m: float = 0.01) -> dict:
    """Compare survey heights against published benchmark values.

    Bias and scatter are reported separately because they mean different
    things: a bias is a datum or systematic error and a scatter is ordinary
    measurement noise, and a single RMSE conflates them.
    """
    d = np.asarray(survey_h, dtype=float) - np.asarray(published_h, dtype=float)
    d = d[np.isfinite(d)]
    if d.size < 2:
        return {"note": "at least two benchmarks are needed"}

    bias = float(np.mean(d))
    scatter = float(np.std(d, ddof=1))
    standard_error = scatter / np.sqrt(d.size)
    combined = float(np.hypot(standard_error, published_sigma_m))

    return {"benchmarks": int(d.size),
            "bias_m": bias, "scatter_m": scatter,
            "bias_uncertainty_m": combined,
            "bias_significant": abs(bias) > 2 * combined,
            "rmse_m": float(np.sqrt(np.mean(d ** 2))),
            "residuals_m": d.tolist()}

Folding the published values’ own uncertainty into the bias uncertainty is what keeps the conclusion honest. A measured bias of 2 cm against benchmarks known to 1 cm, from four points, is not significantly different from zero — and reporting it as a 2 cm datum error would send somebody looking for a fault that is not there.

Interpreting the result

The magnitude of the bias usually names its cause, which makes the diagnosis quick.

def diagnose_bias(bias_m: float, local_undulation_m: float,
                  antenna_height_m: float) -> str:
    """Name the likely cause from the size of the bias."""
    candidates = {
        "vertical datum (geoid undulation applied backwards)": 2 * local_undulation_m,
        "vertical datum (geoid undulation not applied)": local_undulation_m,
        "antenna height applied with the wrong sign": 2 * antenna_height_m,
        "antenna height not applied": antenna_height_m,
    }
    for cause, expected in sorted(candidates.items(), key=lambda kv: -abs(kv[1])):
        if abs(abs(bias_m) - abs(expected)) < max(0.1 * abs(expected), 0.02):
            return f"{cause} — the bias matches {expected:+.3f} m"
    if abs(bias_m) < 0.05:
        return "no significant systematic offset"
    return "systematic offset of unknown origin — investigate the control chain"

The doubled values are worth including because a sign error is at least as common as an omission, and it produces exactly twice the offset. A bias of 96 m in a region with a 48 m undulation is unmistakable once the doubling is in the candidate list and baffling when it is not.

Reading a benchmark comparison: three patterns and their causes Three rows. A constant offset across every benchmark, of a metre or more, indicates a vertical datum problem — the wrong geoid, no geoid, or one applied twice. A small constant offset of a few centimetres indicates an antenna or rod height not accounted for, which is arithmetic rather than datum. An offset that varies smoothly across the site indicates a geoid grid of insufficient resolution for the terrain, or a tilted reconstruction, and the two are distinguished by whether the variation follows the geoid's own gradient or the survey's geometry. a constant offset of metres vertical datum — wrong geoid, no geoid, or applied twice a constant offset of centimetres antenna or rod height — arithmetic, not datum an offset varying across the site grid resolution or a tilted reconstruction The size and the spatial pattern of the offset name the cause between them.

Figure 3 — Three patterns, three quite different remedies.

Edge-case matrix

Situation Effect Handling
Benchmark used as control Agrees by construction Hold it out
Benchmark from the same GNSS session Shares the base station error Use a levelled benchmark
Published value derived through the same geoid Tests nothing about the model Prefer levelled heights
Fewer than four benchmarks Bias uncertain Report the uncertainty with the bias
Benchmark disturbed or destroyed One large residual Robust statistics, and inspect it
Benchmarks clustered Bias measured in one place Distribute across the site
Benchmark network has its own tilt Systematic pattern in residuals Fit a trend before concluding
Survey older than a datum revision Legitimate difference Compare the datum definitions

The clustered-benchmarks row deserves a check of its own, because a bias measured at one end of a site says nothing about the other:

import numpy as np


def benchmark_distribution(positions_xy: np.ndarray, site_bounds: tuple) -> dict:
    """Do the benchmarks cover the site, or one corner of it?"""
    left, bottom, right, top = site_bounds
    site_area = max((right - left) * (top - bottom), 1e-9)

    xy = np.asarray(positions_xy, dtype=float)
    hull_extent = (xy[:, 0].max() - xy[:, 0].min()) * (xy[:, 1].max() - xy[:, 1].min())
    coverage = float(hull_extent / site_area)

    return {"benchmarks": len(xy), "coverage_fraction": coverage,
            "adequate": coverage > 0.3 and len(xy) >= 4,
            "note": ("benchmarks span the site" if coverage > 0.3 else
                     "benchmarks are clustered — the bias is measured where they are "
                     "and says little about the rest of the survey")}

Verification snippet

import numpy as np


def check_for_trend(positions_xy: np.ndarray, residuals_m: np.ndarray) -> dict:
    """Is the difference a constant offset, or does it vary across the site?

    A constant is a datum or antenna error. A trend is a tilt, which points
    at the control distribution or a datum realisation difference rather than
    at the geoid model.
    """
    xy = np.asarray(positions_xy, dtype=float)
    d = np.asarray(residuals_m, dtype=float)
    if d.size < 4:
        return {"note": "need at least four benchmarks to separate a trend"}

    centred = xy - xy.mean(axis=0)
    A = np.column_stack([centred, np.ones(len(d))])
    coeffs, *_ = np.linalg.lstsq(A, d, rcond=None)
    gradient = float(np.hypot(coeffs[0], coeffs[1])) * 1000  # mm per metre
    residual = d - A @ coeffs

    return {"offset_m": float(coeffs[2]),
            "tilt_mm_per_100m": gradient * 100,
            "scatter_after_trend_m": float(np.std(residual, ddof=1)),
            "is_tilt": gradient * 100 > 5.0,
            "note": ("a tilt is present — check control distribution and the datum "
                     "realisation" if gradient * 100 > 5 else
                     "the difference is a constant offset")}
Benchmark residuals for four common causes Four residual patterns across a site. A correct survey shows residuals scattered about zero within two centimetres. An omitted geoid conversion shows every residual offset by the local undulation of forty-eight metres. An antenna height applied with the wrong sign shows every residual offset by twice the antenna height, about three point six metres. A tilted control network shows residuals rising linearly from minus four centimetres at one end to plus four at the other. A note states that the pattern names the cause without further investigation. 0 correct — ±2 cm geoid not applied — 48 m antenna sign — 3.6 m tilt — control or datum The pattern names the cause before any investigation. Which is why bias, tilt and scatter are reported separately rather than as one RMSE.

Figure 2 — Four patterns, four causes, one comparison.

Stating the result

A vertical accuracy statement that survives a client’s own check has four parts: the number of benchmarks, the bias with its uncertainty, the scatter, and the datum and model the comparison was against.

The most common weakness is stating an RMSE alone. An RMSE of 5 cm made of a 4 cm bias and 3 cm scatter is a datum problem; the same RMSE made of no bias and 5 cm scatter is an ordinary survey. A client checking against their own data will find the first and not the second, and a statement that distinguished them would have avoided the conversation.

def vertical_accuracy_statement(comparison: dict, trend: dict,
                                vertical_datum: str, geoid_model: str) -> str:
    """The paragraph that belongs in the survey report."""
    parts = [
        f"Heights are referenced to {vertical_datum}, converted using the "
        f"{geoid_model} geoid model.",
        f"Compared against {comparison['benchmarks']} independent levelled benchmarks, "
        f"the survey shows a mean difference of {comparison['bias_m']:+.3f} m "
        f"(±{comparison['bias_uncertainty_m']:.3f} m) "
        f"with a scatter of {comparison['scatter_m']:.3f} m.",
    ]
    if trend.get("is_tilt"):
        parts.append(f"A residual tilt of {trend['tilt_mm_per_100m']:.0f} mm per 100 m "
                     "is present across the site.")
    return " ".join(parts)

When to escalate

  • The bias is significant and matches no candidate cause. Work back through the control chain: base station height, antenna height, the control survey’s own datum, then the conversion.
  • No independent benchmarks exist near the site. A levelling run from the nearest benchmark is the proper answer. Where that is impractical, state that the vertical datum is unverified rather than implying it was checked.
  • The benchmark network itself is suspect. Older networks have their own distortions. Compare against several benchmarks and treat a consistent pattern as informative rather than as the survey’s fault.

← Geoid Models and Vertical Datum Automation