Change Detection

Change detection is THE core technique in geospatial intelligence. Static features are geography. What CHANGED is intelligence. New construction, vehicle movements, ship activity, vegetation removal, flooding — all are detected by comparing images over time.

Key References

• Singh — “Review article: Digital change detection techniques using remotely-sensed data” (International Journal of Remote Sensing, 1989) — The foundational change detection paper
• Bovolo & Bruzzone — “The Time Variable in Change Detection” (IEEE J. Sel. Top. Appl. Earth Obs., 2017) — Systematic framework for change detection methods
• Copernicus Global Land Service: https://land.copernicus.eu/global/ — Free pan-European biophysical parameter time series
• Earth Engine Data Catalog: https://developers.google.com/earth-engine/datasets/catalog — Comprehensive dataset listing

See Pattern 3: Change Is the Signal.

---

Types of Change and Their Intelligence Value

Change Type	Timescale	Detectable At	Intelligence Meaning
New construction	Weeks-months	10m (Sentinel-2)	Capability investment, permanent posture shift
Earthworks / berms	Days-weeks	10m	Defensive preparation, concealment
Vehicle deployment	Hours-days	<1m (commercial)	Operational readiness, logistics
Vegetation removal	Days-weeks	10m	Site preparation, field clearing
Ship arrivals/departures	Hours	5-10m (SAR), <1m (optical)	Naval operations, trade patterns
Flooding	Hours-days	10m (optical, SAR)	Natural disaster, dam failure
Fire/burn scars	Hours-days	10m-1km	Conflict damage, environmental
Camouflage deployment	Hours	<1m (multispectral)	Concealment effort = something to hide

---

Method 1: Image Differencing

The simplest approach: subtract pixel values between two dates. Where the difference exceeds a threshold, change has occurred.

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import binary_opening, binary_closing, label
 
def image_difference(before, after, threshold=None, auto_threshold_sigma=2):
    """
    Simple image differencing for change detection.
    before, after: 2D numpy arrays (same band, same area, co-registered)
    Returns: difference, change_mask, threshold_used
    """
    diff = after.astype(np.float32) - before.astype(np.float32)
 
    if threshold is None:
        # Automatic threshold: mean +/- N*sigma
        mean = np.nanmean(diff)
        std = np.nanstd(diff)
        threshold = auto_threshold_sigma * std
 
    # Detect significant changes
    increase = diff > threshold    # new bright features
    decrease = diff < -threshold   # removed/darkened features
    change_mask = increase | decrease
 
    return diff, change_mask, threshold
 
 
def clean_change_mask(mask, min_area_pixels=10):
    """
    Morphological cleanup: remove isolated pixels, fill small holes.
    """
    # Close small gaps
    mask = binary_closing(mask, iterations=1)
    # Remove small isolated pixels
    mask = binary_opening(mask, iterations=1)
 
    # Remove connected components smaller than min_area
    labeled, n_features = label(mask)
    for i in range(1, n_features + 1):
        if np.sum(labeled == i) < min_area_pixels:
            mask[labeled == i] = False
 
    return mask

Limitations

• Requires perfect co-registration (sub-pixel alignment)
• Sensitive to radiometric differences (different atmosphere, sun angle)
• Works best for single-band comparison

---

Method 2: NDVI Differencing

Compare vegetation index between dates. More robust than raw band differencing because NDVI normalizes for illumination variations.

import numpy as np
 
def ndvi_change_detection(nir_before, red_before, nir_after, red_after,
                          threshold=0.15):
    """
    Detect vegetation change using NDVI difference.
    Threshold 0.15 is typical for significant change.
    Negative dNDVI = vegetation loss. Positive = vegetation gain.
    """
    ndvi_before = (nir_before - red_before) / (nir_before + red_before + 1e-10)
    ndvi_after = (nir_after - red_after) / (nir_after + red_after + 1e-10)
 
    d_ndvi = ndvi_after - ndvi_before
 
    vegetation_loss = d_ndvi < -threshold   # construction, clearing, fire
    vegetation_gain = d_ndvi > threshold    # regrowth, new planting
 
    return {
        "ndvi_before": ndvi_before,
        "ndvi_after": ndvi_after,
        "d_ndvi": d_ndvi,
        "vegetation_loss": vegetation_loss,
        "vegetation_gain": vegetation_gain,
    }

Intelligence application:

• Vegetation loss at a known facility → new construction, cleared area, earthworks
• Vegetation gain where there was bare ground → possible concealment (planting over disturbed earth)
• Patchy vegetation loss → possible vehicle maneuvering through vegetation

---

Method 3: Change Vector Analysis (CVA)

