Separating Real Change from Noise with Detection Limits

The client looks at the difference map and asks which of the red areas matter. It is a fair question and, as produced, an unanswerable one: the map shows every cell where the two surveys disagreed by any amount, and most of those disagreements are the surveys rather than the site.

The missing piece is a number — the smallest change this pair of surveys could reliably have seen. With it, the map becomes two categories instead of a continuum, and every coloured cell is one the reader can trust. Without it, the reader is being asked to do the statistics by eye.

This page derives that number from the data, extends it to vary across a site where accuracy varies, and covers how to state it so a reviewer can check it. It is the thresholding step of change detection between survey epochs.

Why the limit must come from the data

A sensor datasheet gives a per-point precision under ideal conditions. A survey report gives an RMSE against checkpoints, usually measured on open hardstanding. Neither describes what the pipeline can resolve between two specific flights over one specific site, because the actual limit is set by everything that differed between them: lighting, vegetation state, flight geometry, control quality, processing parameters.

The data already contains the answer. Ground that the site confirms has not changed — a concrete apron, a road, a building roof — should difference to zero. It does not, and the spread of how much it misses by is the noise floor. Everything smaller than a couple of standard deviations of that spread is indistinguishable from the surveys disagreeing.

Two properties make this a good estimator. It is empirical, so it captures causes nobody anticipated. And it is checkable: a reviewer can be handed the stable polygon and reproduce the number.

Three candidate detection limits for the same pair of surveys A histogram of elevation differences over stable ground, with three candidate thresholds drawn on it. The sensor datasheet figure of two centimetres sits well inside the distribution, so using it would report most of the stable ground as change. The checkpoint RMSE of four centimetres sits near the shoulder and would still report about a fifth of stable ground as change. The measured limit of seven point eight centimetres, at ninety-five percent of the stable distribution, sits outside almost all of it. A note records how many stable cells each threshold would falsely flag. datasheet 2 cm would flag 61 % of stable ground checkpoint RMSE 4 cm would flag 21 % measured 7.8 cm flags 5 %, as designed difference over stable ground Only the third threshold was measured on these two flights over this site. The other two describe equipment and conditions, not the comparison being made.

Figure 1 — Three numbers that all sound authoritative. Only one of them is about this pair of surveys.

Minimal reproducible solution

import numpy as np


def detection_limit(diff: np.ndarray, stable_mask: np.ndarray,
                    *, confidence: float = 1.96,
                    min_cells: int = 500) -> dict:
    """Minimum detectable change, measured over ground known not to have moved.

    Robust estimators throughout. A stable polygon drawn from a site walk
    almost always includes a few cells that were not stable — a repainted
    line, a parked vehicle that moved — and a mean and standard deviation
    would let those set the threshold for the whole map.
    """
    r = diff[stable_mask & np.isfinite(diff)]
    if r.size < min_cells:
        raise ValueError(f"{r.size} stable cells is too few for a stable estimate")

    bias = float(np.median(r))
    sigma = float(np.median(np.abs(r - bias)) * 1.4826)
    limit = confidence * sigma

    return {"detection_limit_m": limit,
            "residual_bias_m": bias,
            "residual_sigma_m": sigma,
            "n_stable_cells": int(r.size),
            "confidence": confidence,
            "expected_false_positive_rate": 2 * (1 - _phi(confidence))}


def _phi(z: float) -> float:
    """Standard normal CDF, without pulling in scipy for one value."""
    import math
    return 0.5 * (1 + math.erf(z / math.sqrt(2)))


def apply_limit(diff: np.ndarray, limit: float, bias: float = 0.0) -> np.ndarray:
    """Mask everything the surveys could not distinguish from agreement."""
    centred = diff - bias
    return np.where(np.isfinite(centred) & (np.abs(centred) > limit),
                    centred, np.nan)

Subtracting the residual bias before thresholding is a small step with a large effect. If the two surveys sit 1.5 cm apart on stable ground after registration — which is common and not a failure — then leaving the bias in place makes the threshold asymmetric: change of one sign needs to exceed 6.3 cm and change of the other 9.3 cm. Centring first restores symmetry.

Sweeping the threshold

A single threshold applied once is a choice nobody can check. Sweeping it and plotting the result turns that choice into evidence.

Detected change area against the threshold applied A curve of detected change area against the threshold applied to a difference raster. At a very low threshold almost the whole site registers as changed, because noise dominates. The curve falls steeply, then flattens into a plateau across a broad range of thresholds, then falls again as genuine change starts being excluded. A note states that the plateau is where the answer is insensitive to the exact threshold, that this is the defensible region to report from, and that a threshold chosen on the steep part of the curve gives a number that changes with the threshold rather than with the site. threshold applied to the difference raster detected change area The plateau is where the answer stops depending on the threshold. Report from there.

Figure 3 — Sweep the threshold; the plateau is the defensible region.

Edge-case matrix

Situation Naive handling Correct handling
Stable set includes ground that moved Limit inflated, real change hidden Robust statistics, and re-check the polygon
Stable set on hardstanding only Limit optimistic for vegetated areas Spatially varying limit, or state the restriction
Fewer than a few hundred stable cells Unstable estimate Refuse, and say why
Strongly non-normal residuals ±1.96σ is not 95 % Use empirical percentiles instead
Systematic tilt across the site Single limit misleads De-trend first, or vary the limit spatially
Change smaller than the limit but real Reported as no change Correct, and the report should say what was resolvable
Two surveys from the same processing run Limit implausibly small Flag as a likely duplicate input
Seasonal vegetation difference Huge apparent change Restrict to ground and structures before differencing

