Change Detection in Satellite Imagery Over Time

Subtract one dated scene from another and threshold the result. The arithmetic takes a line; everything that determines whether the output means anything happens before the subtraction.

Co-registration comes first

Differencing assumes pixel (i, j) in both scenes covers the same ground. Where it does not, the difference image shows a bright rim around every high-contrast edge — roads, field boundaries, rooftops — because a road pixel in one scene sits against a soil pixel in the other.

The signature is diagnostic rather than subtle: misregistration artefacts follow edges, are paired bright-and-dark on opposite sides of each edge, and have the same displacement direction everywhere. If your change map looks like an outline drawing of the landscape, the problem is registration and no threshold will fix it.

Modern products are orthorectified to a common grid and are usually co-registered to well under a pixel, but this is a property to verify rather than assume, particularly across sensors or across processing baselines. The standard measurement is phase correlation on overlapping windows, which gives a sub-pixel shift estimate per window; if those shifts are not near zero and not spatially uniform, the terrain model used for orthorectification is the likely cause and the fix is upstream, not in your differencing code.

The sun moved, and that dominates

Two scenes of the same unchanged ground, acquired months apart, do not have the same brightness. The largest single reason is solar geometry. For a roughly Lambertian surface the irradiance on the ground scales with the cosine of the solar zenith angle:

scene 1, June:     solar zenith 30 deg   cos = 0.866
scene 2, October:  solar zenith 55 deg   cos = 0.574

ratio = 0.574 / 0.866 = 0.663

apparent brightness of unchanged ground falls by 33.7%

A third of the signal, with nothing on the ground having changed. Most real land-cover changes produce a smaller difference than that, so an uncorrected difference image is mostly a map of solar geometry.

Three defences, in order of preference. Use surface reflectance products — Sentinel-2 Level-2A, Landsat Collection 2 Level-2 — where the illumination and atmospheric correction has already been applied by the provider. Failing that, apply relative radiometric normalisation: pick pseudo-invariant features that should not change between dates (bare rock, large roofs, deep water, quarry surfaces), fit a per-band linear regression of scene 2 against scene 1 on those pixels only, and apply it to the whole scene. And regardless of both, pair anniversary dates — the same week of the year — so that solar geometry, phenology and shadow length are as close as they can be before any correction runs.

Terrain adds a second layer of this. A north-facing slope and a south-facing slope receive different irradiance, and the difference between them changes with season, so in mountainous terrain uncorrected differencing produces change that follows the topography. Topographic correction is a separate step and is not applied by default in most standard products.

The differencing arithmetic

With surface reflectance in hand, take one pixel — a parcel of vegetation cleared between two dates — in the near-infrared band, where vegetation is bright and bare soil is not:

NIR surface reflectance, same pixel
  2024-06-14:  0.28
  2026-06-09:  0.09
  difference:  -0.19

Differencing a single band is the simplest form. Two refinements are standard. Differencing an index rather than a band — NDVI being the usual choice for vegetation — cancels some of the illumination variation, because a ratio of two bands is less sensitive to a multiplicative brightness change than either band alone. And change vector analysis uses all bands at once:

magnitude = sqrt( sum over bands of (band_2 - band_1)^2 )
direction = the vector of per-band differences

worked, four bands:
  d_blue = -0.01, d_green = +0.02, d_red = +0.06, d_nir = -0.19
  magnitude = sqrt(0.0001 + 0.0004 + 0.0036 + 0.0361)
            = sqrt(0.0402) = 0.2005

The magnitude answers “did this pixel change” and the direction — red up, near-infrared down — answers “in which way”, which here is the signature of vegetation removed. That split is worth preserving: the magnitude is what you threshold, the direction is what you interpret.

Choosing the threshold from the image

A fixed threshold does not transfer between scene pairs. Derive it from the difference image’s own statistics instead. Suppose the NIR difference image across the whole scene has mean 0.002 and standard deviation 0.031:

threshold = mean +/- 2.5 x sd
          = 0.002 +/- 2.5 x 0.031
          = 0.002 +/- 0.0775
          = [-0.0755, +0.0795]

our pixel at -0.19  -> flagged
a pixel at   -0.05  -> not flagged

Now the part that is usually left out. If the unchanged pixels are approximately normal, a two-sided 2.5-sigma cut flags 2 × (1 − Φ(2.5)) = 2 × 0.0062 = 1.24% of them by chance. If genuine change covers 0.3% of the scene — a realistic order of magnitude for a year of a settled landscape — then even with perfect recall:

true positives  = 0.30% of pixels
false positives = 1.24% x 99.7% = 1.24% of pixels
precision       = 0.30 / (0.30 + 1.24) = 19.5%

Four out of five flagged pixels are noise, at a threshold that sounds conservative. That number is derived from stated assumptions — the normality, the change prevalence, the sigma multiple — and every one of them is yours to change, but the shape of the conclusion does not: change detection at the pixel level is inherently precision-poor because change is rare and noise is not.

What actually rescues it is spatial and temporal aggregation. Require a minimum connected-component size — genuine land-cover change is contiguous, sensor noise is not — and a five-pixel minimum removes nearly all of the isolated false positives while costing you only small real changes. Or require the change to persist across a third date, which turns a single noisy comparison into a much stronger test. Otsu’s method on the magnitude histogram is the usual automatic alternative to a sigma rule, and it works well only when the histogram is genuinely bimodal, which it is not when change is rare.

Why post-classification compounds error

The other approach is to classify each scene independently and compare the labels. It has one large advantage — it tells you what changed into what, which differencing never does — and one arithmetic problem.

classifier accuracy per scene: 90%
probability both scenes correct (assuming independence):
  0.90 x 0.90 = 0.81

so up to 19% of pixels can show a label change
that did not happen

Nineteen percent spurious change against a real change rate that might be under one percent. The independence assumption is doing work there and it is pessimistic: two classifications of the same place share error sources, so errors correlate and the real figure is lower. But it is not close to zero, and the errors that do correlate are exactly the ones on hard pixels — mixed pixels at field boundaries, shadowed slopes — which are also where you would most like to trust the answer.

The practical resolution is to use both, in order: difference the scenes to find where change occurred, which is a sensitive test with no classifier in the loop, then classify only within those regions to determine what the change was. That confines the classifier error to the small fraction of the scene it was needed for, and it is much cheaper than classifying two full scenes.

One precondition underlies all of this and belongs before every step: both scenes must be cloud-free at every pixel you compare, because a cloud arriving between two dates is the largest difference in the image and looks exactly like new construction in the visible bands. See cloud masking before analysis.