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.

2025-09-22 usable 40.5%
18 candidate dates; those over 5% clear:
2024-09-03 usable 28.9%, both dates 11.1%
2024-09-13 usable 10.8%, both dates  7.0%
2024-09-18 usable 93.0%, both dates 40.1%
2024-09-20 usable 49.3%, both dates 31.8%
2024-10-03 usable 32.7%, both dates 16.0%
2024-10-05 usable 43.1%, both dates 14.1%

Intersect masks, difference an index, check registration

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 can be much smaller than either mask alone. Over a 10 km window on Snowdon, this compares the 22 September 2025 scene with every Sentinel-2 date in the same month a year earlier:

import os, requests, numpy as np, rasterio
from rasterio.windows import from_bounds

os.environ.update(AWS_NO_SIGN_REQUEST="YES", GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR")
STAC = "https://earth-search.aws.element84.com/v1"
SNOWDON = (422880, 5875440, 432880, 5885440)          # EPSG:32630, tile 30UVD

def clear_mask(item):
    with rasterio.open(item["assets"]["scl"]["href"]) as src:
        scl = src.read(1, window=from_bounds(*SNOWDON, transform=src.transform))
    return np.isin(scl, [4, 5, 6])                    # vegetation, bare soil, water

later = requests.get(f"{STAC}/collections/sentinel-2-l2a/items/S2A_30UVD_20250922_0_L2A",
                     timeout=60).json()
b = clear_mask(later)
print(f"2025-09-22 usable {b.mean():.1%}")
found = requests.post(f"{STAC}/search", timeout=120, json={
    "collections": ["sentinel-2-l2a"], "bbox": [-4.13, 53.02, -3.98, 53.11],
    "datetime": "2024-09-01T00:00:00Z/2024-10-15T23:59:59Z", "limit": 50}).json()["features"]
candidates = sorted((f for f in found if "_30UVD_" in f["id"]),
                    key=lambda f: f["properties"]["datetime"])
print(f"{len(candidates)} candidate dates; those over 5% clear:")
for item in candidates:
    a = clear_mask(item)
    if a.mean() > 0.05:
        print(f"{item['properties']['datetime'][:10]} usable {a.mean():5.1%}, "
              f"both dates {(a & b).mean():5.1%}")

2025-09-22 usable 40.5%
18 candidate dates; those over 5% clear:
2024-09-03 usable 28.9%, both dates 11.1%
2024-09-13 usable 10.8%, both dates  7.0%
2024-09-18 usable 93.0%, both dates 40.1%
2024-09-20 usable 49.3%, both dates 31.8%
2024-10-03 usable 32.7%, both dates 16.0%
2024-10-05 usable 43.1%, both dates 14.1%

Of the 18 candidate dates, most are too cloudy to use at all. Among the rest, 5 October 2024 is 43.1% clear and the 2025 scene is 40.5% clear, yet only 14.1% of cells are clear on both. Report the intersection: a change map covering a seventh 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.

The MAD resists the tails, but not without limit. In the fire example below the burn covers 28% of the window; the MAD of the whole window gives a robust sd of 0.085, against 0.049 for the cells outside the perimeter. Where you know some cells did not change — outside a mapped event, or on land cover that cannot have changed — estimate the noise from those.

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. Over the Airport Fire below, an October-to-October pair flagged 80.4% of the burned area; April to October of the same year flagged 1.3%, because spring green-up had already moved every cell and the threshold rose to match.

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

Change from the Airport Fire, with diagnostics

The data used on this page

The Airport Fire burned 23,519 acres of chaparral in the Santa Ana Mountains, Orange County, California, between 9 September and 6 October 2024. Its final perimeter comes from the WFIGS interagency perimeters service. The imagery is Sentinel-2 L2A from Element 84's Earth Search, tile 11SMT: 20 October 2023 and 19 October 2024 as the anniversary pair, and 2 April 2024 as the out-of-season date. The perimeter gives the examples something no two-date comparison usually has — an independent answer.

import os
import geopandas as gpd
import numpy as np
import rasterio
import requests
from rasterio.enums import Resampling
from rasterio.features import rasterize
from rasterio.windows import from_bounds

os.environ.update(AWS_NO_SIGN_REQUEST="YES", GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR")
WFIGS = ("https://services3.arcgis.com/T4QMspbfLg3qTGWY/arcgis/rest/services/"
         "WFIGS_Interagency_Perimeters/FeatureServer/0/query")
ITEMS = "https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/"

if not os.path.exists("airport_perimeter.geojson"):
    reply = requests.get(WFIGS, timeout=120, params={
        "where": "poly_IncidentName='Airport' AND poly_GISAcres>20000",
        "outFields": "poly_IncidentName,poly_GISAcres", "outSR": 4326, "f": "geojson"})
    open("airport_perimeter.geojson", "w").write(reply.text)
perimeter = gpd.read_file("airport_perimeter.geojson").to_crs(32611)
x0, y0, x1, y1 = perimeter.total_bounds
BOUNDS = (x0 // 20 * 20 - 2000, y0 // 20 * 20 - 2000, -(-x1 // 20) * 20 + 2000, -(-y1 // 20) * 20 + 2000)

def load_ndvi(item_id):
    """NDVI on the 10 m grid and a clear mask from SCL classes 4-6."""
    assets = requests.get(ITEMS + item_id, timeout=60).json()["assets"]
    bands = {}
    for key in ("red", "nir"):
        with rasterio.open(assets[key]["href"]) as src:
            window = from_bounds(*BOUNDS, transform=src.transform)
            bands[key] = src.read(1, window=window).astype("float32") * 1e-4
            transform = src.window_transform(window)
    with rasterio.open(assets["scl"]["href"]) as src:
        scl = src.read(1, window=from_bounds(*BOUNDS, transform=src.transform),
                       out_shape=bands["red"].shape, resampling=Resampling.nearest)
    ndvi = (bands["nir"] - bands["red"]) / (bands["nir"] + bands["red"])
    return ndvi, np.isin(scl, [4, 5, 6]), transform

ndvi_2023, clear_2023, transform = load_ndvi("S2A_11SMT_20231020_0_L2A")
ndvi_2024, clear_2024, _ = load_ndvi("S2B_11SMT_20241019_0_L2A")
ndvi_spring, clear_spring, _ = load_ndvi("S2B_11SMT_20240402_0_L2A")
inside = rasterize(perimeter.geometry, out_shape=ndvi_2023.shape, transform=transform,
                   fill=0, dtype="uint8").astype(bool)
print(f"{ndvi_2023.shape[1]} x {ndvi_2023.shape[0]} cells at 10 m; "
      f"{inside.mean():.1%} inside the perimeter")

def score(change, name):
    """How much of the burn was found, and how much of the change map is the burn."""
    print(f"{name}: {(change & inside).sum() / inside.sum():.1%} of the burned area flagged; "
          f"{(change & inside).sum() / change.sum():.1%} of flagged cells inside the perimeter")

2104 x 1632 cells at 10 m; 27.7% inside the perimeter

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, reference=None):
    """Robust thresholded difference, with the usable area reported.

    reference: optional boolean mask of cells known not to have changed;
    the noise is then estimated from those cells only.
    """
    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)
    sample = difference if reference is None else np.where(reference, difference, np.nan)
    values = sample[np.isfinite(sample)]

    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

change, difference, threshold = detect_change(ndvi_2023, ndvi_2024, clear_2023, clear_2024)
score(change, "anniversary pair")
change_s, _, _ = detect_change(ndvi_spring, ndvi_2024, clear_spring, clear_2024)
score(change_s, "April to October")
change_r, _, _ = detect_change(ndvi_2023, ndvi_2024, clear_2023, clear_2024, reference=~inside)
score(change_r, "anniversary, noise from outside")

  date A usable 100.0%, date B 99.8%, both 99.8%
  difference median -0.0476
  robust sd 0.0847, ordinary sd 0.2036 (2.40x)
  threshold ±0.2118
  removed 8,713 isolated cells in patches under 5 cells
  change 797,200 cells (23.22%): 3,335 increase, 793,865 decrease
anniversary pair: 80.4% of the burned area flagged; 95.9% of flagged cells inside the perimeter
  date A usable 99.1%, date B 99.8%, both 98.9%
  difference median -0.1425
  robust sd 0.2087, ordinary sd 0.2067 (0.99x)
  threshold ±0.5217
  removed 6,749 isolated cells in patches under 5 cells
  change 20,509 cells (0.60%): 5,696 increase, 14,813 decrease
April to October: 1.3% of the burned area flagged; 61.0% of flagged cells inside the perimeter
  date A usable 100.0%, date B 99.8%, both 99.8%
  difference median -0.0235
  robust sd 0.0488, ordinary sd 0.0809 (1.66x)
  threshold ±0.1221
  removed 40,155 isolated cells in patches under 5 cells
  change 967,775 cells (28.18%): 9,639 increase, 958,136 decrease
anniversary, noise from outside: 93.1% of the burned area flagged; 91.5% of flagged cells inside the perimeter

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. The anniversary pair shows 2.40. The April-to-October pair shows 0.99 — not because nothing changed, but because green-up moved the whole scene at once, so there are no tails to find. The signs are elsewhere: a median of −0.14 and a robust sd two and a half times the anniversary pair's.

Estimating the noise from outside the perimeter found 93.1% of the burn instead of 80.4%, and rescued the out-of-season pair from 1.3% to 73.8% when run the same way. The price is more false positives: 91.5% of the flagged cells fall inside the perimeter, against 95.9%. Whether that trade is worth it depends on which mistake costs more.

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 = (int(peak[0]) - max_shift, int(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

check_registration(ndvi_2023, ndvi_2024)
shifted = np.roll(ndvi_2024, 1, axis=1)             # move the later date 10 m east
check_registration(ndvi_2023, shifted)
change_m, _, _ = detect_change(ndvi_2023, shifted, clear_2023, np.roll(clear_2024, 1, axis=1))
score(change_m, "one pixel out")

  peak correlation at offset (0, 0) pixels
  peak correlation at offset (0, -1) pixels
  ! the dates are misregistered — every boundary will show as change
  date A usable 100.0%, date B 99.8%, both 99.8%
  difference median -0.0604
  robust sd 0.1279, ordinary sd 0.2133 (1.67x)
  threshold ±0.3198
  removed 25,578 isolated cells in patches under 5 cells
  change 568,586 cells (16.56%): 3,553 increase, 565,033 decrease
one pixel out: 57.4% of the burned area flagged; 96.0% of flagged cells inside the perimeter

The two dates line up, as scenes from the same relative orbit usually do. Shifting one by a single pixel shows what misregistration does to a robust threshold: the paired edges widen the stable population, the robust sd rises, and the threshold rises with it. The burn is still there; a quarter less of it is found.

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].astype("float64")   # polyfit on float32 is ill-conditioned
    b_values = index_b[stable].astype("float64")
    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())}

