Change Detection and Differencing Workflows

Change detection is the most requested product in remote sensing and the easiest to do badly. Subtracting two dates is one line of code; producing a map where the highlighted pixels correspond to real change on the ground is a different problem, dominated by radiometric comparability, registration and threshold choice. This topic covers that gap, inside the wider context of Satellite Processing Workflows & Index Pipelines.

The specific challenge is that everything varies between two acquisitions, not just the surface. Sun elevation differs, so shadows move. Atmosphere differs, so haze changes contrast. Phenology differs, so healthy vegetation legitimately looks different in June and August. A difference image contains all of that plus the change you care about, and the workflow’s job is to suppress the first three.

Prerequisites

pip install "rioxarray>=0.15" "xarray>=2023.6" "rasterio>=1.3.0" "scikit-image>=0.22" "numpy>=1.23"

Conceptually you need masking from Cloud and Shadow Masking Strategies, alignment from Aligning Two Rasters with reproject_match, and index construction from Spectral Index Calculation Pipelines.

What a difference image actually contains

Before choosing a method, it is worth being explicit about what is in the signal you are about to threshold.

Step-by-step workflow

1. Pair the dates deliberately

The strongest control available is date selection. Anniversary pairs — the same week in successive years — remove most phenological difference for free, and similar sun elevation removes most of the shadow difference. Where single scenes are too cloudy, composite a short window around each anniversary rather than reaching for whatever is clear.

import pystac_client

catalog = pystac_client.Client.open("https://example-stac-api.org/v1")

def anniversary_pair(bbox, tile, year_a: int, year_b: int, month: int = 7):
    """Items from the same month in two years, cloud-filtered, same tile."""
    def search(year):
        return list(catalog.search(
            collections=["sentinel-2-l2a"], bbox=bbox,
            datetime=f"{year}-{month:02d}-01/{year}-{month:02d}-28",
            query={"eo:cloud_cover": {"lt": 20}, "s2:mgrs_tile": {"eq": tile}},
        ).items())
    return search(year_a), search(year_b)

2. Align both dates and take the union of their masks

import rioxarray
import xarray as xr
from rasterio.enums import Resampling


def load_pair(before_path: str, after_path: str, scl_before: str, scl_after: str):
    """Two dates on one grid, with a shared validity mask."""
    before = rioxarray.open_rasterio(before_path, masked=True).squeeze(drop=True)
    after = rioxarray.open_rasterio(after_path, masked=True).squeeze(drop=True)
    after = after.rio.reproject_match(before, resampling=Resampling.bilinear)

    bad = [0, 1, 3, 8, 9, 10, 11]
    m1 = rioxarray.open_rasterio(scl_before).squeeze(drop=True).rio.reproject_match(
        before, resampling=Resampling.nearest)
    m2 = rioxarray.open_rasterio(scl_after).squeeze(drop=True).rio.reproject_match(
        before, resampling=Resampling.nearest)

    valid = ~m1.isin(bad) & ~m2.isin(bad)      # change needs BOTH dates observed
    return before.where(valid), after.where(valid), valid

Taking the union of the two masks is what stops a cloud on one date being reported as change. It also means the output has a smaller valid footprint than either input, which should be recorded rather than hidden.

3. Difference an index, not raw bands

def normalised_difference(a, b):
    denom = a + b
    return xr.where(denom != 0, (a - b) / denom, float("nan"))

ndvi_before = normalised_difference(nir_before, red_before)
ndvi_after = normalised_difference(nir_after, red_after)
d_ndvi = ndvi_after - ndvi_before          # positive = greening, negative = loss

Normalised indices are ratios, so a multiplicative illumination difference between dates largely cancels. Raw band differences do not have that property, which is why a band-difference change map is dominated by sun angle.

4. Choose the threshold from stable areas

import numpy as np

stable = d_ndvi.where(stable_mask)                  # areas known not to have changed
sigma = float(stable.std())
threshold = 2.5 * sigma                             # in the units of the index

loss = d_ndvi < -threshold
gain = d_ndvi > threshold
print(f"sigma over stable ground: {sigma:.3f}, threshold: {threshold:.3f}")

A threshold derived this way is comparable across scenes and seasons, because it is expressed in units of the noise actually present in that pair. A fixed number is not.

