Checkpoint-Based Accuracy Validation in Python

A survey report states a horizontal accuracy of 1.8 cm, derived from the ground control residuals. Those residuals measure how closely the bundle adjustment fitted the points it was told to fit — which, with enough free parameters, can be made arbitrarily small without the survey becoming any more accurate.

The honest measure comes from points the adjustment never saw. Holding a subset of the control out as checkpoints costs nothing in the field and converts a self-assessment into a measurement. This page covers the split, the statistics that follow, and the reporting that makes them defensible. It builds on the placement work in ground control point optimization and coordinate sync.

Audience and prerequisites. Python 3.10+, a survey with more control points than the minimum, and a reconstruction that can be run with a subset of them.

# Prerequisites

# Conceptual architecture

Control points and checkpoints are the same measurements used for different purposes. A control point enters the adjustment as a constraint; a checkpoint is withheld and compared against the result afterwards.

That difference makes their residuals mean opposite things. A small control residual says the solver fitted its constraints, which it will do whether or not the survey is correct. A small checkpoint residual says the survey predicts a measurement it never saw, which is the definition of accuracy.

The corollary is that a survey with all its points used as control has no accuracy measurement at all — only a consistency measurement, which is what the report usually quotes.

Figure 1 — Two residuals with opposite meanings, routinely reported as the same thing.

# Step 1: Split without biasing either set

A random split is the obvious choice and is wrong for a small number of points: it can leave the control clustered in one part of the site, which weakens the adjustment, or the checkpoints clustered, which measures accuracy in one place only.

import numpy as np


def stratified_split(positions_xy: np.ndarray, *, checkpoint_fraction: float = 0.3,
                     seed: int = 0) -> dict:
    """Split control and checkpoints so both sets span the site.

    Points are grouped into spatial cells and the split is made within each
    cell, so neither set can end up clustered — which a random split over a
    dozen points does surprisingly often.
    """
    xy = np.asarray(positions_xy, dtype=float)
    rng = np.random.default_rng(seed)

    cells = max(int(np.sqrt(len(xy) / 3)), 1)
    x_edges = np.quantile(xy[:, 0], np.linspace(0, 1, cells + 1))
    y_edges = np.quantile(xy[:, 1], np.linspace(0, 1, cells + 1))
    ix = np.clip(np.digitize(xy[:, 0], x_edges[1:-1]), 0, cells - 1)
    iy = np.clip(np.digitize(xy[:, 1], y_edges[1:-1]), 0, cells - 1)
    key = ix * cells + iy

    control, checkpoints = [], []
    for cell in np.unique(key):
        members = np.flatnonzero(key == cell)
        rng.shuffle(members)
        n_check = max(int(round(len(members) * checkpoint_fraction)), 0)
        if len(members) - n_check < 1:
            n_check = max(len(members) - 1, 0)
        checkpoints.extend(members[:n_check].tolist())
        control.extend(members[n_check:].tolist())

    return {"control": sorted(control), "checkpoints": sorted(checkpoints),
            "control_count": len(control), "checkpoint_count": len(checkpoints)}

Guaranteeing at least one control point per occupied cell is what keeps the adjustment constrained everywhere. A split that leaves a corner of the site with checkpoints and no control measures the accuracy of an unconstrained extrapolation, which is informative but not what was intended.

# Step 2: Compute the statistics that mean something

import numpy as np


def checkpoint_statistics(residuals_xyz: np.ndarray) -> dict:
    """Horizontal and vertical accuracy statistics from checkpoint residuals.

    Bias and scatter are separated throughout, because a systematic offset
    and random scatter have different causes and different remedies, and the
    conventional RMSE conflates them.
    """
    r = np.asarray(residuals_xyz, dtype=float)
    r = r[np.isfinite(r).all(axis=1)]
    if len(r) < 3:
        return {"note": "at least three checkpoints are needed"}

    horizontal = np.linalg.norm(r[:, :2], axis=1)
    vertical = r[:, 2]

    return {
        "checkpoints": int(len(r)),
        "horizontal": {
            "rmse_m": float(np.sqrt(np.mean(horizontal ** 2))),
            "mean_m": float(np.mean(horizontal)),
            "max_m": float(np.max(horizontal)),
            "bias_east_m": float(np.mean(r[:, 0])),
            "bias_north_m": float(np.mean(r[:, 1])),
        },
        "vertical": {
            "rmse_m": float(np.sqrt(np.mean(vertical ** 2))),
            "bias_m": float(np.mean(vertical)),
            "scatter_m": float(np.std(vertical, ddof=1)),
            "max_abs_m": float(np.max(np.abs(vertical))),
        },
        "systematic_fraction": float(
            np.linalg.norm([np.mean(r[:, 0]), np.mean(r[:, 1]), np.mean(r[:, 2])])
            / max(np.sqrt(np.mean((r ** 2).sum(axis=1))), 1e-9)),
    }

The systematic_fraction is the summary worth watching across a programme: it is the share of the total error that is a common offset rather than scatter. A value near zero is an ordinary survey; a value near one means almost everything wrong with the survey is a single shift, which is a datum or control problem rather than a photogrammetric one.

