Separating Real Change from Illumination Differences

Difference normalised indices, not raw bands, and normalise the two dates on pixels that did not change:

import numpy as np

def pif_normalise(target, reference, pif_mask):
    """Fit target -> reference on pseudo-invariant pixels and apply to all pixels."""
    a, b = np.polyfit(target[pif_mask], reference[pif_mask], 1)
    return a * target + b

nir_t2 = pif_normalise(nir_t2, nir_t1, pif)
red_t2 = pif_normalise(red_t2, red_t1, pif)
dndvi = (nir_t2 - red_t2) / (nir_t2 + red_t2) - (nir_t1 - red_t1) / (nir_t1 + red_t1)

A difference image records every reason two dates differ — land change, but also a lower sun, longer shadows, a hazier atmosphere and a different view angle. This page belongs to change detection and differencing workflows in Satellite Processing Workflows & Index Pipelines.

Where False Change Comes From

The cheapest remedy of all is date selection. Comparing anniversary dates — the same week in different years — equalises sun elevation and phenology at once, and removes most of the problem before any processing. Where the question forces dates from different seasons, the remaining remedies become essential rather than optional.

Environment & Setup

pip install "numpy>=1.23" "rasterio>=1.3.0" "scipy>=1.11"

Complete Working Example

import numpy as np
from scipy import stats


def find_pifs(t1: dict, t2: dict, *, quantile: float = 0.1) -> np.ndarray:
    """Pseudo-invariant features: pixels whose NIR/red ratio barely changed."""
    r1 = t1["nir"] / np.maximum(t1["red"], 1e-4)
    r2 = t2["nir"] / np.maximum(t2["red"], 1e-4)
    change = np.abs(np.log(r2 / r1))
    valid = np.isfinite(change)
    thresh = np.nanquantile(change[valid], quantile)
    return valid & (change <= thresh)


def normalise(t1: dict, t2: dict, pif: np.ndarray) -> dict:
    out = {}
    for band in t2:
        x, y = t2[band][pif], t1[band][pif]
        fit = stats.theilslopes(y, x)              # robust to leftover outliers
        out[band] = fit.slope * t2[band] + fit.intercept
    return out


def ndvi(d: dict) -> np.ndarray:
    return (d["nir"] - d["red"]) / (d["nir"] + d["red"])


def illumination_robust_dndvi(t1: dict, t2: dict) -> np.ndarray:
    pif = find_pifs(t1, t2)
    t2n = normalise(t1, t2, pif)
    return ndvi(t2n) - ndvi(t1)

Pseudo-invariant features — rooftops, tarmac, deep water, bare rock — are surfaces whose reflectance should not change between dates. Any difference observed on them is illumination and atmosphere, so a regression fitted on them describes the non-change part of the relationship and can be removed. Selecting them automatically from the least-changed decile of a ratio works well; a hand-curated mask of known stable surfaces works better where one exists.

Why Ratios Help and Where They Stop Helping

A normalised difference cancels any factor that multiplies both bands equally, which is what a change in sun elevation mostly is. It does not cancel additive effects: haze adds path radiance to both bands, which shifts the ratio, and at low sun angles canopy self-shadowing changes red and NIR differently. That is why ratio indices on their own are a good first defence but not a complete one, and why relative normalisation on invariant pixels remains worthwhile.

Terrain Makes It Local

On flat land, illumination differences are nearly uniform and a single regression per band removes them. In mountains the effect depends on slope and aspect: a north-facing slope that was lit in June is in shadow in December, and no global correction can fix that. Applying a topographic correction to each date first, as described in applying a terrain illumination correction, and then masking pixels in cast shadow on either date, is the only reliable approach. Where both are impractical, at least stratify the change threshold by illumination class so that steep shaded slopes need a larger difference before they are flagged.

Setting a Change Threshold that Respects Noise

After normalisation, the difference image over stable land should be centred on zero with a spread that reflects residual noise. That spread is the natural yardstick for a threshold: estimate its standard deviation on the invariant pixels, or on a robust measure such as the median absolute deviation over the whole scene, and flag change only where the difference exceeds two to three times that value. Because the spread depends on the dates, a threshold derived this way adapts automatically, where a fixed number tuned on one pair will be too strict for a clean pair and too loose for a hazy one. The follow-on steps of cleaning the binary map are covered in thresholding change maps and removing noise.

Keeping the Evidence

Save the invariant-pixel mask and the fitted coefficients for every date pair next to the change map. When someone questions a flagged area months later, those few numbers show at once whether the normalisation was sound.

Verification

import numpy as np

pif = find_pifs(t1, t2)
before = ndvi(t2) - ndvi(t1)
after = illumination_robust_dndvi(t1, t2)
print("before: mean", np.nanmean(before[pif]), "std", np.nanstd(before[pif]))
print("after:  mean", np.nanmean(after[pif]), "std", np.nanstd(after[pif]))
assert abs(np.nanmean(after[pif])) < 0.01

Also check a hand-picked set of known stable sites that were not used to fit the normalisation; if they show a residual shift, the automatic invariant set was contaminated by real change.

Common Errors

Whole scene appears to have changed

The dates differ in sun angle or haze and were differenced raw. Use surface reflectance, ratio indices and PIF normalisation.

North-facing slopes flagged as vegetation loss

Seasonal terrain shadow. Correct topography and mask cast shadow on both dates.

Normalisation makes things worse

The invariant set includes changed pixels. Tighten the quantile or use a curated stable-surface mask.

The fit is dominated by a few bright pixels

Ordinary least squares is sensitive to outliers. Use Theil–Sen or another robust regression.

Frequently Asked Questions

Q: What are pseudo-invariant features? Surfaces whose reflectance should not change between dates — roofs, roads, deep water, bare rock. Any difference measured on them is due to illumination and atmosphere, so they are used to fit and remove that difference.

Q: Do normalised indices remove illumination effects? They cancel effects that multiply both bands equally, which covers much of the sun-angle difference. They do not cancel additive haze or terrain-dependent shadowing, which need atmospheric correction and topographic correction.

Q: Why compare anniversary dates? Because the same week in different years has nearly the same sun elevation and phenology, removing the two largest sources of false change before any processing.

Q: Is surface reflectance enough on its own? It removes most atmospheric difference, but residual errors between dates remain. Relative normalisation on invariant pixels removes what is left.

Related

• Change Detection and Differencing Workflows — the parent topic.
• Computing NDVI Difference between Two Dates — the basic differencing step.
• Harmonising Landsat and Sentinel-2 Reflectance — cross-sensor consistency before differencing.
• Thresholding Change Maps and Removing Noise — turning the difference into a change map.