Uses multiple bands simultaneously. For each pixel, compute the magnitude and direction of change in spectral space.

import numpy as np
 
def change_vector_analysis(bands_before, bands_after, threshold_percentile=95):
    """
    Change Vector Analysis using multiple bands.
    bands_before, bands_after: dicts of {name: 2D_array}
    Returns: magnitude (how much changed), direction (what kind of change)
    """
    # Stack bands into multi-dimensional arrays
    band_names = sorted(bands_before.keys())
    stack_before = np.stack([bands_before[b] for b in band_names], axis=0).astype(np.float32)
    stack_after = np.stack([bands_after[b] for b in band_names], axis=0).astype(np.float32)
 
    # Change vector = difference in each band
    diff = stack_after - stack_before
 
    # Magnitude: Euclidean distance in spectral space
    magnitude = np.sqrt(np.sum(diff**2, axis=0))
 
    # Direction: angle of change vector (first two principal components)
    # For 2-band simplification (NIR, Red):
    if "B08" in bands_before and "B04" in bands_before:
        nir_idx = band_names.index("B08")
        red_idx = band_names.index("B04")
        direction = np.arctan2(diff[nir_idx], diff[red_idx])
    else:
        direction = np.arctan2(diff[0], diff[1])
 
    # Threshold
    thresh = np.nanpercentile(magnitude, threshold_percentile)
    change_mask = magnitude > thresh
 
    return {
        "magnitude": magnitude,
        "direction": direction,
        "change_mask": change_mask,
        "threshold": thresh,
    }

Advantage over single-band methods: CVA detects changes that affect multiple bands simultaneously (e.g., construction changes both NIR and SWIR), and the direction vector tells you WHAT TYPE of change occurred (vegetation→urban has a different spectral direction than vegetation→water).

---

Method 4: Object-Based Change Detection

Instead of comparing pixels, segment the image into objects (regions of similar pixels) and compare object properties.

from scipy.ndimage import label, find_objects
import numpy as np
 
def segment_and_compare(before, after, segment_before=True):
    """
    Simple object-based change detection.
    1. Segment the 'before' image into objects
    2. Compare mean values of each object between before and after
    3. Flag objects with significant change
    """
    # Simple segmentation by thresholding + connected components
    # (In practice, use more sophisticated segmentation: SLIC, watershed, etc.)
    threshold = np.mean(before)
    binary = before > threshold
    labeled, n_objects = label(binary)
 
    changes = []
    for obj_id in range(1, n_objects + 1):
        obj_mask = labeled == obj_id
        area = np.sum(obj_mask)
        if area < 20:  # skip tiny objects
            continue
 
        mean_before = np.mean(before[obj_mask])
        mean_after = np.mean(after[obj_mask])
        pct_change = (mean_after - mean_before) / (mean_before + 1e-10) * 100
 
        if abs(pct_change) > 20:  # significant change threshold
            centroid = np.mean(np.argwhere(obj_mask), axis=0)
            changes.append({
                "id": obj_id,
                "centroid": centroid,
                "area_pixels": area,
                "pct_change": pct_change,
                "direction": "increase" if pct_change > 0 else "decrease",
            })
 
    return changes

---

Method 5: Time Series Analysis

For monitoring over many dates, analyze the full time series to detect trends and break points.

import numpy as np
import matplotlib.pyplot as plt
 
def detect_breaks_in_timeseries(dates, values, min_change=0.15, window=3):
    """
    Simple break detection in a time series of NDVI values.
    Detects sudden drops that exceed min_change threshold.
 
    For production use, consider BFAST or LandTrendr algorithms.
    """
    values = np.array(values, dtype=float)
    dates = np.array(dates)
 
    # Moving average as baseline
    baseline = np.convolve(values, np.ones(window)/window, mode="same")
 
    # Deviation from baseline
    deviation = values - baseline
 
    # Detect significant negative breaks
    breaks = []
    for i in range(window, len(values) - 1):
        if deviation[i] < -min_change and deviation[i-1] >= -min_change:
            breaks.append({
                "date": dates[i],
                "value": values[i],
                "baseline": baseline[i],
                "drop": deviation[i],
            })
 
    return breaks, baseline
 
 
# Simulate NDVI time series for a location
np.random.seed(42)
n_dates = 36  # 3 years, monthly
dates = np.arange(n_dates)
months = dates % 12
 
# Seasonal NDVI pattern (peaks in summer, low in winter)
seasonal = 0.3 + 0.25 * np.sin(2 * np.pi * (months - 3) / 12)
 
# Add noise
noise = np.random.normal(0, 0.03, n_dates)
 
