Change Detection Between Survey Epochs

Monitoring is where drone surveying earns repeat work: a quarry face checked monthly, a landfill cell tracked through filling, an embankment watched for movement after heavy rain. In every case the deliverable is not a survey but a difference, and differences are far more sensitive to processing errors than the surveys they come from — because the signal is small and the errors are not.

A two-centimetre disagreement between two surveys is excellent survey work. On a site where the real movement being watched for is five centimetres, it is also half the signal. That ratio is what makes change detection a discipline of its own rather than a subtraction.

This page covers the four things that separate a defensible change map from a noisy one: co-registering the epochs on ground that did not move, deriving a detection limit from the data rather than asserting one, choosing between raster differencing and cloud-to-cloud distances, and reporting the result so a reviewer can check it. It builds on the classified clouds from classifying point clouds with PDAL and Python.

Audience and prerequisites. Python 3.10+, two surveys of the same site in the same projected CRS with the same vertical datum, and some ground that is known not to have changed. That last item is a requirement, not a convenience: without stable ground there is nothing to co-register against and no way to measure the detection limit.

# Prerequisites

# Conceptual architecture

The pipeline has a fixed shape. Both epochs are filtered to comparable content — the same classes, the same extent — then co-registered on stable ground, then differenced, then thresholded against a limit derived from the residual disagreement, then summarised into areas and volumes.

Two decisions inside that shape do the work. The first is what counts as stable: a polygon of hardstanding, roofs and roads that the site confirms has not been touched. The co-registration and the detection limit both come from it, so a poorly chosen stable set corrupts everything downstream. The second is raster or cloud: a DEM of difference is fast, easy to read and measures only vertical change, while a cloud-to-cloud distance measures change in whatever direction it occurred and is the only correct choice on a vertical face.

Figure 1 — Everything hangs off the stable-ground polygon. It is the one input a site visit should establish, not a technician.

# Step 1: Co-register on ground that did not move

The two epochs are independently georeferenced, which means they carry independent errors. Differencing them without reconciling those errors turns a systematic offset into apparent site-wide movement — the single most common way a change map ends up entirely coloured.

import numpy as np
import open3d as o3d


def coregister_on_stable(cloud_a: np.ndarray, cloud_b: np.ndarray,
                         stable_mask_a: np.ndarray, stable_mask_b: np.ndarray,
                         *, max_distance: float = 0.5) -> dict:
    """Rigid-body fit of epoch B onto epoch A, using stable ground only.

    Fitting on the whole site is the classic mistake: real change then drags
    the transform, and the result under-reports exactly the movement being
    looked for. Only stable ground may inform the registration.
    """
    pa = o3d.geometry.PointCloud()
    pa.points = o3d.utility.Vector3dVector(cloud_a[stable_mask_a])
    pb = o3d.geometry.PointCloud()
    pb.points = o3d.utility.Vector3dVector(cloud_b[stable_mask_b])
    pa.estimate_normals()

    result = o3d.pipelines.registration.registration_icp(
        pb, pa, max_distance, np.eye(4),
        o3d.pipelines.registration.TransformationEstimationPointToPlane())

    t = result.transformation
    shift = np.linalg.norm(t[:3, 3])
    rot_deg = np.degrees(np.arccos(np.clip((np.trace(t[:3, :3]) - 1) / 2, -1, 1)))
    return {"transform": t, "fitness": result.fitness,
            "inlier_rmse": result.inlier_rmse,
            "translation_m": float(shift), "rotation_deg": float(rot_deg)}

Inspect the transform before applying it. A translation of a few centimetres is ordinary; one of half a metre means the surveys disagree badly enough that something upstream is wrong, and applying the fit would hide it. A rotation of more than a few hundredths of a degree on a site a few hundred metres across is similarly a red flag rather than a correction. Aligning two epochs with ICP before differencing covers the diagnostics in full.

# Step 2: Derive the detection limit from the residual

After co-registration the two surveys still disagree on stable ground, and that residual disagreement is the noise floor of the whole exercise. Measuring it is what turns a colourful difference raster into a map with a defensible threshold.

import numpy as np


def detection_limit(stable_diff: np.ndarray, *, confidence: float = 1.96) -> dict:
    """Level of change below which the two surveys cannot be distinguished.

    Robust statistics throughout: a median and a scaled median absolute
    deviation, because the stable set almost always contains a few cells
    that were not as stable as the site believed.
    """
    r = stable_diff[np.isfinite(stable_diff)]
    if r.size < 500:
        raise ValueError(f"only {r.size} stable cells — the limit would be unstable")

    bias = float(np.median(r))
    sigma = float(np.median(np.abs(r - bias)) * 1.4826)
    return {"residual_bias_m": bias, "residual_sigma_m": sigma,
            "detection_limit_m": confidence * sigma,
            "n_stable_cells": int(r.size)}