# Step 3: Detect blunders before computing anything

One mis-identified checkpoint drags every statistic, and it is detectable because a blunder is different in kind from measurement noise.

import numpy as np


def flag_blunders(residuals_xyz: np.ndarray, *, threshold: float = 3.5) -> dict:
    """Identify checkpoints whose residual is inconsistent with the rest.

    A modified z-score on the robust scale, rather than a standard deviation,
    because the blunder being detected would inflate a standard deviation
    enough to hide itself.
    """
    magnitude = np.linalg.norm(np.asarray(residuals_xyz, dtype=float), axis=1)
    median = float(np.median(magnitude))
    mad = float(np.median(np.abs(magnitude - median)))
    if mad < 1e-9:
        return {"blunders": [], "note": "residuals are identical; nothing to flag"}

    score = 0.6745 * (magnitude - median) / mad
    blunders = np.flatnonzero(score > threshold).tolist()
    return {"blunders": blunders, "scores": score.tolist(),
            "clean_indices": [i for i in range(len(magnitude)) if i not in blunders],
            "note": (f"{len(blunders)} checkpoint(s) inconsistent with the rest — "
                     "inspect before excluding" if blunders else "no blunders detected")}

Flagging rather than automatically excluding is deliberate. A checkpoint with a large residual may be a mis-identified target, a disturbed marker, or the one place where the survey genuinely is wrong — and only inspection distinguishes them.

# Step 4: Decide how many points a survey needs

The split above assumes there are enough points to divide. Deciding how many to collect is a field-planning question with a defensible answer, and it is driven by what the statistics need rather than by convention.

The adjustment needs enough control to constrain the solution: a minimum of three for a rigid placement, and in practice five or more distributed across the interior as well as the perimeter to suppress the deformation modes described in the camera calibration work. The accuracy statement needs enough checkpoints that its own uncertainty is small relative to the tolerance being claimed — which, for a normal distribution, means the standard error of the RMSE falls roughly as one over the square root of twice the checkpoint count.

import numpy as np


def required_points(tolerance_m: float, expected_accuracy_m: float,
                    *, confidence: float = 0.95) -> dict:
    """How many checkpoints make an accuracy claim statistically meaningful.

    The claim is that the true accuracy is below the tolerance. The number of
    checkpoints needed grows as the expected accuracy approaches the
    tolerance, which is why a survey claiming to just meet a specification
    needs far more validation than one comfortably inside it.
    """
    margin = max(tolerance_m - expected_accuracy_m, 1e-6) / expected_accuracy_m
    z = 1.96 if confidence >= 0.95 else 1.64
    needed = int(np.ceil(0.5 * (z / margin) ** 2)) + 1

    return {"checkpoints_needed": max(needed, 5),
            "control_needed": max(5, int(np.ceil(needed * 0.7 / 0.3 * 0.3))),
            "margin_fraction": margin,
            "note": ("comfortable margin; five checkpoints suffice" if margin > 0.8 else
                     "the claimed tolerance is close to the expected accuracy — "
                     "more checkpoints are needed to support it")}

The result is often uncomfortable. A survey expecting 3 cm accuracy and claiming a 3.5 cm tolerance needs dozens of checkpoints to support the claim at 95 % confidence; the same survey claiming 6 cm needs five. That is not a statistical technicality — it is the reason a tight tolerance costs field time, and it is better known before the flight.

# Step 5: Report the result so it cannot be misread

Three conventions make an accuracy statement robust against the ways it is usually misread.

State the confidence level. An RMSE is a one-sigma figure covering about 68 % of outcomes, and readers routinely treat it as a bound. Where a specification says “within X”, it almost always means 95 %, which is roughly two RMSE for a well-behaved distribution.

State the checkpoint count. An RMSE from five points and one from fifty are different claims, and only the count distinguishes them.

State what was held out. “Accuracy assessed against eight checkpoints withheld from the adjustment” is a different sentence from “control residual RMSE”, and clients who have been given the second when they wanted the first are the reason this page exists.

def accuracy_statement(stats: dict, *, confidence: float = 0.95) -> str:
    """The paragraph that belongs in a survey report."""
    k = 1.96 if confidence >= 0.95 else 1.0
    h = stats["horizontal"]["rmse_m"]
    v = stats["vertical"]["rmse_m"]
    return (
        f"Accuracy was assessed against {stats['checkpoints']} checkpoints withheld "
        f"from the adjustment. Horizontal RMSE is {h:.3f} m and vertical RMSE is "
        f"{v:.3f} m, corresponding to approximately {k * h:.3f} m horizontal and "
        f"{k * v:.3f} m vertical at {confidence:.0%} confidence. "
        f"Vertical bias is {stats['vertical']['bias_m']:+.3f} m."
    )

Generating the sentence from the computed statistics rather than typing it keeps it true, which is the same argument made for every other generated statement in this pipeline: a number written by hand is a claim, and one written by code is a record.

# Step 6: Track accuracy across a programme

A single survey’s accuracy is a number; a programme’s accuracy history is a diagnostic. Three patterns are worth watching over a season of flights.

