Detecting Distribution Drift in Spatial Datasets

A building-footprint segmenter trained on cloud-free summer imagery ships to production and scores well on the held-out test set. Six months later the same model starts missing rooftops — not because the weights changed, but because the incoming tiles did. Winter acquisitions carry snow cover, low sun angles, and long shadows; a newly onboarded sensor delivers a slightly different spectral response and a coarser ground sample distance. None of this raises an exception. The pipeline keeps producing predictions, the dashboards stay green, and accuracy quietly decays until a downstream consumer notices the maps are wrong. By then weeks of degraded inference have already propagated into derived products.

This is distribution drift, and it is the failure mode that a monitoring layer inside an active learning loop exists to catch. The goal is to detect when the imagery or the label distribution feeding your model has moved away from what it was trained on, quantify how far it has moved, and raise a trigger that routes the divergent tiles back to annotators before model performance collapses. This guide builds that monitor end to end: per-band spectral histograms, Population Stability Index and Kolmogorov-Smirnov covariate scoring, class-prior comparison, and a banded threshold gate that emits a re-labeling signal.

# Prerequisites & Monitoring Contract

A drift monitor is only trustworthy if its reference is fixed and its inputs are normalized. Before computing a single statistic, agree on what “the distribution the model was trained on” means and freeze it. That reference is the frozen training set the deployed checkpoint was fit on — not last week’s data, not a rolling window that silently absorbs the very drift you want to catch.

Required Python packages (pinned):

pip install numpy==1.26.4 pandas==2.2.2 scipy==1.13.1 rasterio==1.3.10 scikit-learn==1.5.1

Monitoring contract to establish first:

• A frozen reference dataset, versioned so every drift report cites the exact snapshot it compared against.
• A canonical coordinate reference system and ground sample distance that every batch is resampled to before statistics are computed.
• Fixed histogram bin edges derived once from the reference, reused for every batch — comparing histograms with different edges is meaningless.
• A minimum sample size per batch below which alerts are held as provisional.
• A destination for the re-label trigger: the review queue that the active learning loop drains.

What drift monitoring assumes about the data:

• Every tile carries a valid CRS and geotransform, so pixel statistics reflect scene content rather than a projection artefact.
• Band ordering is consistent across sensors, or an explicit band-mapping table normalizes it.
• Class labels use a stable label taxonomy so a “class-prior shift” is a real change in scene composition, not a renamed category.

With those fixed, the monitor reduces to three questions asked on every batch: has the imagery changed (covariate drift), has the class mix changed (prior drift), and is the change large enough to act on (thresholding).

# Core Drift-Monitoring Workflow

# Compute per-band spectral histograms

The covariate in a geospatial pipeline is the imagery itself. Rather than track every pixel, summarize each band as a normalized histogram over fixed bin edges. Edges are computed once from the reference and frozen; every subsequent batch is binned against those same edges so the two histograms are directly comparable.

from __future__ import annotations

from pathlib import Path

import numpy as np
import rasterio


def band_edges(reference_paths: list[str], n_bins: int = 20) -> np.ndarray:
    """Derive fixed histogram bin edges from the reference set (per-band shared edges)."""
    samples: list[np.ndarray] = []
    for path in reference_paths:
        with rasterio.open(path) as src:
            arr = src.read().astype(np.float64)  # (bands, h, w)
        samples.append(arr.reshape(arr.shape[0], -1))
    stacked = np.concatenate(samples, axis=1)  # (bands, pixels)
    lo, hi = np.percentile(stacked, [1.0, 99.0])
    return np.linspace(float(lo), float(hi), n_bins + 1)


def band_histograms(path: str, edges: np.ndarray) -> np.ndarray:
    """Return an (n_bands, n_bins) array of normalized per-band histograms."""
    with rasterio.open(path) as src:
        arr = src.read(masked=True).astype(np.float64)
    n_bands = arr.shape[0]
    n_bins = len(edges) - 1
    out = np.zeros((n_bands, n_bins), dtype=np.float64)
    for b in range(n_bands):
        pixels = arr[b].compressed()  # drops nodata via the mask
        counts, _ = np.histogram(pixels, bins=edges)
        total = counts.sum()
        out[b] = counts / total if total else counts
    return out


def batch_histograms(paths: list[str], edges: np.ndarray) -> np.ndarray:
    """Average per-band histograms across every tile in a batch."""
    per_tile = [band_histograms(p, edges) for p in paths]
    return np.mean(np.stack(per_tile, axis=0), axis=0)

