Detecting Blunders in Checkpoint Residuals

A checkpoint set of twenty produces nineteen residuals under four centimetres and one at thirty-one. The quoted RMSE is 8.2 cm, which fails the specification. Delete the outlier and it becomes 3.1 cm, which passes comfortably.

Both numbers are defensible and only one is honest, and which one depends entirely on why that point is at thirty-one centimetres. This page sets out how to answer that question with a procedure rather than an instinct, so that the decision can be written down and defended. It belongs to checkpoint-based accuracy validation in Python.

# Blunder versus genuine error

The distinction is about mechanism, not magnitude.

A blunder is a point whose measured value is wrong: a transposed northing, a rod height not subtracted, a target photographed after it was moved, a point observed in the wrong datum, two points whose identifiers were swapped in the field book. The reconstruction is fine; the reference is not. Including it measures the note-taking, not the survey.

A genuine large error is a point where the reconstruction really is that far out: a corner with no control nearby, a patch of water or moving vegetation, a steep face seen from one direction only. Excluding it hides exactly the weakness the checkpoint existed to expose.

The statistics cannot tell these apart. Only evidence can — which is why a blunder test is a way of nominating suspects, never a way of removing them.

Figure 1 — Two mechanisms, one number.

# Minimal reproducible solution

The standard test is the normalised residual against a robust scale. Robust matters: the mean and standard deviation are themselves pulled by the blunder you are trying to find, so a point large enough to matter masks itself.

import numpy as np


def robust_scale(values: np.ndarray) -> float:
    """Median absolute deviation, scaled to be comparable with a standard deviation.

    The 1.4826 factor makes the MAD of normally distributed data equal its
    standard deviation, so the thresholds below keep their usual meaning while
    the estimate itself stays immune to a handful of gross outliers.
    """
    median = np.median(values)
    mad = np.median(np.abs(values - median))
    return float(1.4826 * mad)


def nominate_blunders(residuals: np.ndarray, ids: list[str],
                      *, threshold: float = 3.5) -> list[dict]:
    """Nominate suspect checkpoints. Nominating is not excluding."""
    scale = robust_scale(residuals)
    if scale == 0.0:                       # identical residuals, degenerate set
        return []
    centre = float(np.median(residuals))
    scores = np.abs(residuals - centre) / scale
    return [
        {"id": ids[i], "residual_m": float(residuals[i]),
         "score": float(scores[i]), "threshold": threshold}
        for i in np.argsort(-scores) if scores[i] > threshold
    ]

A threshold of 3.5 on the robust score is the usual starting point. It flags roughly one point in two thousand from clean normal data, which on a twenty-point set means a flag is far more likely to be real than to be chance.

# Why the robust scale is not optional

Substituting the ordinary standard deviation produces a test that fails exactly when it is needed. With nineteen residuals near 3 cm and one at 31 cm, the sample standard deviation is inflated by the outlier itself to around 6.5 cm, so the offender scores about 4.3 — near the threshold, and one more blunder would push it safely under. This is masking, and it is why a rule of “anything over three sigma” quietly misses pairs of blunders.

def compare_scales(residuals: np.ndarray) -> dict:
    """Show how much the outlier inflates a non-robust scale estimate."""
    classical = float(np.std(residuals, ddof=1))
    robust = robust_scale(residuals)
    worst = float(np.max(np.abs(residuals - np.median(residuals))))
    return {
        "classical_sd_m": classical,
        "robust_scale_m": robust,
        "inflation_factor": classical / robust if robust else float("inf"),
        "worst_score_classical": worst / classical if classical else float("inf"),
        "worst_score_robust": worst / robust if robust else float("inf"),
    }

Run this on any set where a nomination sits near the threshold. An inflation factor above about 1.5 says the classical estimate has already been captured by the very points you are testing.

# Investigating a nomination

A nomination is a task, and there are five checks that between them explain almost every real blunder.

Compare against the field record. A transposition shows up immediately when the manifest value is set beside the field book. This finds more blunders than the other four combined.