Publishing the limit alongside the map is what makes the product reviewable. “No significant change” means nothing on its own; “no change exceeding ±7 cm, which is the 95 % limit measured over 18,400 stable cells” is a statement a reviewer can accept or challenge.

# Step 3: Difference, and choose the right comparison

A DEM of difference is a subtraction of two aligned rasters. It is fast, trivially mappable, and measures only the vertical component of any change — which on a quarry face or a retaining wall is close to none of it.

import numpy as np
import rasterio


def dem_of_difference(path_a: str, path_b: str, out_path: str,
                      limit: float) -> dict:
    """Vertical difference between two aligned surfaces, thresholded."""
    with rasterio.open(path_a) as a, rasterio.open(path_b) as b:
        if a.transform != b.transform or a.shape != b.shape:
            raise ValueError("rasters are not on the same grid — resample first")
        za = a.read(1, masked=True).filled(np.nan)
        zb = b.read(1, masked=True).filled(np.nan)
        profile = a.profile

    diff = zb - za
    significant = np.where(np.abs(diff) > limit, diff, np.nan)

    profile.update(dtype="float32", nodata=np.nan, count=1)
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(significant.astype("float32"), 1)
        dst.update_tags(1, DETECTION_LIMIT_M=f"{limit:.4f}")

    cell = abs(profile["transform"].a * profile["transform"].e)
    return {"changed_area_m2": float(np.isfinite(significant).sum() * cell),
            "net_volume_m3": float(np.nansum(significant) * cell),
            "gain_m3": float(np.nansum(np.where(significant > 0, significant, 0)) * cell),
            "loss_m3": float(np.nansum(np.where(significant < 0, significant, 0)) * cell)}

Writing the detection limit into the raster’s own tags means the threshold travels with the product. A change raster whose threshold is recorded only in a report is one that will eventually be re-interpreted with the wrong assumption.

Where the geometry is not broadly horizontal, the vertical difference understates the change by the cosine of the surface angle — and on a vertical face it reports zero regardless of how much material moved. The correct tool there is a cloud-to-cloud distance measured along the local surface normal, covered in computing M3C2 distances on drone point clouds.

Figure 2 — Why a DEM of difference is the wrong instrument on a face. It is not imprecise; it is measuring a different quantity.

# Step 4: Summarise into something a client reads

A change raster is an intermediate. What a monitoring report needs is a handful of numbers and a small number of named regions, and producing them automatically stops every month’s report being a manual exercise.

import numpy as np
from scipy import ndimage


def summarise_change(diff: np.ndarray, cell: float, limit: float,
                     min_area_m2: float = 5.0) -> dict:
    """Connected change regions above a minimum area, with per-region totals."""
    signif = np.isfinite(diff) & (np.abs(diff) > limit)
    labelled, n = ndimage.label(signif)

    regions = []
    for i in range(1, n + 1):
        m = labelled == i
        area = float(m.sum() * cell ** 2)
        if area < min_area_m2:
            continue
        vol = float(np.nansum(diff[m]) * cell ** 2)
        ys, xs = np.nonzero(m)
        regions.append({"id": i, "area_m2": area, "volume_m3": vol,
                        "mean_change_m": float(np.nanmean(diff[m])),
                        "centroid_px": (float(ys.mean()), float(xs.mean())),
                        "direction": "gain" if vol > 0 else "loss"})

    regions.sort(key=lambda r: -abs(r["volume_m3"]))
    return {"regions": regions[:20], "region_count": len(regions),
            "total_changed_area_m2": sum(r["area_m2"] for r in regions),
            "net_volume_m3": sum(r["volume_m3"] for r in regions)}

The minimum-area filter is doing real work. Below about five square metres, a “region” is a handful of cells that cleared the threshold by chance, and a report listing four hundred of them buries the three that matter.

# Step 5: Build a monitoring series, not a pair of surveys

A single difference answers “what changed since last time”. A monitoring contract asks a harder question: “is this moving, and how fast”. Answering it needs the pipeline to treat the surveys as a series rather than as independent pairs, and three practices make that work.

Register everything to the first epoch, not to the previous one. Chaining registrations accumulates error: each fit carries a small residual, and after twelve months the twelfth survey is several centimetres from the first through nothing but drift. Registering every epoch directly onto a single reference — usually the first, or a dedicated control survey — keeps every difference comparable to every other.

Track a fixed set of monitoring points. A time series of elevation at named locations is far more useful to an engineer than a series of maps, because it can be plotted, trended and compared against a trigger level. The points cost nothing to define and are sampled from each epoch in the same way.

import numpy as np