Reading with masked=True matters: nodata borders from reprojected or clipped tiles would otherwise pile into an edge bin and manufacture drift where none exists. Averaging per-tile histograms rather than pooling all pixels keeps a single unusually large tile from dominating the batch summary.

# Measure covariate drift with PSI and the KS test

With reference and batch histograms in hand, quantify how far they diverge. Two complementary statistics do the job. The Population Stability Index is a binned, symmetric measure that sums per-bin contributions and maps cleanly onto action bands. The two-sample Kolmogorov-Smirnov test works on the raw pixel samples and returns a distribution-free significance value, catching shape changes that a coarse binning might blur.

from __future__ import annotations

import numpy as np
import rasterio
from scipy import stats


def population_stability_index(
    reference: np.ndarray,
    batch: np.ndarray,
    epsilon: float = 1e-6,
) -> float:
    """PSI between two normalized histograms over identical bin edges."""
    ref = np.clip(reference, epsilon, None)
    cur = np.clip(batch, epsilon, None)
    ref = ref / ref.sum()
    cur = cur / cur.sum()
    return float(np.sum((cur - ref) * np.log(cur / ref)))


def band_psi(reference_hist: np.ndarray, batch_hist: np.ndarray) -> dict[int, float]:
    """PSI per band, given (n_bands, n_bins) histogram matrices."""
    return {
        b: population_stability_index(reference_hist[b], batch_hist[b])
        for b in range(reference_hist.shape[0])
    }


def band_ks(
    reference_paths: list[str],
    batch_paths: list[str],
    band: int,
    max_pixels: int = 200_000,
) -> tuple[float, float]:
    """Two-sample KS statistic and p-value for one band, on subsampled pixels."""
    rng = np.random.default_rng(42)

    def sample(paths: list[str]) -> np.ndarray:
        vals: list[np.ndarray] = []
        for path in paths:
            with rasterio.open(path) as src:
                arr = src.read(band + 1, masked=True).astype(np.float64)
            vals.append(arr.compressed())
        pool = np.concatenate(vals)
        if pool.size > max_pixels:
            pool = rng.choice(pool, size=max_pixels, replace=False)
        return pool

    ref = sample(reference_paths)
    cur = sample(batch_paths)
    result = stats.ks_2samp(ref, cur)
    return float(result.statistic), float(result.pvalue)

PSI answers “how much has this band moved” on a scale you can threshold; the KS p-value answers “is the move statistically real given the sample size.” Reporting both prevents two opposite mistakes: acting on a large PSI that is pure small-batch noise, and ignoring a modest but highly significant shift on a large batch.

# Measure label and prior drift

Covariate drift describes the imagery; prior drift describes the labels. If summer batches were 8% impervious surface and winter batches are 20% because snow reshapes what annotators mark, the model’s learned class priors no longer match the incoming scene composition even when every band histogram looks stable. Track per-class annotation frequency and score it with the same PSI machinery.

from __future__ import annotations

from collections import Counter

import numpy as np
import pandas as pd


def class_frequencies(labels: list[str], classes: list[str]) -> np.ndarray:
    """Normalized frequency vector over a fixed class ordering."""
    counts = Counter(labels)
    total = sum(counts.values())
    if total == 0:
        return np.zeros(len(classes), dtype=np.float64)
    return np.array([counts.get(c, 0) / total for c in classes], dtype=np.float64)


def prior_drift(
    reference_labels: list[str],
    batch_labels: list[str],
    classes: list[str],
) -> pd.DataFrame:
    """Per-class reference vs batch frequency plus a PSI contribution table."""
    ref = class_frequencies(reference_labels, classes)
    cur = class_frequencies(batch_labels, classes)
    eps = 1e-6
    ref_c = np.clip(ref, eps, None)
    cur_c = np.clip(cur, eps, None)
    contribution = (cur_c - ref_c) * np.log(cur_c / ref_c)
    return pd.DataFrame(
        {
            "class": classes,
            "reference_freq": ref,
            "batch_freq": cur,
            "psi_contribution": contribution,
        }
    ).sort_values("psi_contribution", ascending=False, ignore_index=True)

The returned frame does double duty: its psi_contribution column sums to the overall class-prior PSI, and sorting by it surfaces exactly which classes drove the shift. A large contribution from a single rare class often signals a real emerging phenomenon — a new construction site, seasonal flooding — that is precisely what you want annotators to look at.

