Thresholding Change Maps and Removing Noise

To turn a continuous difference into a binary change map, measure the noise over stable ground, cut at a multiple of it, then remove objects below the minimum mapping unit:

import numpy as np
from skimage.morphology import remove_small_objects

sigma = float(d_ndvi.where(stable_mask).std())
loss = (d_ndvi < -2.5 * sigma).values

loss_clean = remove_small_objects(loss, min_size=10)     # 10 px at 10 m ≈ 0.1 ha
print(f"sigma={sigma:.3f}  threshold={-2.5 * sigma:.3f}  patches kept={loss_clean.sum()} px")

This is the step that turns the difference layer from Computing an NDVI Difference Between Two Dates into a product.

Why This Arises in Remote Sensing Workflows

A continuous difference layer is honest but unusable: nobody can act on “this pixel changed by −0.07”. Turning it into a map means committing to a boundary between change and noise, and that commitment is where most change products acquire their errors.

Two habits cause most of the damage. The first is a fixed threshold copied from a paper — 0.1, 0.2, 0.15 — applied regardless of sensor, season or land cover. Such a value encodes the noise level of the pair it was derived from, and applied elsewhere it is either far too tight or far too loose. The second is skipping the cleanup, which leaves a map speckled with isolated pixels that no field visit could confirm and that inflate the reported area of change dramatically.

Both are avoidable with a few lines of code, and both are visible in the output if anyone looks: an unclean map has a characteristic salt-and-pepper texture, and a badly chosen threshold either flags nothing or flags a quarter of the tile.

Threshold choice, drawn on the map At one sigma the map is dominated by speckle covering nearly a third of the tile. At two sigma real patches emerge with some speckle around them. At three sigma only the largest, strongest patches survive, and small genuine changes are lost. The middle option plus a minimum mapping unit is usually the defensible choice. 1σ — 31% flagged 2.5σ — 2.4% flagged 3.5σ — 0.4% flagged real patches buried in speckle patches clear, little speckle left patches shrunk; small change lost The same data in all three panels — only the cut moved. This is why the threshold must be reported.

Environment & Setup

Package Version Why
scikit-image ≥0.22 Morphology, remove_small_objects, label
numpy ≥1.23 Array handling and statistics
rioxarray ≥0.15 Reading the difference layer and writing the map
scipy ≥1.11 Optional: distance transforms and labelling 
pip install "scikit-image>=0.22" "scipy>=1.11" "rioxarray>=0.15"

Complete Working Example

This function takes a continuous difference and a stable-ground mask and returns a cleaned, signed change map along with the parameters used to make it.

Expected flagged fraction by threshold multiple Assuming approximately normal noise over stable ground, a threshold at two sigma flags about 4.6 percent of unchanged pixels, at two and a half about 1.2 percent, and at three about 0.3 percent. A minimum mapping unit removes most of what remains because noise is isolated. What a threshold costs in false positives threshold stable pixels flagged after a 10-pixel MMU 1.5σ 13.4% ~2.1% 2.0σ 4.6% ~0.4% 2.5σ 1.2% ~0.05% 3.0σ 0.3% ~0.01% Noise is isolated and real change is contiguous — which is exactly why the MMU works so well. 
import numpy as np
import rioxarray  # noqa: F401
import xarray as xr
from skimage.morphology import binary_opening, disk, remove_small_objects


def threshold_change(
    diff: xr.DataArray,
    stable_mask: xr.DataArray,
    *,
    k: float = 2.5,
    mmu_pixels: int = 10,
    opening_radius: int = 1,
) -> tuple[xr.DataArray, dict]:
    """Signed change map: −1 loss, +1 gain, 0 no change, NaN unassessed."""
    stable = diff.where(stable_mask)
    sigma = float(stable.std(skipna=True))
    centre = float(stable.median(skipna=True))       # residual offset between dates
    if not np.isfinite(sigma) or sigma == 0:
        raise ValueError("stable-ground sigma is undefined — check the stable mask")

    lo, hi = centre - k * sigma, centre + k * sigma

    loss = (diff < lo).fillna(False).values
    gain = (diff > hi).fillna(False).values

    # Opening first removes single-pixel noise; MMU then drops sub-threshold patches
    selem = disk(opening_radius)
    loss = remove_small_objects(binary_opening(loss, selem), min_size=mmu_pixels)
    gain = remove_small_objects(binary_opening(gain, selem), min_size=mmu_pixels)

    out = xr.zeros_like(diff, dtype="float32")
    out = out.where(diff.notnull())                  # unassessed stays NaN
    out.values[loss] = -1.0
    out.values[gain] = 1.0

    params = {
        "sigma_stable": round(sigma, 4),
        "centre_offset": round(centre, 4),
        "k": k,
        "threshold_low": round(lo, 4),
        "threshold_high": round(hi, 4),
        "mmu_pixels": mmu_pixels,
        "opening_radius": opening_radius,
        "loss_pixels": int(loss.sum()),
        "gain_pixels": int(gain.sum()),
    }
    out.attrs.update(params)
    return out, params