def sample_monitoring_points(surface: np.ndarray, transform,
                             points_xy: list[tuple[float, float]],
                             window: int = 3) -> list[float]:
    """Median elevation in a small window at each monitoring point.

    A single-cell sample is noisy; a three-by-three median is stable and
    still local enough to track a real movement of a few centimetres.
    """
    out = []
    inv = ~transform
    h, w = surface.shape
    for x, y in points_xy:
        col, row = (int(v) for v in inv * (x, y))
        r0, r1 = max(row - window // 2, 0), min(row + window // 2 + 1, h)
        c0, c1 = max(col - window // 2, 0), min(col + window // 2 + 1, w)
        patch = surface[r0:r1, c0:c1]
        finite = patch[np.isfinite(patch)]
        out.append(float(np.median(finite)) if finite.size else float("nan"))
    return out

Report rates, with their own uncertainty. A point that has moved 4 cm over six surveys has a trend that can be fitted, and the standard error of that trend is usually far smaller than the detection limit of any single pair — which means a series can detect movement that no individual difference could. This is the strongest argument for monitoring at a regular interval rather than on demand: the statistics improve with every flight.

import numpy as np


def movement_rate(times_days: np.ndarray, elevations_m: np.ndarray) -> dict:
    """Linear trend through a monitoring point's history, with its uncertainty."""
    ok = np.isfinite(elevations_m)
    t, z = times_days[ok], elevations_m[ok]
    if t.size < 4:
        return {"verdict": "too few epochs to fit a trend"}

    A = np.column_stack([t, np.ones_like(t)])
    (slope, intercept), residuals, *_ = np.linalg.lstsq(A, z, rcond=None)
    resid = z - (slope * t + intercept)
    se = float(np.sqrt(np.sum(resid ** 2) / (t.size - 2)
                       / np.sum((t - t.mean()) ** 2)))

    return {"rate_mm_per_year": float(slope * 365.25 * 1000),
            "rate_se_mm_per_year": float(se * 365.25 * 1000),
            "significant": bool(abs(slope) > 1.96 * se),
            "epochs": int(t.size)}

A trend that is significant at the 95 % level after eight flights, each individually unable to resolve the movement, is the product a monitoring contract is really buying. Reporting it alongside the maps changes the conversation from “did anything change this month” to “this is moving at 11 ± 3 mm per year”, which is a statement an engineer can design against.

# Parameter deep-dive

# Verification and output inspection

The strongest check is one the data provides for free: the stable set should show no change after processing.

import numpy as np


def verify_stable_shows_nothing(diff: np.ndarray, stable_mask: np.ndarray,
                                limit: float, tol: float = 0.02) -> dict:
    """Stable ground must come out clean, or the pipeline is reporting noise."""
    s = diff[stable_mask & np.isfinite(diff)]
    exceed = float(np.count_nonzero(np.abs(s) > limit) / max(s.size, 1))
    bias = float(np.median(s))

    problems = []
    if exceed > 0.05:
        problems.append(f"{exceed:.1%} of stable cells exceed the limit — "
                        "the limit is too tight or the registration failed")
    if abs(bias) > tol:
        problems.append(f"stable ground shows a {bias:+.3f} m offset after "
                        "co-registration — check the vertical datum")
    return {"stable_exceed_fraction": exceed, "stable_bias_m": bias,
            "problems": problems}

By construction about five percent of stable cells should exceed a 95 % limit. Substantially more means something is unaccounted for; substantially fewer means the limit is conservative and real change is being suppressed. Both are worth knowing before the map is published.

Figure 3 — Every change figure carries this chain with it.

# Troubleshooting

The whole site shows uniform change. A datum or georeferencing offset, not movement. Check the median difference over stable ground before anything else; co-registration should remove it.

Change appears as paired gain and loss on opposite slopes. A horizontal misalignment between epochs. The pattern is a dipole and sums to near zero. Re-run the registration with a tighter correspondence distance.

The change map is entirely coloured. No detection limit was applied, or it was asserted rather than measured. Derive it from the stable residual.

A vertical face shows no change after a visible rockfall. A DEM of difference measures only the vertical component. Use a cloud-to-cloud distance along the surface normal.

The registration reports excellent fitness and the result is wrong. The stable set probably included ground that moved, so the transform absorbed real change. Shrink the stable set to surfaces the site can confirm.

Change is reported where vegetation grew. Expected, and usually not what the client means by change. Filter both epochs to ground and structures before differencing, using the classes from the classification stage.

# Related

• Computing DEM of Difference Rasters in Python
• Separating Real Change from Noise with Detection Limits
• Computing M3C2 Distances on Drone Point Clouds
• Aligning Two Epochs with ICP Before Differencing
• Reconciling Volume Differences Between Flights

← Point Cloud Processing & 3D Deliverables