The non-normal row deserves a note. Residuals from photogrammetric surfaces frequently have heavier tails than a Gaussian, so the 1.96σ figure under-covers. Where that matters, take the empirical 2.5th and 97.5th percentiles of the stable residual and use those directly as the limits — the arithmetic is simpler and makes no distributional assumption.

import numpy as np


def empirical_limits(diff: np.ndarray, stable_mask: np.ndarray,
                     lower_pct: float = 2.5, upper_pct: float = 97.5) -> dict:
    """Asymmetric limits straight from the stable distribution, no assumptions."""
    r = diff[stable_mask & np.isfinite(diff)]
    lo, hi = np.percentile(r, [lower_pct, upper_pct])
    return {"lower_m": float(lo), "upper_m": float(hi),
            "width_m": float(hi - lo), "n": int(r.size)}

A limit that varies across the site

A single number assumes the survey is equally good everywhere, and it never is. Accuracy degrades with distance from control, in areas of poor texture, on steep slopes and under partial vegetation. On a large site the difference between the best and worst areas is routinely a factor of three, which means one global limit is simultaneously too strict in the good areas and too loose in the bad.

The practical approach is to model the limit as a function of an available covariate — most usefully the local point density, which is high where reconstruction went well and low where it did not.

import numpy as np


def spatial_detection_limit(diff: np.ndarray, stable_mask: np.ndarray,
                            density: np.ndarray, *, bins: int = 5,
                            confidence: float = 1.96) -> np.ndarray:
    """Per-cell detection limit, binned by reconstruction point density.

    Cells are grouped into density quantiles; each group's own stable
    residual sets the limit for every cell in that group. Sparse areas
    correctly get a looser threshold than dense ones.
    """
    valid = np.isfinite(diff) & np.isfinite(density)
    edges = np.nanpercentile(density[valid], np.linspace(0, 100, bins + 1))
    out = np.full(diff.shape, np.nan, dtype="float32")

    for i in range(bins):
        in_bin = valid & (density >= edges[i]) & (density <= edges[i + 1])
        r = diff[in_bin & stable_mask]
        if r.size < 200:
            continue                      # not enough stable cells in this bin
        bias = np.median(r)
        sigma = np.median(np.abs(r - bias)) * 1.4826
        out[in_bin] = confidence * sigma
    # Fall back to the global limit wherever a bin had too little stable ground.
    global_r = diff[stable_mask & np.isfinite(diff)]
    gb = np.median(global_r)
    out[np.isnan(out) & valid] = confidence * np.median(np.abs(global_r - gb)) * 1.4826
    return out

The payoff is concrete: on a site where the yard reconstructs at 400 points per square metre and the vegetated margin at 30, the yard gets a 4 cm limit and the margin a 12 cm one. Real 6 cm movement in the yard is then reported, where a global 12 cm limit would have suppressed it; and noise in the margin is not reported, where a global 4 cm limit would have filled it with false change.

Global against density-binned detection limits on one site Five density bins across a site, from thirty points per square metre in the vegetated margin to four hundred in the yard. A horizontal line marks a single global limit of eight centimetres. The per-bin limits fall from twelve centimetres in the sparsest bin to four centimetres in the densest. Two annotations mark the consequences of the global line: real six centimetre movement in the dense yard is suppressed, and noise in the sparse margin is reported as change. 30 80 150 260 400 pts/m² reconstruction point density detection limit global 8 cm 12 cm 4 cm global line reports noise here and hides real change here One number for a site whose accuracy varies by a factor of three is wrong in both directions at once.

Figure 2 — The cost of a single threshold, in the two directions it fails.

Verification snippet

import numpy as np


def verify_limit(diff: np.ndarray, stable_mask: np.ndarray,
                 limit: float, expected_rate: float = 0.05) -> dict:
    """The limit must flag stable ground at close to the design rate."""
    s = diff[stable_mask & np.isfinite(diff)]
    rate = float(np.count_nonzero(np.abs(s) > limit) / max(s.size, 1))
    ratio = rate / expected_rate if expected_rate else float("inf")
    verdict = ("as designed" if 0.5 < ratio < 2.0 else
               "limit too tight — stable ground is being reported as change"
               if ratio >= 2.0 else
               "limit too loose — real change may be suppressed")
    return {"observed_false_positive_rate": rate, "expected": expected_rate,
            "verdict": verdict}

When to escalate

  • No stable ground exists. On a site under complete reconstruction there is nothing to measure against, and the honest answer is that the comparison has no measurable detection limit. Establish a control area — even a small concrete pad — for future flights.
  • The limit exceeds the change the client cares about. The survey cannot answer their question. Say so before delivering a map; more flights at the same quality will not help, but better control or a lower flying height might.
  • Stable ground shows a spatial pattern rather than scatter. A tilt or a dome is a georeferencing problem, and thresholding around it will systematically mis-report one half of the site. Fix the survey before setting any limit.

← Change Detection Between Survey Epochs