# Add a change event at month 18: construction cleared vegetation
values = seasonal + noise
values[18:] -= 0.25  # permanent NDVI drop from construction
values = np.clip(values, -0.1, 0.9)
 
# Detect
breaks, baseline = detect_breaks_in_timeseries(dates, values, min_change=0.12)
 
# Plot
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(dates, values, "ko-", markersize=4, label="NDVI observations")
ax.plot(dates, baseline, "b--", linewidth=1, label="Moving average baseline")
 
for b in breaks:
    ax.axvline(b["date"], color="red", linestyle="--", alpha=0.7)
    ax.annotate(f"Break detected\nNDVI drop: {b['drop']:.2f}",
                xy=(b["date"], b["value"]),
                xytext=(b["date"]+2, b["value"]+0.15),
                arrowprops=dict(arrowstyle="->", color="red"),
                fontsize=9, color="red")
 
# Mark seasons
for y in range(3):
    ax.axvspan(y*12+5, y*12+8, alpha=0.05, color="green")  # summer
 
ax.set_xlabel("Month")
ax.set_ylabel("NDVI")
ax.set_title("NDVI Time Series — Construction Event Detection")
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xticks(range(0, 36, 6))
ax.set_xticklabels([f"M{m}" for m in range(0, 36, 6)])
plt.tight_layout()
plt.savefig("ndvi_timeseries_break.png", dpi=150)
plt.show()

---

Critical Challenge: Radiometric Normalization

Two images of the same location on different dates will have different pixel values even if NOTHING changed on the ground, because of:

• Different atmospheric conditions (haze, aerosols, water vapor)
• Different sun angle (time of year, time of day)
• Different sensor calibration
• Different viewing geometry (off-nadir angle)

Solutions:

1. Use L2A products (surface reflectance) — atmospheric correction already applied
2. Use indices (NDVI, NDWI) — ratio-based, partially self-normalizing
3. Relative normalization: find stable pseudo-invariant features (PIF) in both images, compute linear transform to match
4. Histogram matching: adjust after image histogram to match before image

import numpy as np
 
def relative_radiometric_normalization(reference, target, n_samples=10000):
    """
    Pseudo-Invariant Feature normalization.
    Assumes most pixels didn't change (true for most areas).
    Fits a linear model between reference and target pixel values.
    """
    # Sample random pixels
    valid = ~(np.isnan(reference) | np.isnan(target))
    valid_idx = np.argwhere(valid)
    if len(valid_idx) > n_samples:
        sample_idx = valid_idx[np.random.choice(len(valid_idx), n_samples, replace=False)]
    else:
        sample_idx = valid_idx
 
    ref_vals = reference[sample_idx[:, 0], sample_idx[:, 1]]
    tgt_vals = target[sample_idx[:, 0], sample_idx[:, 1]]
 
    # Fit linear model: target_corrected = a * target + b
    # such that corrected values match reference statistics
    a = np.std(ref_vals) / (np.std(tgt_vals) + 1e-10)
    b = np.mean(ref_vals) - a * np.mean(tgt_vals)
 
    corrected = a * target + b
    return corrected

---

Critical Challenge: Cloud Masking

Clouds are the enemy of optical change detection. A cloud in one image looks like a massive change.

import numpy as np
 
def apply_cloud_mask(data, scl_band, valid_classes=None):
    """
    Mask clouds using Sentinel-2 Scene Classification Layer.
    scl_band: SCL from the Sentinel-2 product (20m, resample to match data)
 
    SCL values: 0=NoData, 1=Saturated, 2=Dark, 3=CloudShadow,
    4=Vegetation, 5=Bare, 6=Water, 7=Unclassified,
    8=CloudMedium, 9=CloudHigh, 10=Cirrus, 11=Snow
    """
    if valid_classes is None:
        valid_classes = [4, 5, 6, 7, 11]  # veg, bare, water, unclass, snow
 
    cloud_free = np.isin(scl_band.astype(int), valid_classes)
    masked = np.where(cloud_free, data, np.nan)
    return masked, cloud_free

---

Complete Change Detection Pipeline

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import binary_opening, binary_closing, label
 