5. Remove speckle with a minimum mapping unit

import numpy as np
from skimage.morphology import remove_small_objects

MMU_PIXELS = 10          # 10 pixels at 10 m ≈ 0.1 ha

loss_clean = remove_small_objects(loss.values.astype(bool), min_size=MMU_PIXELS)

The minimum mapping unit is a stated product property, not a tuning knob: publish it alongside the map, because a change map with an unstated MMU cannot be compared with any other.

Parameter reference

Verification and testing

Three checks separate a plausible change map from a defensible one.

The first is a stable-site check. Pick areas that certainly did not change — established forest, water bodies, paved surfaces — and confirm that fewer than a small percentage of their pixels are flagged. A stable site with five percent change flagged is telling you the threshold is too tight or the registration is poor.

The second is a known-change check. If any change is documented — a harvested block, a new development, a fire perimeter from an independent source — confirm it is detected and that its mapped extent is sensible. Detection without extent accuracy is common when the threshold is too loose.

The third is a direction check. Change maps should be signed, and the sign should make sense: vegetation loss must be negative in an NDVI difference and positive in an NBR difference of burn severity. A map where both directions appear scattered together usually indicates a registration or masking problem rather than a busy landscape.

Choosing between two-date and time-series methods

Two-date differencing answers “what changed between these dates”. It is cheap, easy to explain, and appropriate for abrupt high-magnitude events: fire, clear-felling, flooding, construction. Its weakness is that a single pair cannot distinguish a real trend from the noise of two particular acquisitions, so gradual change — thinning, degradation, slow encroachment — is either invisible or indistinguishable from phenology.

Time-series methods answer “how is this pixel behaving”. Fitting a seasonal model over many dates and flagging departures from it separates trend from cycle, detects gradual change, and reports a date of change rather than an interval. The cost is data volume, compute and complexity — and a much greater sensitivity to gaps, which is where the interpolation questions in Temporal Aggregation and Time-Series Analysis become decisive.

A practical middle path is a composite pair: build a median composite over a month around each anniversary rather than picking single scenes. That removes most cloud and phenological noise for a fraction of the effort of a full time-series method, and it keeps the two-date interpretation that stakeholders find easy to read.

Pairing strategies, compared

How the two dates are chosen determines more of the result than any later parameter, and there are three practical strategies with different costs.

For an operational monitoring product the composite pair is usually the right default, because it is the only option whose behaviour does not depend on the weather in a particular week. The exception is rapid response — mapping a fire or a flood within days — where the event date is the point and a composite would average it away. The compositing mechanics are the same as in Creating Monthly NDVI Composites with xarray resample, and the choice of reducer matters as much here as it does there.

Reporting a change map so it can be used

A change map is a claim about the ground, and it travels further than the person who made it. Five pieces of metadata determine whether a downstream user can act on it.

The dates, obviously — but both the acquisition dates and, if composited, the window and the number of contributing observations per pixel. A change detected between a four-observation composite and a one-observation composite is not as reliable as one between two well-populated composites, and only the count reveals that.

The index and the threshold, in the index’s own units, together with how the threshold was derived. “NDVI difference below −0.15, being 2.5σ over stable forest in this pair” is reproducible. “significant vegetation loss” is not.

The minimum mapping unit, because it sets the smallest feature the map can contain and therefore what its absence means. A user searching for 0.05 ha clearings in a product with a 0.1 ha MMU will find nothing and may conclude, wrongly, that there is nothing there.

The valid footprint, which is the intersection of both dates’ clear-sky areas rather than the scene extent. Areas excluded because one date was cloudy must be distinguishable from areas assessed and found unchanged — a distinction that is invisible in a two-colour map and essential in a table.

The direction convention. Sign errors are common and easy to make when an index is inverted; state explicitly that negative values mean loss, and include a legend that a non-specialist can read.

Carrying those five fields costs a few tags in the output and turns a picture into evidence, which is the same discipline argued for in Writing and Validating Cloud-Optimized GeoTIFFs.

---

Running change detection at scale

A single pair is a script. A national product is thousands of pairs across hundreds of tiles, and three decisions determine whether it finishes.

