Computing RMSE and Accuracy Statistics from Checkpoints

A set of checkpoint residuals reduces to a handful of numbers, and which numbers are chosen determines whether the accuracy statement is useful or misleading. An RMSE alone hides a systematic offset. A mean hides the scatter. A combined three-dimensional figure hides the fact that vertical accuracy is always worse than horizontal and is usually what a client cares about.

This page covers computing the statistics that survive scrutiny, from residuals produced by the split described in checkpoint-based accuracy validation in Python.

Why horizontal and vertical stay separate

Photogrammetric error is not isotropic. Vertical accuracy is typically one and a half to two times worse than horizontal, for a structural reason: depth is recovered from the intersection of rays whose convergence angle is limited by the overlap, while horizontal position comes from the image plane directly.

Combining the three axes into one figure therefore reports a number that is dominated by the vertical and is compared against specifications that are usually stated separately. It also hides the diagnostic value of the ratio: a vertical-to-horizontal ratio near two is normal, and a ratio of five says something specific about the geometry or the datum.

Vertical-to-horizontal accuracy ratio and what each value indicates A scale of the ratio between vertical and horizontal checkpoint accuracy. Below one point two is unusually good and suggests the vertical was over-constrained by control. Between one point five and two is the normal range for photogrammetry, set by the ray convergence angle the overlap allows. Between two and three suggests weak vertical geometry, typically from insufficient overlap or a flat flight. Above three indicates a datum or calibration problem rather than a geometric one. A note states that a combined three-dimensional figure hides all of this. < 1.2 over-constrained 1.5 – 2.0 normal for photogrammetry 2.0 – 3.0 weak vertical geometry > 3.0 datum or calibration vertical RMSE ÷ horizontal RMSE The ratio is a diagnostic, and a combined figure discards it. Specifications are stated separately, so the statistics should be too. Vertical is worse because depth comes from ray convergence, which the overlap limits.

Figure 1 — The ratio, and what each range says about the survey.

Minimal reproducible solution

import numpy as np


def accuracy_statistics(residuals_xyz: np.ndarray, *, confidence: float = 0.95) -> dict:
    """Horizontal and vertical accuracy statistics with a stated confidence.

    Bias is reported alongside RMSE throughout. An RMSE of 5 cm made of a
    4 cm bias and 3 cm scatter is a systematic problem, and the same RMSE
    made of pure scatter is an ordinary survey — and only the bias
    distinguishes them.
    """
    r = np.asarray(residuals_xyz, dtype=float)
    r = r[np.isfinite(r).all(axis=1)]
    n = len(r)
    if n < 3:
        return {"note": "at least three checkpoints are needed"}

    horizontal = np.linalg.norm(r[:, :2], axis=1)
    vertical = r[:, 2]
    k = 1.96 if confidence >= 0.95 else 1.0

    def summarise(values, signed=None):
        rmse = float(np.sqrt(np.mean(values ** 2)))
        out = {"rmse_m": rmse, "max_m": float(np.max(np.abs(values))),
               "at_confidence_m": k * rmse}
        if signed is not None:
            out["bias_m"] = float(np.mean(signed))
            out["scatter_m"] = float(np.std(signed, ddof=1))
            out["bias_standard_error_m"] = out["scatter_m"] / np.sqrt(n)
        return out

    h = summarise(horizontal)
    h["bias_east_m"] = float(np.mean(r[:, 0]))
    h["bias_north_m"] = float(np.mean(r[:, 1]))

    return {"checkpoints": n, "confidence": confidence,
            "horizontal": h, "vertical": summarise(np.abs(vertical), vertical),
            "vertical_horizontal_ratio": (summarise(np.abs(vertical))["rmse_m"]
                                          / max(h["rmse_m"], 1e-9))}

What a small sample can support

Accuracy statistics from a handful of points carry substantial uncertainty of their own, and reporting a figure without it invites a client’s own check to disagree for purely statistical reasons.

import numpy as np
from scipy import stats