if __name__ == "__main__":
    diff = rioxarray.open_rasterio("dndvi.tif", masked=True).squeeze(drop=True)
    stable = rioxarray.open_rasterio("stable_forest_mask.tif").squeeze(drop=True).astype(bool)

    change, params = threshold_change(diff, stable)
    print(params)
    change.rio.write_nodata(np.nan, inplace=True)
    change.rio.to_raster("change_map.tif", driver="COG", compress="DEFLATE", blocksize=512)

Two design choices are worth calling out. Centring the thresholds on the stable-ground median rather than on zero absorbs any residual radiometric offset between the dates, so a pair with a small systematic difference does not report half the tile as change. And returning the parameters rather than only the map means they can be written into the output and compared between runs.

Variant Patterns

1. Histogram-based thresholds, and when they apply

Otsu’s method finds the cut that best separates a bimodal distribution, which is exactly right for a fire scar covering a third of a scene and exactly wrong for a scene where two percent changed.

import numpy as np
from skimage.filters import threshold_otsu

values = diff.values[np.isfinite(diff.values)]
changed_fraction_guess = float((values < -0.05).mean())

if changed_fraction_guess > 0.15:
    cut = threshold_otsu(values)                 # genuinely bimodal — Otsu is appropriate
else:
    cut = float(np.nanmedian(values) - 2.5 * np.nanstd(values))

The guard is the point. Applying Otsu unconditionally produces a threshold that varies wildly between tiles of the same product, which is how a mosaicked change map acquires tile-shaped discontinuities.

2. Two thresholds: seed and grow

A single cut forces one trade-off between missing weak change and admitting noise. Hysteresis uses a strict threshold to find confident cores and a looser one to grow them, which recovers the fringe of a real patch without admitting isolated noise.

import numpy as np
from skimage.morphology import reconstruction

strict = (diff < centre - 3.0 * sigma).fillna(False).values
loose = (diff < centre - 1.8 * sigma).fillna(False).values

# Grow the strict seeds inside the loose mask only
seeds = np.where(strict, 1.0, 0.0)
grown = reconstruction(seeds, np.where(loose, 1.0, 0.0), method="dilation") > 0

This is usually the single largest quality improvement available over a plain threshold, and it costs one extra parameter to report.

3. Cleaning without destroying real shapes

What each cleanup step removes Opening with a small structuring element removes isolated pixels and one-pixel-wide connections but also erodes thin linear features such as firebreaks and tracks. A minimum mapping unit removes small objects regardless of shape, preserving thin features that are long enough. Applying both in that order keeps compact patches intact. Three feature types, two cleanup operations feature opening (radius 1) MMU (10 px) isolated single pixels removed removed thin linear feature (track, firebreak) eroded or destroyed kept if long enough compact patch (clearing, burn scar) kept, edges smoothed kept If linear change matters to the product, skip the opening and rely on the minimum mapping unit alone. The MMU should be stated in hectares as well as pixels, since consumers think in area.

Validating the Map

Once thresholded, the map can be checked against things you already know, which is easier than validating a continuous layer.

Over stable reference areas, the flagged fraction should be small and roughly predictable from the threshold: about 1.2 percent at 2.5σ if the noise is normal, and much less after the minimum mapping unit removes isolated pixels. A stable area with ten percent flagged means the noise is not what the stable mask suggested, usually because the mask includes something that genuinely changed.

Over known change, the patch should be detected and its area should be sensible. Systematic under-measurement of area points to a threshold that is too strict or an over-aggressive opening; systematic over-measurement usually means the mask fringe was not buffered, so cloud edges are being counted as change — the buffering point from Masking Clouds with the Sentinel-2 SCL Band.

Across tiles, the flagged fraction should vary smoothly. A histogram of per-tile change fraction with a long right tail identifies tiles where something went wrong — a bad date pair, a registration failure, an unmasked cloud — far faster than looking at maps.

Common Errors

Half the tile is flagged as change

The threshold was derived from a stable mask that is not stable, or the difference has a large offset. Check the stable-ground median first; if it is far from zero, fix the radiometry before thresholding.

The map is empty despite obvious change

k is too large, or the sign is inverted — a loss test applied to a difference computed as before minus after finds nothing. Print the min and max of the difference before thresholding.

remove_small_objects raises on a float array

It expects a boolean or labelled integer array. Cast with .astype(bool) after filling NaN, as in the example.

Frequently Asked Questions

Q: Is Otsu’s method a good default for change maps? Only when change occupies a substantial share of the scene. Otsu assumes a bimodal histogram; in a typical scene where under two percent of pixels changed, the histogram has one mode and Otsu places the cut arbitrarily inside the noise.

Q: Should I clean before or after thresholding? After. Smoothing the continuous difference before thresholding blurs real edges and shifts the apparent boundary of a change patch. Morphology on the binary map removes speckle without moving edges.

Q: What minimum mapping unit should I use? Whatever the application can act on, stated explicitly. A forestry product that reports at 0.5 ha should not emit 3-pixel patches; a rapid-response fire product may want everything above one pixel. There is no correct value, only a declared one.