# Threshold and raise a re-label trigger

The final stage collapses every per-band and per-class score into a single decision. Banded PSI thresholds turn continuous scores into stable, moderate, and significant verdicts; the gate raises a re-label trigger when any covariate band or the class-prior distribution crosses into significant, and opens a softer review ticket at moderate.

from __future__ import annotations

from dataclasses import dataclass


@dataclass(frozen=True)
class DriftBands:
    stable_max: float = 0.10
    moderate_max: float = 0.25


def classify(psi: float, bands: DriftBands = DriftBands()) -> str:
    if psi < bands.stable_max:
        return "stable"
    if psi < bands.moderate_max:
        return "moderate"
    return "significant"


@dataclass(frozen=True)
class DriftReport:
    band_psi: dict[int, float]
    prior_psi: float
    verdict: str
    trigger_relabel: bool


def evaluate_batch(
    band_psi_scores: dict[int, float],
    prior_psi: float,
    min_samples_met: bool,
    bands: DriftBands = DriftBands(),
) -> DriftReport:
    scores = list(band_psi_scores.values()) + [prior_psi]
    verdicts = [classify(s, bands) for s in scores]
    if "significant" in verdicts:
        overall = "significant"
    elif "moderate" in verdicts:
        overall = "moderate"
    else:
        overall = "stable"
    # Small-batch guard: never fire a hard trigger on an undersized batch.
    trigger = overall == "significant" and min_samples_met
    return DriftReport(
        band_psi=band_psi_scores,
        prior_psi=prior_psi,
        verdict=overall,
        trigger_relabel=trigger,
    )

The min_samples_met guard is deliberate: an undersized batch can produce a significant verdict on sampling noise alone, so the gate downgrades it to a provisional review rather than spending annotation budget. When trigger_relabel is True, the report is the payload you hand to the annotation queue.

# Drift Metrics & Threshold Reference

Treat the numeric bands as starting points, not universals. A high-cadence daily-revisit constellation with tight radiometric calibration can tolerate a lower significant threshold; a heterogeneous archive stitched from several sensors may need a higher one to avoid alert fatigue. Calibrate against a labeled drift event if you have one — replay a known summer-to-winter transition and tune the bands so it trips exactly once.

# Edge Cases & Gotchas

# Seasonal cycles masquerading as drift

Vegetation greenness, snow, sun angle, and soil moisture all swing on an annual cycle. Compared against a single summer reference, every winter batch trips the gate — but that is expected variation the model should already handle, not a reason to re-label. Build the reference from a full annual cycle, or maintain season-matched references and compare like against like. True drift accumulates and does not revert across cycles; seasonal signal returns to baseline every year, and a rolling annual reference absorbs it.

# Sensor calibration shifts

Onboarding a new satellite or drone camera, or a vendor re-processing an archive to a new radiometric baseline, moves band histograms wholesale even though the scenes are unchanged. This is genuine covariate drift from the model’s perspective, but the fix is often a radiometric harmonization step rather than re-labeling. Tag every tile with its sensor and processing-baseline identifier, and stratify drift reports by sensor so a calibration shift is diagnosed as such instead of being blamed on scene content.

# Small-batch noise

PSI is unstable on small samples: with a handful of tiles, an empty bin clipped to epsilon can dominate the sum and inflate the score. Enforce a minimum sample size, accumulate tiles until it is met, and mark any alert on an undersized batch as provisional — exactly what the min_samples_met guard in the gate encodes. Bootstrapping a confidence interval on PSI across resampled tiles is a cheap way to see whether a score is robust or an artefact of one outlier tile.

# CRS and resolution mismatch posing as drift

A batch delivered in a different coordinate reference system or at a coarser ground sample distance produces pixel statistics that diverge from the reference for reasons that have nothing to do with the scene. Reprojection resamples pixel values; a resolution change alters texture and edge density; both shift histograms and can trip the covariate gate on a pipeline configuration difference. Normalize CRS and GSD before any histogram is computed, and assert both in the monitor so a mismatch raises a configuration error rather than a false drift alert. A batch whose only “drift” is a projection artefact will waste an entire re-labeling cycle if it slips through.

# Multimodal bands and bin starvation