def rmse_confidence_interval(residuals: np.ndarray,
                             *, confidence: float = 0.95) -> dict:
    """Confidence interval on an RMSE, from the chi-squared distribution.

    An RMSE from five checkpoints has an interval spanning roughly a factor
    of two, which is worth knowing before quoting it to three decimal places.
    """
    r = np.asarray(residuals, dtype=float)
    r = r[np.isfinite(r)]
    n = len(r)
    if n < 3:
        return {"note": "too few points for an interval"}

    rmse = float(np.sqrt(np.mean(r ** 2)))
    alpha = 1 - confidence
    lower = rmse * np.sqrt(n / stats.chi2.ppf(1 - alpha / 2, n))
    upper = rmse * np.sqrt(n / stats.chi2.ppf(alpha / 2, n))

    return {"rmse_m": rmse, "lower_m": float(lower), "upper_m": float(upper),
            "n": n, "relative_width": float((upper - lower) / max(rmse, 1e-9)),
            "note": (f"with {n} checkpoints the RMSE is known to within a factor of "
                     f"{upper / max(lower, 1e-9):.1f}")}

Reporting the interval, or at least the checkpoint count, is what stops a client’s ten-point check disagreeing with a five-point statement and both being right.

Separating the parts a single RMSE figure conflates Three rows. Bias is the mean of the residuals, and a non-zero vertical bias is almost always a datum or geoid problem rather than a precision problem, which makes it the most diagnostic of the three and the one most often invisible inside an RMSE. Scatter is the standard deviation about that mean, and is what most people mean when they say accuracy. RMSE combines the two, which is correct for a specification but useless for diagnosis, because a survey with three centimetres of bias and one of scatter reports the same figure as one with the reverse. bias the mean residual — a non-zero vertical bias is a datum problem, not a precision one scatter the standard deviation about that mean — what most people mean by accuracy RMSE both combined: right for a specification, useless for diagnosis Report all three. The two that RMSE hides are the two that say what to fix.

Figure 3 — One number, two causes, and no way to tell them apart from the number.

Edge-case matrix

Situation Effect Handling
RMSE quoted alone Bias hidden Report bias and scatter separately
Axes combined Vertical dominates, ratio lost Keep horizontal and vertical apart
Confidence unstated One-sigma read as a bound State the level and the multiplier
Fewer than five checkpoints Wide interval Report the count and the interval
A blunder included All statistics inflated Flag and inspect before computing
Residual sign convention unstated Bias sign ambiguous Define survey minus reference
Checkpoints clustered Accuracy of one area Report the spatial coverage
Non-normal residuals Multiplier inexact Use empirical percentiles instead

The sign-convention row matters more than it looks. A vertical bias of +4 cm means different things depending on whether the residual is survey minus reference or the reverse, and a report that does not say has forced the reader to guess which direction to adjust.

Verification snippet

import numpy as np
from scipy import stats


def residual_normality(residuals: np.ndarray) -> dict:
    """Are the residuals consistent with a normal distribution?

    The confidence multipliers assume normality. Heavy tails or skew mean the
    multiplier understates the interval, and empirical percentiles should be
    quoted instead.
    """
    r = np.asarray(residuals, dtype=float)
    r = r[np.isfinite(r)]
    if r.size < 8:
        return {"note": "too few points to assess the distribution"}

    statistic, p_value = stats.shapiro(r)
    skew = float(stats.skew(r))
    kurtosis = float(stats.kurtosis(r))

    return {"n": int(r.size), "p_value": float(p_value),
            "skew": skew, "excess_kurtosis": kurtosis,
            "normal": p_value > 0.05 and abs(skew) < 1.0,
            "note": ("residuals are consistent with normal; the multiplier applies"
                     if p_value > 0.05 else
                     "residuals are not normal — quote empirical percentiles rather "
                     "than a multiplied RMSE")}