A gradual rise in checkpoint RMSE usually means something in the workflow has drifted — a camera that needs re-calibrating, control markers that have degraded, or a processing setting that changed. It is invisible in any single survey and obvious in a series.

A step change points at a specific event: a firmware update, a new camera, a changed geoid model, a different operator. Correlating the step against the change log usually names it immediately.

A rising vertical bias with stable scatter is the signature of a datum problem creeping in, and it is the one most worth alerting on because it is the largest error available and the least visible.

import numpy as np


def accuracy_trend(history: list[dict]) -> dict:
    """Trend in checkpoint accuracy across a programme's surveys."""
    if len(history) < 4:
        return {"note": "too few surveys to establish a trend"}

    t = np.arange(len(history), dtype=float)
    h = np.array([s["horizontal"]["rmse_m"] for s in history])
    bias = np.array([s["vertical"]["bias_m"] for s in history])

    h_slope = float(np.polyfit(t, h, 1)[0])
    bias_slope = float(np.polyfit(t, bias, 1)[0])
    steps = np.abs(np.diff(h))
    step_index = int(np.argmax(steps)) if steps.size else -1

    findings = []
    if h_slope > 0.002:
        findings.append(f"horizontal RMSE rising by {h_slope * 1000:.1f} mm per survey")
    if abs(bias_slope) > 0.002:
        findings.append(f"vertical bias drifting by {bias_slope * 1000:+.1f} mm per survey "
                        "— check the vertical datum chain")
    if steps.size and steps.max() > 3 * np.median(steps):
        findings.append(f"a step change at survey {step_index + 1} — correlate against "
                        "the change log")

    return {"surveys": len(history), "horizontal_trend_mm_per_survey": h_slope * 1000,
            "bias_trend_mm_per_survey": bias_slope * 1000, "findings": findings}

# The objection, and the answer to it

The standard objection to holding points out is that it wastes them: a point used as a checkpoint is a point not constraining the adjustment, so the survey is less accurate than it could have been.

That is true and it is the wrong trade to optimise. A survey constrained by ten control points and validated by none has an unknown accuracy; one constrained by seven and validated by three has a slightly worse accuracy that is known. For any deliverable that carries an accuracy statement, the second is worth more, because the statement on the first is not supported by anything.

There is also a practical answer for projects where every point genuinely matters: run the adjustment twice. Once with the split, to measure the accuracy, and once with every point as control, to produce the deliverable. The measured accuracy is then a slight overestimate of the delivered product’s error, which is the conservative direction and is exactly what an accuracy statement should be.

The cost is one extra adjustment run, which on any modern pipeline is minutes.

# Parameter deep-dive

# Verification and output inspection

import numpy as np


def compare_control_and_checkpoints(control_residuals: np.ndarray,
                                    checkpoint_residuals: np.ndarray) -> dict:
    """The ratio that reveals whether a solution is fitting or measuring."""
    c = np.linalg.norm(np.asarray(control_residuals), axis=1)
    k = np.linalg.norm(np.asarray(checkpoint_residuals), axis=1)
    if c.size < 3 or k.size < 3:
        return {"note": "need at least three of each"}

    rms_c = float(np.sqrt(np.mean(c ** 2)))
    rms_k = float(np.sqrt(np.mean(k ** 2)))
    ratio = rms_k / max(rms_c, 1e-9)
    return {"control_rmse_m": rms_c, "checkpoint_rmse_m": rms_k, "ratio": ratio,
            "verdict": ("checkpoints and control agree — the accuracy statement is sound"
                        if ratio < 2 else
                        "checkpoints are much worse than control — the solution is "
                        "fitting its constraints rather than the site")}

A ratio near one is the healthy state and is what a well-constrained survey produces. A ratio above about three is the doming signature described in diagnosing doming and bowl effect in reconstructions, and it means the quoted accuracy is not the survey’s.

Figure 3 — Two numbers that look alike and answer different questions.

Figure 4 — Five stages, one of which cannot be reordered.

# Troubleshooting

Control residuals are excellent and checkpoints are poor. The solution absorbed error into its parameters. Check the intrinsics policy and the control distribution.

One checkpoint is far worse than the rest. A blunder. Inspect the marker and the identification before excluding it.

Checkpoint accuracy varies across the site. Control is unevenly distributed. The poor area is the one furthest from a control point.

Vertical accuracy is much worse than horizontal. Normal to a degree — vertical is typically one and a half to two times horizontal — but a factor of five suggests a datum or calibration problem.

The statistics change a lot with the split. Too few points. Report the uncertainty, or collect more control.

Accuracy is reported and no checkpoints were held out. The figure is a consistency measure. Say so, or hold points out and re-run.

# Related

• Computing RMSE and Accuracy Statistics from Checkpoints
• Splitting Control and Checkpoints Without Bias
• Reporting Accuracy in a Machine-Readable Manifest
• Detecting Blunders in Checkpoint Residuals
• Setting Accuracy Thresholds for Survey Projects

← Ground Control Point Optimization & Coordinate Sync