⚙️ How-To Satellite imagery Raster & rasterio Satellite Imagery · step 15 of 16 Updated

How to Detect Change Between Two Satellite Images in Python

Problem statement

Differencing two dates is one line and produces a change map full of things that did not change. The signal is buried under four sources of apparent change, and all four are larger than most real changes:

  • cloud and shadow — 51% of a single Sentinel-2 date over Snowdonia
  • illumination — different sun angle, different date, different atmosphere
  • misregistration — a half-pixel shift puts a boundary in a different place
  • phenology — the same field in June and September genuinely differs

Only the first is easy to remove. The rest are what separates a change map from a difference map.

Quick answer

Mask both dates, difference an index rather than a band, and threshold from the data:

import numpy as np

usable = mask_a & mask_b                       # both dates clear
difference = np.where(usable, ndvi_b - ndvi_a, np.nan)

stable = difference[np.isfinite(difference)]
threshold = 2.5 * np.median(np.abs(stable - np.median(stable))) * 1.4826
print(f"robust sd {threshold / 2.5:.4f}, threshold ±{threshold:.4f}")

change = np.abs(difference) > threshold
print(f"{np.nansum(change):,} cells flagged ({np.nanmean(change):.2%})")

A robust scale estimate — the median absolute deviation — is essential. Using the ordinary standard deviation lets the change itself inflate the threshold that is supposed to detect it.

Four sources of apparent change: cloud and shadow, illumination differences, misregistration and phenology, all larger than typical real change.
Differencing detects all four. Only the first has a clean fix.

Step-by-step solution

1. Mask both dates, and intersect the masks

A cell is comparable only if it was clear on both dates. The usable area is the intersection, which is much smaller than either mask alone.

On a pair of Snowdonia scenes with 49% and 62% usable individually, the intersection could be as low as 11% in the worst case. Report it — a change map covering a third of the area is a different product from one covering all of it.

2. Difference an index, not a raw band

A raw band difference is dominated by illumination. A normalised index divides out anything scaling all bands equally — sun angle, slope, thin haze — which is most of the nuisance variation.

NDVI for vegetation change, NBR for burn severity, NDWI for water extent. Which index depends on the change you are looking for, and the choice determines what you can detect.

3. Check the co-registration

A half-pixel shift produces a bright and dark pair along every boundary — the classic misregistration signature. It looks like change along every field edge and road.

from scipy.signal import correlate2d
shift = np.unravel_index(np.argmax(correlate2d(a - a.mean(), b - b.mean(),
                                               mode="same")), a.shape)

Sentinel-2 tiles from the same relative orbit are well co-registered. Across orbits, sensors or providers, check.

4. Set the threshold from the stable population

Most of the scene did not change, so the distribution of differences is a narrow peak around zero with tails.

The median absolute deviation, scaled by 1.4826, estimates the standard deviation of that peak without being inflated by the tails. Anything beyond about 2.5 of those is a candidate.

Using the ordinary standard deviation instead is circular: a large real change widens the distribution and raises the threshold that should have detected it.

5. Separate change from phenology

Two dates in different seasons differ everywhere. Comparing the same season across years — an anniversary pair — removes most of it.

Where that is not possible, compare each date against a seasonal baseline built from several years rather than against each other.

A half-pixel shift producing paired bright and dark edges along every boundary, mimicking change.
Misregistration produces a bright-dark pair along every edge. It is unmistakable once you know the pattern.

Code examples

Example 1 — change detection with the diagnostics 

import numpy as np


def detect_change(index_a, index_b, mask_a, mask_b, sigmas=2.5,
                  min_patch=5):
    """Robust thresholded difference, with the usable area reported."""
    usable = mask_a & mask_b
    print(f"  date A usable {mask_a.mean():.1%}, date B {mask_b.mean():.1%}, "
          f"both {usable.mean():.1%}")
    if usable.mean() < 0.2:
        print("  ! under a fifth of the scene is comparable")

    difference = np.where(usable, index_b - index_a, np.nan)
    values = difference[np.isfinite(difference)]

    median = float(np.median(values))
    mad = float(np.median(np.abs(values - median)))
    robust_sd = mad * 1.4826
    ordinary_sd = float(values.std())
    threshold = sigmas * robust_sd

    print(f"  difference median {median:+.4f}")
    print(f"  robust sd {robust_sd:.4f}, ordinary sd {ordinary_sd:.4f} "
          f"({ordinary_sd / max(robust_sd, 1e-9):.2f}x)")
    print(f"  threshold ±{threshold:.4f}")
    if abs(median) > robust_sd:
        print("  ! the whole scene has shifted — check calibration, "
              "illumination or season")

    change = np.abs(difference - median) > threshold
    increase = (difference - median) > threshold
    decrease = (difference - median) < -threshold

    if min_patch > 1:
        from scipy.ndimage import label
        labels, n = label(change)
        sizes = np.bincount(labels.ravel())
        small = np.isin(labels, np.nonzero(sizes < min_patch)[0])
        removed = int((change & small).sum())
        change = change & ~small
        print(f"  removed {removed:,} isolated cells in patches under "
              f"{min_patch} cells")

    print(f"  change {np.nansum(change):,} cells ({np.nanmean(change):.2%}): "
          f"{np.nansum(increase & change):,} increase, "
          f"{np.nansum(decrease & change):,} decrease")
    return change, difference, threshold

The ratio of ordinary to robust standard deviation is a free diagnostic. A ratio near 1 means little real change; a ratio of 2 or more means the tails are heavy, which is either real change or unmasked cloud.

Example 2 — checking the co-registration 