def full_change_detection_pipeline(
    nir_before, red_before, scl_before,
    nir_after, red_after, scl_after,
    ndvi_threshold=0.15, min_area_m2=1000, pixel_size_m=10,
):
    """
    Complete change detection pipeline for Sentinel-2 data.
 
    Steps:
    1. Cloud mask both dates
    2. Compute NDVI for both dates
    3. Difference NDVI
    4. Threshold
    5. Morphological cleanup
    6. Compute change statistics
 
    Returns dict with all intermediate and final products.
    """
    # Step 1: Cloud mask
    valid_before = np.isin(scl_before.astype(int), [4, 5, 6, 7, 11])
    valid_after = np.isin(scl_after.astype(int), [4, 5, 6, 7, 11])
    both_valid = valid_before & valid_after
 
    # Step 2: Compute NDVI
    ndvi_before = np.where(
        both_valid,
        (nir_before - red_before) / (nir_before + red_before + 1e-10),
        np.nan,
    )
    ndvi_after = np.where(
        both_valid,
        (nir_after - red_after) / (nir_after + red_after + 1e-10),
        np.nan,
    )
 
    # Step 3: Difference
    d_ndvi = ndvi_after - ndvi_before
 
    # Step 4: Threshold
    loss_mask = d_ndvi < -ndvi_threshold
    gain_mask = d_ndvi > ndvi_threshold
 
    # Step 5: Cleanup
    min_pixels = int(min_area_m2 / (pixel_size_m ** 2))
    loss_clean = _cleanup_mask(loss_mask, min_pixels)
    gain_clean = _cleanup_mask(gain_mask, min_pixels)
 
    # Step 6: Statistics
    pixel_area_m2 = pixel_size_m ** 2
    loss_area = np.sum(loss_clean) * pixel_area_m2
    gain_area = np.sum(gain_clean) * pixel_area_m2
 
    # Label individual change regions
    loss_labeled, n_loss = label(loss_clean)
    gain_labeled, n_gain = label(gain_clean)
 
    loss_regions = []
    for i in range(1, n_loss + 1):
        region_mask = loss_labeled == i
        area = np.sum(region_mask) * pixel_area_m2
        centroid = np.mean(np.argwhere(region_mask), axis=0)
        mean_change = np.nanmean(d_ndvi[region_mask])
        loss_regions.append({
            "id": i,
            "area_m2": area,
            "centroid_rowcol": centroid,
            "mean_dndvi": mean_change,
        })
 
    return {
        "ndvi_before": ndvi_before,
        "ndvi_after": ndvi_after,
        "d_ndvi": d_ndvi,
        "loss_mask": loss_clean,
        "gain_mask": gain_clean,
        "loss_area_m2": loss_area,
        "gain_area_m2": gain_area,
        "n_loss_regions": n_loss,
        "n_gain_regions": n_gain,
        "loss_regions": sorted(loss_regions, key=lambda r: -r["area_m2"]),
        "valid_mask": both_valid,
    }
 
 
def _cleanup_mask(mask, min_pixels):
    """Morphological cleanup of binary mask."""
    mask = mask.copy()
    mask = binary_closing(mask, iterations=1)
    mask = binary_opening(mask, iterations=1)
    labeled, n = label(mask)
    for i in range(1, n + 1):
        if np.sum(labeled == i) < min_pixels:
            mask[labeled == i] = False
    return mask
 
 
def plot_change_detection(results, save_path=None):
    """Visualize change detection results."""
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
 
    # NDVI before
    im = axes[0, 0].imshow(results["ndvi_before"], cmap="RdYlGn", vmin=-0.1, vmax=0.8)
    axes[0, 0].set_title("NDVI Before")
    plt.colorbar(im, ax=axes[0, 0], fraction=0.046)
 
    # NDVI after
    im = axes[0, 1].imshow(results["ndvi_after"], cmap="RdYlGn", vmin=-0.1, vmax=0.8)
    axes[0, 1].set_title("NDVI After")
    plt.colorbar(im, ax=axes[0, 1], fraction=0.046)
 
    # Difference
    im = axes[0, 2].imshow(results["d_ndvi"], cmap="RdBu", vmin=-0.5, vmax=0.5)
    axes[0, 2].set_title("NDVI Difference (after - before)")
    plt.colorbar(im, ax=axes[0, 2], fraction=0.046)
 
    # Valid mask
    axes[1, 0].imshow(results["valid_mask"], cmap="RdYlGn")
    axes[1, 0].set_title("Valid Pixels (cloud-free both dates)")
 
    # Loss mask
    axes[1, 1].imshow(results["loss_mask"], cmap="Reds")
    axes[1, 1].set_title(f"Vegetation Loss ({results['n_loss_regions']} regions, "
                          f"{results['loss_area_m2']/1e4:.1f} ha)")
 
    # Combined overlay
    overlay = np.zeros((*results["d_ndvi"].shape, 3))
    overlay[results["loss_mask"]] = [1, 0, 0]    # red = loss
    overlay[results["gain_mask"]] = [0, 1, 0]    # green = gain
    overlay[~results["valid_mask"]] = [0.5, 0.5, 0.5]  # gray = clouds
    axes[1, 2].imshow(overlay)
    axes[1, 2].set_title("Change Overlay (red=loss, green=gain)")
 
    for ax in axes.flat:
        ax.axis("off")
 
    plt.tight_layout()
    if save_path:
        plt.savefig(save_path, dpi=150, bbox_inches="tight")
    plt.show()
 
    # Print top change regions
    print("\nTop vegetation loss regions:")
    for r in results["loss_regions"][:5]:
        print(f"  Region {r['id']}: {r['area_m2']:.0f} m^2, "
              f"mean dNDVI={r['mean_dndvi']:.3f}, "
              f"center=({r['centroid_rowcol'][0]:.0f}, {r['centroid_rowcol'][1]:.0f})")