Check the residual’s direction. A rod height produces a residual that is almost purely vertical and close to a round number — 1.5 m, 2.0 m, or the instrument height. A datum mistake produces a large vertical offset shared by every point, not one.

Check the pattern across neighbours. If the three nearest checkpoints are also elevated, the problem is regional weakness, not a bad point. A single high residual surrounded by clean ones points at the point itself.

Look at the imagery. Open the photographs covering the target. A target that is partly obscured, badly lit, or visibly in a different place between flight lines answers the question outright.

Check for a swap. Two points whose residuals are large and roughly opposite are usually the same two identifiers, exchanged.

def characterise(residual: dict, *, tolerance: float = 0.05) -> list[str]:
    """Cheap hypotheses for a nominated point, to direct the investigation."""
    notes = []
    de, dn, dh = residual["d_east"], residual["d_north"], residual["d_height"]
    horizontal = float(np.hypot(de, dn))

    if abs(dh) > 4 * max(horizontal, 1e-6):
        notes.append("almost purely vertical — check rod or antenna height")
    for common in (1.0, 1.5, 1.8, 2.0):
        if abs(abs(dh) - common) < tolerance:
            notes.append(f"vertical offset near {common} m — likely an instrument height")
    if horizontal > 4 * max(abs(dh), 1e-6):
        notes.append("almost purely horizontal — check identification or transposition")
    if abs(abs(de) - abs(dn)) < tolerance and horizontal > 0.2:
        notes.append("easting and northing offsets similar — possible transposition")
    return notes or ["no characteristic signature — inspect imagery and field record"]

Figure 3 — Why the scale estimate has to be robust.

# Edge-case matrix

The fifth row is the one that matters most. An exclusion without a recorded reason turns an accuracy statement into an assertion, because nothing in the manifest lets a reader check whether the removal was justified.

# Verification snippet

def resolve_nomination(nomination: dict, *, verdict: str, evidence: str) -> dict:
    """Close out a nomination with an explicit, recorded decision.

    Every nomination ends in exactly one of three states. The function refuses
    anything else, which is what stops a survey from quietly shipping with a
    point that was removed because it was inconvenient.
    """
    allowed = {"blunder_excluded", "genuine_retained", "corrected_and_reprocessed"}
    if verdict not in allowed:
        raise ValueError(f"verdict must be one of {sorted(allowed)}")
    if verdict == "blunder_excluded" and len(evidence.strip()) < 20:
        raise ValueError("an exclusion needs recorded evidence, not a label")
    return {
        "checkpoint": nomination["id"],
        "residual_m": nomination["residual_m"],
        "score": nomination["score"],
        "verdict": verdict,
        "evidence": evidence,
    }


def audit_exclusions(resolutions: list[dict], total_checkpoints: int) -> dict:
    """Is the exclusion rate itself defensible?"""
    excluded = [r for r in resolutions if r["verdict"] == "blunder_excluded"]
    rate = len(excluded) / total_checkpoints if total_checkpoints else 0.0
    return {
        "excluded": len(excluded),
        "rate": rate,
        # Above roughly one in ten, the problem is the field procedure rather
        # than a run of bad luck, and the survey should be re-observed.
        "acceptable": rate <= 0.10,
        "ids": [r["checkpoint"] for r in excluded],
    }

Figure 2 — The three endings a nomination is allowed to have.

# When to escalate

• More than one checkpoint in ten is nominated. The field procedure is producing the blunders, and re-observation is cheaper than arguing about which removals are justified.
• Every residual is large. This is not a blunder problem. Check the datum, the geoid model and the scale before touching any individual point.
• A nomination cannot be resolved either way. Retain it, and say in the report that it was retained without explanation. A slightly worse figure that is true is worth more than a good one that is not.

# Related

• Checkpoint-Based Accuracy Validation in Python
• Computing RMSE and Accuracy Statistics from Checkpoints
• Reporting Accuracy in a Machine-Readable Manifest

← Checkpoint-Based Accuracy Validation in Python