corrected, fit = normalise_to_stable(ndvi_2023, ndvi_2024, clear_2023 & clear_2024)
change_n, _, _ = detect_change(ndvi_2023, corrected, clear_2023, clear_2024)
score(change_n, "after normalisation")

  59.9% of cells treated as stable
  b ≈ 0.9527 * a + -0.0083
  date A usable 100.0%, date B 99.8%, both 99.8%
  difference median -0.0143
  robust sd 0.0830, ordinary sd 0.2101 (2.53x)
  threshold ±0.2075
  removed 8,564 isolated cells in patches under 5 cells
  change 827,624 cells (24.10%): 3,371 increase, 824,253 decrease
after normalisation: 83.3% of the burned area flagged; 95.7% of flagged cells inside the perimeter

The cast to float64 is not tidiness. With two million float32 values np.polyfit warns that the fit is poorly conditioned and returns b ≈ 0.4688·a + 0.2591; correcting with that line raised the robust sd of the difference from 0.085 to 0.401, and the change map built from it found 0.4% of the burn. The same cells in float64 give a slope of 0.953.

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.

Why index differences beat band differences

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 resists the tails, so the threshold mostly reflects the noise among the unchanged cells — which is exactly what it should be measuring. "Mostly" matters: its breakdown point is half the sample, and well before that a large change shifts the median and widens the MAD. With the burn at 28% of the window, it gave a threshold 1.7 times the one estimated from outside the perimeter.

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. Over the chaparral outside the Airport Fire, a one-pixel shift doubled the mean absolute difference on the steepest-gradient tenth of cells, from 0.065 to 0.134, and moved it only from 0.051 to 0.062 on the flattest half.

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.

Robust thresholds and co-registration

• 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.

Internal links

• How to calculate NDVI from Sentinel-2 in Python — the index being differenced
• Cloud masking explained — why the masks matter most
• How to mask clouds in Sentinel-2 imagery with Python — building the masks
• Spectral indices explained — choosing the index for the change
• How to extract a vegetation index time series for a polygon — the stronger design
• Radiometric levels explained — a processing change looks like an event
• How to resample satellite bands to a common grid — co-registration
• How to build a cloud-free composite from many satellite scenes — comparing composites instead
• How to detect change between two Sentinel-1 images — change detection with the log-ratio
• dNBR shows burn severity outside the fire perimeter — when the change is real but not the change you meant

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, or two dates from different seasons. Normalise the dates against stable cells before thresholding, and compare the same season across years.

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.