import numpy as np
from scipy.signal import fftconvolve


def check_registration(a, b, max_shift=5):
    """Cross-correlate two dates to find a sub-pixel offset."""
    a = np.nan_to_num(a - np.nanmean(a))
    b = np.nan_to_num(b - np.nanmean(b))

    correlation = fftconvolve(a, b[::-1, ::-1], mode="same")
    centre = np.array(correlation.shape) // 2
    window = correlation[centre[0] - max_shift:centre[0] + max_shift + 1,
                         centre[1] - max_shift:centre[1] + max_shift + 1]
    peak = np.unravel_index(np.argmax(window), window.shape)
    offset = (peak[0] - max_shift, peak[1] - max_shift)

    print(f"  peak correlation at offset {offset} pixels")
    if offset != (0, 0):
        print("  ! the dates are misregistered — every boundary will show "
              "as change")
    return offset

Sub-pixel misregistration will not show as a whole-pixel offset here, and it still produces edge artefacts. If the change map shows a bright-dark pair along every field boundary, that is the diagnosis regardless of what the correlation says.

Example 3 — a stable-area normalisation 

import numpy as np


def normalise_to_stable(index_a, index_b, usable, quantile=0.6):
    """Align the two dates using the cells least likely to have changed."""
    difference = np.where(usable, index_b - index_a, np.nan)
    magnitude = np.abs(difference - np.nanmedian(difference))
    cutoff = np.nanquantile(magnitude, quantile)
    stable = usable & (magnitude <= cutoff)

    a_values = index_a[stable]
    b_values = index_b[stable]
    slope, intercept = np.polyfit(a_values, b_values, 1)

    print(f"  {stable.mean():.1%} of cells treated as stable")
    print(f"  b ≈ {slope:.4f} * a + {intercept:+.4f}")
    if abs(slope - 1) > 0.1 or abs(intercept) > 0.05:
        print("  ! substantial systematic difference between the dates")

    corrected = (index_b - intercept) / slope
    return corrected, {"slope": float(slope), "intercept": float(intercept),
                       "stable_fraction": float(stable.mean())}

Relative radiometric normalisation using pseudo-invariant features is the standard remedy for a systematic offset between dates. Fitting on the most stable 60% of cells and applying the correction to everything removes a global shift without erasing real change.

It cannot fix a spatially varying difference — a haze gradient, a partial cloud shadow — which is why masking comes first.

Explanation

Why differencing an index beats differencing a band

A raw band difference includes every multiplicative effect: sun angle, atmospheric transmission, sensor calibration, topographic illumination.

A normalised difference index divides by the sum of its bands, which cancels anything scaling both equally. That removes most of the illumination variation at the cost of some sensitivity.

The residual difference is then dominated by things that changed the ratio between bands — which is what a change in surface type does.

Why the robust threshold matters

The distribution of differences over a scene is a narrow peak of unchanged cells with tails of real change.

The ordinary standard deviation is computed over the whole distribution, so the tails inflate it. A scene with substantial change produces a large standard deviation, a large threshold, and a change map that misses the change.

The median absolute deviation is unaffected by the tails, so the threshold reflects the noise among the unchanged cells — which is exactly what it should be measuring.

Why misregistration looks like change along every edge

A half-pixel shift means each boundary pixel contains a different mixture of the two surfaces on each date. Along one side of a boundary the value rises; along the other it falls.

The change map then shows a paired bright and dark line along every field edge, road and coastline — a pattern nothing physical produces.

Once recognised it is unmistakable, and it is why every change detection workflow should start by looking at the difference image along a few known boundaries.

Why two dates is the weakest possible design

Two observations can tell you that something differs. They cannot distinguish a permanent change from a seasonal fluctuation, a temporary flood, or a single bad observation.

A time series can. Fitting a seasonal model and looking for departures from it separates trend from cycle and gives every detection a date rather than an interval.

Where only two dates are available, an anniversary pair — the same time of year, different years — removes most of the seasonal component and is much stronger than two arbitrary dates.

An ordinary standard deviation inflated by the change against a robust median absolute deviation that is not.
Using the ordinary standard deviation is circular: the change raises its own threshold.

Edge cases or notes

  • Intersect the masks; a cell is comparable only if clear on both dates.
  • Difference an index, not a band.
  • Use a robust scale estimate; the ordinary standard deviation is inflated by the change itself.
  • Check co-registration. Paired bright and dark edges are the signature.
  • Use anniversary pairs to remove seasonality.
  • Remove isolated single-cell detections; they are almost all noise.
  • Report the comparable area — often far less than either date alone.
  • A non-zero median difference means a scene-wide shift, not change.

FAQ

How do I detect change between two satellite images?

Mask both dates, intersect the masks, difference a normalised index rather than a raw band, and threshold using a robust scale estimate from the unchanged cells.

Why is my whole scene flagged as changed?

Usually an illumination or calibration difference producing a non-zero median. Normalise the dates against stable cells before thresholding.

What threshold should I use?

About 2.5 times a robust standard deviation — the median absolute deviation times 1.4826 — computed from the difference image itself.

Why does change follow every field boundary?

Misregistration. A sub-pixel shift makes each boundary pixel a different mixture on each date, producing paired bright and dark lines.

Should I difference bands or indices?

Indices. A normalised difference cancels effects that scale all bands equally, which is most of the illumination variation.

How do I avoid detecting seasonal change?

Use an anniversary pair — the same time of year in different years — or compare each date against a multi-year seasonal baseline.

Is two dates enough?

It is the weakest design. Two observations cannot separate permanent change from a seasonal fluctuation or a single bad observation.