Bands with genuinely multimodal distributions — water versus land in a near-infrared band, for instance — can leave interior bins near-empty, making PSI hypersensitive there. Derive bin edges from percentiles rather than a linear span, or use quantile bins so each carries comparable mass. Empty reference bins clipped to epsilon are the single most common source of spuriously large PSI.

# Integration & Automation Hooks

# Feed drift alerts into the active learning loop

A drift report is only useful if it changes what gets annotated. When trigger_relabel is True, push the batch identifier and the top-contributing bands and classes into the review queue that the active learning loop consumes. Drift monitoring and uncertainty sampling are complementary signals: uncertainty finds tiles the current model is unsure about, while drift finds tiles that are unlike anything the model has seen. A tile flagged by both is the highest-value annotation target in the batch, so union the two queues and prioritize the intersection.

# Version the reference alongside the dataset with DVC

The reference distribution must move in lockstep with the model. Track it as an artefact under DVC versioning so every drift report cites the exact reference snapshot, and rebuild the reference as a pipeline stage immediately after each retraining. This keeps a subtle bug at bay: comparing new batches against a reference from two model generations ago will report drift that the current model has already absorbed.

# dvc.yaml
stages:
  build_drift_reference:
    cmd: python scripts/build_reference.py --train data/train/ --out artifacts/reference.npz
    deps:
      - scripts/build_reference.py
      - data/train/
    outs:
      - artifacts/reference.npz
  drift_gate:
    cmd: python scripts/drift_gate.py --reference artifacts/reference.npz --batch data/incoming/
    deps:
      - scripts/drift_gate.py
      - artifacts/reference.npz
      - data/incoming/
    metrics:
      - reports/drift.json:
          cache: false

Emitting the drift verdict as a DVC metric means the gate result is diffable across runs, and a significant verdict can fail the pipeline stage so a divergent batch never reaches training unreviewed.

# Validation & Testing

A drift monitor that never fires is as dangerous as one that always fires. Validate it against synthetic ground truth before trusting it in production: construct a batch you know is drifted and confirm the gate trips, and construct an in-distribution batch and confirm it stays quiet.

from __future__ import annotations

import numpy as np


def test_identical_batch_is_stable() -> None:
    rng = np.random.default_rng(0)
    hist = rng.dirichlet(np.ones(20))
    assert population_stability_index(hist, hist) < 1e-9


def test_shifted_batch_trips_significant() -> None:
    edges = np.linspace(0.0, 1.0, 21)
    ref_pixels = np.random.default_rng(1).normal(0.4, 0.05, 50_000)
    cur_pixels = np.random.default_rng(2).normal(0.7, 0.05, 50_000)
    ref_hist, _ = np.histogram(ref_pixels, bins=edges, density=True)
    cur_hist, _ = np.histogram(cur_pixels, bins=edges, density=True)
    psi = population_stability_index(ref_hist, cur_hist)
    assert classify(psi) == "significant"


def test_small_batch_holds_trigger() -> None:
    report = evaluate_batch(
        band_psi_scores={0: 0.9},
        prior_psi=0.02,
        min_samples_met=False,
    )
    assert report.verdict == "significant"
    assert report.trigger_relabel is False

Three properties are worth asserting in CI: a batch identical to the reference scores near-zero PSI, a deliberately shifted batch classifies as significant, and an undersized batch never raises a hard trigger regardless of its score. Add a replay test over a historical drift event if one exists in your archive — it is the closest thing to a real-world integration check the monitor will get.

Beyond unit tests, log every batch’s full drift report even when it stays stable. The stream of scores over time is itself a diagnostic: a slow upward creep in band PSI that never quite crosses 0.25 is early warning that the reference is aging and a retrain is due, long before any single batch trips the gate.

# Frequently Asked Questions

# Related

• Monitoring Class-Balance Drift Across Image Tiles — a focused walkthrough of the class-prior PSI monitor with a full pandas implementation
• Uncertainty Sampling for Geospatial Active Learning — the complementary signal that ranks tiles by model uncertainty rather than distributional distance
• Closing the Loop with Automated Model Retraining — where a confirmed drift trigger ultimately lands: an evaluation-gated retraining pipeline
• Implementing DVC for Geospatial Training Data — versioning the frozen reference so every drift report is reproducible
• Coordinate Reference Systems in Annotation Pipelines — normalizing CRS and resolution so projection artefacts do not masquerade as drift

This guide is part of the broader Active Learning & Model Feedback Loops for Geospatial Annotation topic area that keeps annotation effort focused where the model is weakest.