Confidence interval on an RMSE against the checkpoint count A band showing the ninety-five percent confidence interval on a measured RMSE against the number of checkpoints, from five to fifty. At five checkpoints the interval spans from about zero point six to two point two times the measured value. At ten it spans zero point seven to one point six. At twenty it spans zero point seven six to one point four. At fifty it spans zero point eight four to one point two. A note states that a five-point RMSE quoted to three decimal places implies a precision the sample cannot support. 1× 510 203550 checkpoints 0.6× to 2.2× at five checkpoints 0.84× to 1.2× at fifty A five-point RMSE quoted to three decimal places implies precision the sample cannot support.

Figure 2 — The interval the sample size actually supports.

Reporting in a form clients can check

The most useful accuracy report is one whose numbers a client can reproduce from the same residuals. That means publishing the residuals, not only the summary.

A table of checkpoint identifier, surveyed coordinate, survey coordinate and the three residual components lets anybody recompute every statistic on this page. It also makes a disagreement productive: instead of two RMSE figures that differ, there are two computations over the same data and the difference is traceable to a method rather than to a mystery.

import pandas as pd


def residual_table(ids: list[str], reference_xyz, survey_xyz) -> pd.DataFrame:
    """The table that should accompany every accuracy statement."""
    import numpy as np
    ref = np.asarray(reference_xyz, dtype=float)
    sur = np.asarray(survey_xyz, dtype=float)
    d = sur - ref
    return pd.DataFrame({
        "checkpoint": ids,
        "ref_e": ref[:, 0], "ref_n": ref[:, 1], "ref_h": ref[:, 2],
        "survey_e": sur[:, 0], "survey_n": sur[:, 1], "survey_h": sur[:, 2],
        "d_east": d[:, 0], "d_north": d[:, 1], "d_height": d[:, 2],
        "d_horizontal": np.linalg.norm(d[:, :2], axis=1),
    })

Matching the statistics to a specification

Accuracy specifications come in several forms and they are not interchangeable. Converting between them correctly is the last step before a claim can be made.

An RMSE specification compares directly against the computed RMSE with no multiplier. A 95 % specification for horizontal accuracy conventionally uses a multiplier of about 1.73 on the horizontal RMSE — not 1.96, because horizontal error combines two axes and follows a different distribution. A 95 % vertical specification uses 1.96 on the vertical RMSE, because that is a single axis. And a maximum error specification is not a statistical statement at all: it is a bound that any single checkpoint can violate.

def against_specification(stats: dict, spec: dict) -> dict:
    """Compare computed statistics against a specification in its own terms."""
    h_rmse = stats["horizontal"]["rmse_m"]
    v_rmse = stats["vertical"]["rmse_m"]
    results = {}

    if "horizontal_rmse_m" in spec:
        results["horizontal_rmse"] = h_rmse <= spec["horizontal_rmse_m"]
    if "horizontal_95_m" in spec:
        results["horizontal_95"] = 1.7308 * h_rmse <= spec["horizontal_95_m"]
    if "vertical_rmse_m" in spec:
        results["vertical_rmse"] = v_rmse <= spec["vertical_rmse_m"]
    if "vertical_95_m" in spec:
        results["vertical_95"] = 1.96 * v_rmse <= spec["vertical_95_m"]
    if "max_error_m" in spec:
        results["max_error"] = (stats["horizontal"]["max_m"] <= spec["max_error_m"]
                                and stats["vertical"]["max_m"] <= spec["max_error_m"])

    return {"results": results, "meets_all": all(results.values()) if results else None}

Using 1.73 rather than 1.96 for the horizontal is the detail most often got wrong, and it matters: applying the single-axis multiplier to a two-axis quantity overstates the error by about thirteen percent, which is enough to fail a survey that met its specification.

When to escalate

  • The vertical-to-horizontal ratio exceeds three. That is not geometry. Check the vertical datum and the camera calibration before adjusting anything about the control.
  • The interval is wider than the tolerance. The claim cannot be supported by this sample. Collect more checkpoints or claim a looser tolerance.
  • Residuals are strongly non-normal. Quote empirical percentiles, and look for a blunder or a spatial pattern rather than treating it as noise.

← Checkpoint-Based Accuracy Validation in Python