Synthetic Demo

import numpy as np
 
np.random.seed(42)
size = 300
 
# Before: mixed land cover
nir_before = np.random.normal(0.45, 0.05, (size, size))  # vegetated
red_before = np.random.normal(0.05, 0.01, (size, size))
# Existing built-up area
nir_before[120:160, 80:140] = np.random.normal(0.20, 0.03, (40, 60))
red_before[120:160, 80:140] = np.random.normal(0.12, 0.02, (40, 60))
# Water
nir_before[220:260, 50:120] = np.random.normal(0.01, 0.005, (40, 70))
red_before[220:260, 50:120] = np.random.normal(0.03, 0.005, (40, 70))
 
scl_before = np.full((size, size), 4, dtype=np.uint8)  # all valid vegetation
scl_before[220:260, 50:120] = 6  # water
 
# After: construction event + some vegetation loss
nir_after = nir_before.copy()
red_after = red_before.copy()
 
# New construction: vegetation cleared, replaced with bare/built
nir_after[40:70, 180:240] = np.random.normal(0.18, 0.03, (30, 60))
red_after[40:70, 180:240] = np.random.normal(0.14, 0.02, (30, 60))
 
# Road construction
nir_after[150:155, 140:250] = np.random.normal(0.15, 0.02, (5, 110))
red_after[150:155, 140:250] = np.random.normal(0.13, 0.02, (5, 110))
 
# Cloud in after image (small patch)
scl_after = scl_before.copy()
scl_after[0:30, 0:50] = 9  # cloud
 
# Run pipeline
results = full_change_detection_pipeline(
    nir_before.astype(np.float32), red_before.astype(np.float32), scl_before,
    nir_after.astype(np.float32), red_after.astype(np.float32), scl_after,
    ndvi_threshold=0.15,
    min_area_m2=500,
    pixel_size_m=10,
)
 
plot_change_detection(results, save_path="change_detection_pipeline.png")

---

Connection to Other Techniques

• SAR change detection: Uses amplitude ratio instead of differencing. Advantage: works through clouds. See SAR Fundamentals and Analysis.
• Object detection on changes: Run object detection only on changed areas (tip and cue).
• Forensics: Change detection is evidence of activity over time. Document with timestamps, coordinates, and before/after pairs.
• OSINT corroboration: Always correlate imagery changes with open-source intelligence. See Multi-Source Intelligence Fusion.

---

Exercises

Exercise 1: Detect Construction at a Facility

1. Download two Sentinel-2 scenes of a location where you know construction has occurred (e.g., a new neighborhood being built)
2. Choose scenes 6-12 months apart, with <20% cloud cover
3. Run the full pipeline above
4. Identify and describe each change region

Exercise 2: Map Flood Extent

1. Find a flood event (search news for recent flooding)
2. Download Sentinel-2 before and during/after the flood
3. Use NDWI differencing instead of NDVI: areas that became water
4. Compute flooded area in km
5. Bonus: use Sentinel-1 SAR for the same flood (SAR sees through clouds that accompany floods)

Exercise 3: Track Deforestation Over 12 Months

1. Download 12 Sentinel-2 scenes (one per month) for a tropical area with deforestation (Rondonia, Brazil works well with Landsat)
2. Compute NDVI for each date
3. Create a time series for 5 sample points
4. Apply the break detection function to find when deforestation occurred
5. Create a “date of change” map

---

Self-Test Questions

1. Why is NDVI differencing more robust than raw band differencing?
2. You detect a large NDVI decrease in winter between two years. Is this necessarily real change? What else could cause it?
3. What is the minimum detectable change size at 10m resolution?
4. Why must you cloud-mask BOTH images before computing change?
5. A facility shows no change in optical imagery but SAR coherence dropped. What might have happened?

---

See also: Multispectral Analysis | SAR Fundamentals and Analysis | Case Study - Monitoring Military Installations Next: SAR Fundamentals and Analysis