The unit of work should be the tile-pair, not the scene and not the country. Each tile-pair is independent — it reads two dates for one tile, writes one change layer, and shares nothing with its neighbours — which makes it a natural fit for the container-per-unit model described in Processing COGs on AWS Batch with Docker. Anything that needs cross-tile context, such as an area statistic, belongs in a second pass over the outputs rather than inside the first.

The threshold has to be decided consistently across tiles, and this is where naive scaling breaks. Deriving σ per tile means neighbouring tiles get different thresholds, which produces visible seams along tile boundaries in the mosaicked product — the same class of artefact treated in Seamless Mosaicking and Edge Blending. The fix is to derive σ once over a representative sample spanning the whole area and pass it in as a parameter, so every tile applies the same rule.

Intermediates are usually worth keeping. Writing the continuous difference layer as well as the thresholded map costs storage but means a threshold change does not require re-reading the imagery — a rerun becomes a cheap pass over the differences instead of an expensive pass over the archive. Given how often thresholds are revised after the first review, that trade is almost always favourable.

Two smaller points save real time. Compute the difference at the coarsest resolution the product requires rather than at the finest available: a 20 m change map derived from 20 m inputs costs a quarter of the 10 m version and, for most reporting units, says the same thing. And validate on a sample of tiles before launching the full run — a systematic registration or masking problem is identical in every tile, and finding it after ten thousand outputs have been written costs the whole run.

---

Sensor differences when the pair spans missions

Occasionally the before and after dates come from different sensors — a Landsat archive image paired with a recent Sentinel-2 acquisition. This is legitimate for coarse, high-magnitude change, and it carries two extra caveats worth stating in the product.

Band centres differ, so the same index is not numerically identical between missions even over unchanged ground. The offset is small for NDVI and larger for indices that use narrow bands, and it biases the difference in a constant direction. Estimating that offset over stable ground and subtracting it before thresholding removes most of the problem.

Resolutions differ, so one date must be resampled onto the other’s grid, and the resampled date is smoother. Fine-scale change near the resolution limit of the coarser sensor will be detected inconsistently, which argues for a larger minimum mapping unit than a single-mission pair would need.

Troubleshooting

Every field boundary is flagged

Sub-pixel misregistration. Co-register the pair before differencing, and check by differencing a date against itself resampled — the result should be zero everywhere.

The whole scene shifts in one direction

A radiometric offset between dates, usually different processing baselines or an atmospheric difference. Normalise using stable targets, or difference an index rather than reflectance.

Change appears only in shadowed terrain

Sun elevation differs substantially between the dates. Prefer anniversary pairs, and consider masking steep slopes where illumination varies most.

The map is dominated by isolated single pixels

The threshold is too loose, the MMU filter is missing, or both. Tighten to 3σ and apply a minimum mapping unit before publishing.

Results change when the pair is re-run

Non-deterministic masking — often a cloud probability threshold applied to a model that was re-run — or an unpinned resampling default. Pin both, and record them in the output metadata.

Frequently Asked Questions

Q: Why does my change map light up along every field boundary? Sub-pixel misregistration between the two dates. Edges are where the gradient is steepest, so a fraction of a pixel of shift produces a large difference there and almost none in homogeneous interiors. Co-register before differencing.

Q: Should I difference reflectance or an index? An index, in almost every case. Normalised indices cancel much of the illumination and atmospheric difference between dates, whereas raw band differences respond strongly to sun angle and haze.

Q: How do I choose the change threshold? From the data. Compute the difference over areas known to be stable, take its standard deviation, and set the threshold at two or three of those units. A round number like 0.2 has no meaning across sensors, seasons or land covers.

Q: Is one image pair enough? For abrupt, high-magnitude events such as fire or clear-felling, often yes. For gradual change, a two-date difference cannot separate trend from noise, and a time-series approach over many dates is the honest method.

---

Related

• Computing an NDVI Difference Between Two Dates — the minimal correct implementation of this workflow.
• Thresholding Change Maps and Removing Noise — turning a continuous difference into a defensible map.
• Mapping Burned Area with NBR Differencing — the best-established application of the method.
• Cloud and Shadow Masking Strategies — the masking both dates depend on.
• Temporal Aggregation and Time-Series Analysis — the alternative when change is gradual rather than abrupt.