Computing DEM of Difference Rasters in Python

The difference raster came out covered in noise, with a faint diagonal corduroy pattern across the whole site and a rim of extreme values around every edge. The two input surfaces are both clean. Nothing is wrong with the survey; what is wrong is the subtraction, and specifically three things about it that a one-line b - a does not handle.

This page covers producing a difference raster that is worth looking at: putting both surfaces on genuinely the same grid, handling NoData so it propagates rather than poisons, attaching the uncertainty that makes the values interpretable, and writing the result so the threshold and provenance travel with the file. It is the raster half of change detection between survey epochs.

# The three faults hiding in b - a

Grid mismatch. Two surveys gridded independently rarely share an origin. Their transforms differ by a fraction of a cell, and resampling one onto the other to make the subtraction possible introduces an interpolation error that is largest exactly where the terrain is steepest. The corduroy pattern is the resampling kernel beating against the grid offset.

NoData poisoning. If either surface uses a sentinel such as −9999 and the subtraction is done on raw arrays, every cell where one survey has no data becomes a difference of thousands of metres. Those cells then dominate every statistic computed afterward, including the ones used to choose a colour ramp, which is why the map often renders as two colours.

Edge effects. The two surveys almost never cover exactly the same extent. Where one has data and the other does not, the difference is undefined — but a naive subtraction over a union extent produces a rim of nonsense, and over an intersection extent silently drops real coverage.

Figure 1 — What a one-line subtraction quietly gets wrong.

# Minimal reproducible solution

Regrid both surfaces onto one explicitly defined grid, mask properly, and subtract.

import numpy as np
import rasterio
from rasterio.warp import reproject, Resampling


def common_grid_difference(path_a: str, path_b: str, out_path: str,
                           *, resampling=Resampling.bilinear) -> dict:
    """Difference two surfaces on a shared grid, with honest NoData.

    Epoch A's grid is adopted as the reference rather than inventing a third,
    because it keeps the result comparable with any earlier difference that
    used the same reference.
    """
    with rasterio.open(path_a) as a, rasterio.open(path_b) as b:
        if a.crs != b.crs:
            raise ValueError(f"CRS mismatch: {a.crs} vs {b.crs}")
        za = a.read(1, masked=True).filled(np.nan).astype("float32")

        zb = np.full(za.shape, np.nan, dtype="float32")
        reproject(source=rasterio.band(b, 1), destination=zb,
                  src_transform=b.transform, src_crs=b.crs,
                  dst_transform=a.transform, dst_crs=a.crs,
                  src_nodata=b.nodata, dst_nodata=np.nan,
                  resampling=resampling)
        profile = a.profile

    both = np.isfinite(za) & np.isfinite(zb)
    diff = np.where(both, zb - za, np.nan).astype("float32")

    profile.update(dtype="float32", nodata=np.nan, count=1, compress="deflate")
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(diff, 1)
        dst.update_tags(1,
                        EPOCH_A=path_a, EPOCH_B=path_b,
                        RESAMPLING=resampling.name,
                        COVERAGE_FRACTION=f"{both.mean():.4f}")

    return {"coverage_fraction": float(both.mean()),
            "a_only_cells": int(np.count_nonzero(np.isfinite(za) & ~np.isfinite(zb))),
            "b_only_cells": int(np.count_nonzero(np.isfinite(zb) & ~np.isfinite(za)))}

Using NaN as the nodata value throughout, rather than a sentinel, is what removes the poisoning class of bug entirely: NaN propagates through arithmetic by definition, and np.nanmean and friends ignore it without being told. The one cost is that integer rasters cannot carry it, which is not a constraint for elevation data.

Reporting a_only_cells and b_only_cells separately matters because they mean different things. Cells present in the older survey and missing from the newer one usually indicate a reconstruction that got worse — new vegetation, worse light. The reverse usually means coverage improved. Either way, a difference raster that silently lost a quarter of the site should say so.

# What is actually in the difference raster

A DEM of difference is not a map of change. It is a map of change plus three other things, and separating them is the whole job.

Figure 3 — Four contributions, one raster.

# Edge-case matrix

# Verification snippet

The check that matters is whether the difference behaves like a difference: near zero over stable ground, with a spread consistent with the two inputs’ own accuracy.

import numpy as np


def sanity_check_difference(diff: np.ndarray, stable_mask: np.ndarray,
                            expected_sigma: float) -> dict:
    """Does this difference raster look like two good surveys disagreeing?"""
    s = diff[stable_mask & np.isfinite(diff)]
    if s.size < 200:
        raise ValueError("not enough stable cells to judge the difference")

    bias = float(np.median(s))
    sigma = float(np.median(np.abs(s - bias)) * 1.4826)
    problems = []
    if abs(bias) > 0.02:
        problems.append(f"stable ground offset by {bias:+.3f} m — datum or registration")
    if sigma > 2.5 * expected_sigma:
        problems.append(f"stable scatter {sigma:.3f} m is far above the expected "
                        f"{expected_sigma:.3f} m — check the resampling and grids")
    if sigma < 0.3 * expected_sigma:
        problems.append("stable scatter is suspiciously low — the two surfaces "
                        "may share a source or be the same file")
    return {"stable_bias_m": bias, "stable_sigma_m": sigma, "problems": problems}

The suspiciously-low branch catches a real and embarrassing failure: differencing a survey against itself, or against a copy that was never actually reprocessed. It produces a beautiful, clean, entirely meaningless map, and without the check nobody notices until a client asks why nothing has changed in four months.

Figure 2 — Coverage is part of the result. A shrinking differenced area is often the most important thing a monitoring run has to say.

# Choosing a resampling method that does not invent change

Resampling one surface onto the other’s grid is unavoidable, and the method chosen has a direct effect on the apparent change. Bilinear is the right default for elevation: it is continuous, cheap, and its error is proportional to local curvature. Cubic is smoother and overshoots at breaks, which manufactures paired positive and negative change along every kerb and bench edge — precisely the features a monitoring run cares about.

Nearest neighbour has one legitimate use here and it is worth knowing: when the two grids are offset by exactly zero and the resolutions match, nearest is a no-op and therefore lossless. Asserting that condition and taking the fast path avoids resampling error entirely on the common case of two surveys gridded from the same template.

import numpy as np


def grids_are_identical(a, b, tol: float = 1e-9) -> bool:
    """True when two rasters share a grid exactly, so no resampling is needed."""
    return (a.crs == b.crs and a.shape == b.shape
            and all(abs(x - y) < tol for x, y in zip(a.transform, b.transform)))

Making every survey of a site use one fixed grid template — origin, resolution and extent stored with the site, exactly as the base surface is in computing volumes and stockpiles in Python — turns this from a judgement into a check. Every subsequent difference is then exact, every monitoring point samples the same cells, and the corduroy artefact cannot occur.

# When to escalate

• The stable-ground scatter is far above the surveys’ own accuracy. Something other than resampling is wrong: most likely the two epochs used different control, which is a georeferencing problem rather than a raster one.
• Coverage dropped sharply between epochs. The newer flight is worse. Investigate the flight rather than patching the raster; interpolating across the loss produces a surface that differences cleanly and means nothing.
• The difference must be measured on steep faces. A vertical difference is the wrong quantity there; move to a cloud-to-cloud distance.

# Related

• Change Detection Between Survey Epochs
• Separating Real Change from Noise with Detection Limits
• Fixing NoData and Alpha Band Confusion in Exports

← Change Detection Between Survey Epochs