Model Drift Detection for Geospatial Inference

Q: What PSI threshold indicates real drift in spectral bands?

The conventional bands are PSI below 0.1 for a stable distribution, 0.1 to 0.25 for a moderate shift worth investigating, and above 0.25 for significant drift that warrants retraining. For satellite bands, calibrate these against a seasonal baseline because vegetation phenology alone can push a healthy band to 0.1 to 0.2.

Q: Should I use PSI or the KS test for drift detection?

Use both. PSI is a binned, symmetric divergence that is robust to sample size and easy to threshold, while the Kolmogorov-Smirnov statistic is a distribution-free test sensitive to shifts in the cumulative distribution. PSI is better for routine monitoring dashboards; KS adds a p-value for statistical confirmation when PSI flags a band.

Q: How do I stop seasonal change from triggering false drift alerts?

Build a per-season reference distribution rather than a single global baseline, compare each incoming tile against the matching season, and use spectral indices that are more stable than raw bands. Persistent drift across multiple seasons is real; a spike that reverts on the next seasonal cycle is phenology.

Q: Does input drift always mean the model is degrading?

No. Covariate shift in the inputs raises the probability of degraded predictions but does not confirm it. Confirm concept drift by monitoring prediction distributions and, where ground truth arrives later, by tracking rolling accuracy. Treat input drift as an early warning that gates a more expensive label-based evaluation.

A land-cover classifier that scored 0.92 macro-F1 on last year’s validation set is not guaranteed to score anything close to that on the tiles it processed this morning. Sensors degrade and are recalibrated, atmospheric correction pipelines change, new satellites join a constellation with slightly different spectral response functions, and the landscape itself changes through seasons, fires, and construction. Any of these silently shifts the distribution of pixel values feeding your model, and because inference pipelines rarely receive immediate ground truth, the damage is invisible until someone audits a badly wrong map weeks later. Drift detection is the discipline of catching that shift from the inputs alone, before it becomes a wrong decision.

This is part of Geospatial MLOps and Model Deployment, which covers the full operational lifecycle from packaging models through monitoring them in production. This guide focuses on the monitoring half: how to establish a reference distribution from your training data, quantify covariate shift per spectral band with the population stability index (PSI) and the Kolmogorov-Smirnov (KS) test, schedule those checks over a stream of incoming raster tiles, and wire the results into an alert that can trigger retraining. The companion guide on detecting covariate shift in satellite input features drills into the spectral-feature specifics.

Problem Framing

Drift comes in two flavours, and confusing them wastes effort. Data drift (also called covariate shift) is a change in the distribution of the model’s inputs — the pixel values, spectral indices, or engineered features — while the relationship between inputs and target stays fixed. Concept drift is a change in that relationship itself: the same input now maps to a different correct output, for example when a new crop variety changes the spectral signature of “healthy vegetation.” In deployed geospatial systems, data drift is the one you can detect immediately and cheaply, because it needs no labels; you only need the incoming rasters and a memory of what the training data looked like.

The reason this matters so much for satellite and aerial inference is that the input space is high-dimensional, continuous, and physically coupled to the environment. A model trained on summer Sentinel-2 imagery over temperate farmland sees a near-different planet in winter: near-infrared reflectance collapses as canopies senesce, shortwave-infrared responds to soil moisture, and cloud-adjacency effects change the tail of every band. Some of that is genuine, recoverable seasonality; some of it is real drift from a recalibrated sensor or an expanded geographic footprint. The job of a drift monitor is to quantify the shift objectively so a human — or an automated policy — can decide whether the model is still operating inside the distribution it was trained for.

Two nuances make geospatial drift detection different from tabular monitoring. First, nodata and cloud-masked pixels are pervasive and, if not excluded, they inject a spike at a sentinel value (often 0 or -9999) that dominates every histogram and manufactures fake drift. Second, spatial autocorrelation means a single tile is not a random sample of the landscape, so per-tile statistics are noisy; you aggregate across tiles to get a stable read. Both concerns are handled explicitly in the code below. If your upstream normalization statistics are themselves computed per-batch rather than pinned, drift and stale scaling can masquerade as one another — see feature scaling for geospatial inputs for why the reference statistics must be frozen alongside the model.

Prerequisites & Environment Setup

Drift metrics are only comparable if the environment that computes them is pinned. A scipy minor-version change can alter the KS statistic’s tie-handling, and a numpy histogram edge-case change can shift PSI at the third decimal — enough to flip a borderline band across a threshold.

# Pinned requirements for geospatial drift detection
evidently==0.4.30
scipy==1.13.1
numpy==1.26.4
rasterio==1.3.10
pandas==2.2.2

Install with:

pip install "evidently==0.4.30" "scipy==1.13.1" "numpy==1.26.4" \
            "rasterio==1.3.10" "pandas==2.2.2"

GDAL system dependencies (Ubuntu/Debian):

sudo apt-get install -y gdal-bin libgdal-dev

The reference must be immutable and versioned

The single most important operational rule is that the reference distribution is an artefact, not a computation you re-run. It is derived once from the exact training data that produced the deployed model, serialised, and versioned alongside the model weights. If you recompute the reference from “recent” data you will define drift away — the monitor will always report that today looks like today. Store the reference JSON in the same versioned store as your training data; the spatial dataset versioning with DVC and lakeFS guide covers how to keep that artefact reproducible.

import rasterio

# A projected, consistent grid is assumed for all tiles compared against one reference.
with rasterio.open("training_scene.tif") as src:
    assert src.count >= 4, "Expected at least 4 spectral bands for the reference"
    assert src.crs is not None, "Reference raster has no CRS — cannot guarantee tile alignment"
    print(f"Reference: {src.count} bands, CRS {src.crs}, nodata {src.nodata}")

Step-by-Step Implementation

Step 1 — Compute a reference feature distribution from training data

The reference captures, per band, the bin edges and the normalised histogram of the training pixels. Using fixed quantile-based edges (rather than uniform edges) makes PSI stable for skewed spectral bands, where most reflectance values cluster in a narrow range with a long bright tail. Critically, nodata pixels are masked before any statistic is computed.

import numpy as np
import rasterio
import json


def masked_band_values(src: rasterio.io.DatasetReader, band: int) -> np.ndarray:
    """Read one band and return finite, non-nodata pixel values as a 1-D array."""
    arr = src.read(band).astype("float64")
    mask = np.ones(arr.shape, dtype=bool)
    if src.nodata is not None:
        mask &= arr != src.nodata
    mask &= np.isfinite(arr)
    values = arr[mask]
    if values.size == 0:
        raise ValueError(f"Band {band} has no valid pixels after masking nodata")
    return values


def build_reference(raster_path: str, n_bins: int = 10) -> dict:
    """Build a per-band reference: quantile bin edges + normalised histogram.

    Args:
        raster_path: Path to the training raster (multi-band).
        n_bins: Number of quantile bins per band.
    Returns:
        A JSON-serialisable dict keyed by band index.
    """
    reference = {"n_bins": n_bins, "bands": {}}
    with rasterio.open(raster_path) as src:
        for band in range(1, src.count + 1):
            values = masked_band_values(src, band)
            # Quantile edges; unique() guards against ties in low-variance bands.
            quantiles = np.linspace(0, 1, n_bins + 1)
            edges = np.unique(np.quantile(values, quantiles))
            if edges.size < 3:
                raise ValueError(f"Band {band} is near-constant; cannot bin")
            counts, _ = np.histogram(values, bins=edges)
            proportions = counts / counts.sum()
            reference["bands"][str(band)] = {
                "edges": edges.tolist(),
                "proportions": proportions.tolist(),
                "mean": float(values.mean()),
                "std": float(values.std()),
                "n": int(values.size),
            }
    return reference


reference = build_reference("training_scene.tif", n_bins=10)
with open("drift_reference.json", "w") as f:
    json.dump(reference, f)

# Inline validation: proportions per band must sum to 1.
for band, stats in reference["bands"].items():
    total = sum(stats["proportions"])
    assert abs(total - 1.0) < 1e-9, f"Band {band} proportions sum to {total}, not 1"
print(f"Reference built for {len(reference['bands'])} bands")

Step 2 — Compute PSI and the KS test per spectral band

The population stability index compares two binned distributions with a symmetric divergence. For each bin i it computes (live_i - ref_i) * ln(live_i / ref_i) and sums over bins. Empty bins would send the logarithm to infinity, so both proportions are floored with a small epsilon. The KS test is computed on the raw values and returns a distribution-free statistic plus a p-value.

import numpy as np
from scipy.stats import ks_2samp

_EPS = 1e-6


def population_stability_index(
    ref_proportions: np.ndarray,
    live_values: np.ndarray,
    edges: np.ndarray,
) -> float:
    """PSI of live_values against a reference histogram defined by edges.

    Args:
        ref_proportions: Reference bin proportions (sum to 1).
        live_values: 1-D array of live pixel values (nodata already removed).
        edges: Bin edges used to build the reference.
    Returns:
        Non-negative PSI. Higher means greater divergence.
    """
    live_counts, _ = np.histogram(live_values, bins=edges)
    live_prop = live_counts / max(live_counts.sum(), 1)
    ref = np.clip(ref_proportions, _EPS, None)
    live = np.clip(live_prop, _EPS, None)
    return float(np.sum((live - ref) * np.log(live / ref)))


def band_drift(reference_band: dict, live_values: np.ndarray,
               ref_sample: np.ndarray | None = None) -> dict:
    """Compute PSI and a KS statistic for one band."""
    edges = np.asarray(reference_band["edges"])
    ref_prop = np.asarray(reference_band["proportions"])
    psi = population_stability_index(ref_prop, live_values, edges)

    result = {"psi": psi}
    if ref_sample is not None and ref_sample.size and live_values.size:
        ks_stat, p_value = ks_2samp(ref_sample, live_values)
        result["ks_stat"] = float(ks_stat)
        result["ks_pvalue"] = float(p_value)
    return result

A quick assertion confirms the metric behaves: identical distributions must score zero PSI, and a shifted one must score higher.

rng = np.random.default_rng(0)
base = rng.normal(0.3, 0.05, 20000)
edges = np.unique(np.quantile(base, np.linspace(0, 1, 11)))
ref_prop = np.histogram(base, bins=edges)[0] / base.size

same = population_stability_index(ref_prop, base, edges)
shifted = population_stability_index(ref_prop, base + 0.08, edges)
assert same < 0.01, f"Identical data should give near-zero PSI, got {same}"
assert shifted > 0.25, f"Shifted data should exceed the drift threshold, got {shifted}"
print(f"PSI sanity check passed: same={same:.4f}, shifted={shifted:.3f}")

Step 3 — Schedule drift checks over incoming raster tiles

Production tiles arrive continuously, so the check runs as a batch job over a catalogue of new tiles. Each tile is read band by band with the same masking as the reference, values are pooled across all tiles in the batch to defeat per-tile spatial autocorrelation, and one drift report is emitted per batch with a timestamp. Reading a modest reference sample once lets the KS test run without holding the entire training scene in memory.

import numpy as np
import rasterio
import pandas as pd
from datetime import datetime, timezone
from pathlib import Path


def load_reference_sample(raster_path: str, per_band: int = 20000) -> dict:
    """Draw a fixed random sample of valid pixels per band for KS testing."""
    rng = np.random.default_rng(42)
    samples = {}
    with rasterio.open(raster_path) as src:
        for band in range(1, src.count + 1):
            values = masked_band_values(src, band)
            take = min(per_band, values.size)
            samples[str(band)] = rng.choice(values, size=take, replace=False)
    return samples


def run_drift_batch(tile_paths: list[str], reference: dict,
                    ref_samples: dict) -> pd.DataFrame:
    """Pool pixels across a batch of tiles and score drift per band."""
    if not tile_paths:
        raise ValueError("No tiles supplied for the drift batch")

    pooled: dict[str, list[np.ndarray]] = {b: [] for b in reference["bands"]}
    for path in tile_paths:
        with rasterio.open(path) as src:
            available = {str(b) for b in range(1, src.count + 1)}
            missing = set(reference["bands"]) - available
            if missing:
                raise ValueError(f"{path} missing bands {sorted(missing)}")
            for band in reference["bands"]:
                pooled[band].append(masked_band_values(src, int(band)))

    ts = datetime.now(timezone.utc).isoformat()
    rows = []
    for band, ref_band in reference["bands"].items():
        live_values = np.concatenate(pooled[band])
        metrics = band_drift(ref_band, live_values, ref_samples.get(band))
        rows.append({"timestamp": ts, "band": int(band), **metrics,
                     "n_live": int(live_values.size)})
    return pd.DataFrame(rows)


tiles = sorted(str(p) for p in Path("incoming/").glob("*.tif"))
ref_samples = load_reference_sample("training_scene.tif")
report = run_drift_batch(tiles, reference, ref_samples)
report.to_parquet("drift_reports/2026-07-10.parquet", index=False)
print(report[["band", "psi", "ks_stat", "ks_pvalue"]].round(4))

Step 4 — Set an alert threshold and trigger retraining

A raw PSI number is not an action. The final step maps each band to a status using the conventional PSI bands — stable below 0.1, warning from 0.1 to 0.25, drifted above 0.25 — cross-checked against the KS p-value so a large PSI on a tiny sample does not fire spuriously. The policy escalates only when a band the model actually depends on drifts, and emits a structured payload that a scheduler or message queue can act on.

import pandas as pd

PSI_WARN, PSI_ALERT, KS_ALPHA = 0.10, 0.25, 0.01


def classify(row: pd.Series) -> str:
    if row["psi"] >= PSI_ALERT and row.get("ks_pvalue", 0.0) < KS_ALPHA:
        return "drifted"
    if row["psi"] >= PSI_WARN:
        return "warning"
    return "stable"


def evaluate_alert(report: pd.DataFrame, critical_bands: set[int]) -> dict:
    """Return an alert decision and per-band statuses for a drift report."""
    report = report.copy()
    report["status"] = report.apply(classify, axis=1)
    drifted = report[report["status"] == "drifted"]
    critical_hit = drifted[drifted["band"].isin(critical_bands)]

    decision = {
        "timestamp": report["timestamp"].iloc[0],
        "alert": not critical_hit.empty,
        "trigger_retraining": not critical_hit.empty,
        "drifted_bands": drifted["band"].tolist(),
        "max_psi": float(report["psi"].max()),
        "statuses": dict(zip(report["band"], report["status"])),
    }
    return decision


# NIR (4) and SWIR (5, 6) usually dominate a land-cover model's decisions.
decision = evaluate_alert(report, critical_bands={4, 5, 6})
if decision["alert"]:
    print(f"DRIFT ALERT: bands {decision['drifted_bands']} "
          f"(max PSI {decision['max_psi']:.3f}) — retraining triggered")
else:
    print(f"No drift; max PSI {decision['max_psi']:.3f}")

The trigger_retraining flag is the integration point: publish decision to the queue that drives your retraining pipeline, which lives in the same containerized environment as inference so the rebuilt model ships through the identical path — see containerizing geospatial inference pipelines.

Verification & Testing

Trust in a drift monitor comes from proving it fires when it should and stays quiet when it should not. Two synthetic tests give that assurance: a null test where live data is drawn from the same distribution as the reference (expect all bands stable), and an injection test where a known shift is applied to one band (expect exactly that band to alert).

import numpy as np
import pandas as pd


def synthetic_report(reference: dict, shift: dict[int, float],
                     n: int = 30000, seed: int = 7) -> pd.DataFrame:
    """Generate a drift report from reference stats with optional per-band shifts."""
    rng = np.random.default_rng(seed)
    rows = []
    ts = "2026-07-10T00:00:00+00:00"
    for band, stats in reference["bands"].items():
        live = rng.normal(stats["mean"] + shift.get(int(band), 0.0),
                          stats["std"], n)
        ref_sample = rng.normal(stats["mean"], stats["std"], n)
        metrics = band_drift(stats, live, ref_sample)
        rows.append({"timestamp": ts, "band": int(band), **metrics, "n_live": n})
    return pd.DataFrame(rows)


# Null test — no shift, nothing should alert.
null_report = synthetic_report(reference, shift={})
null_decision = evaluate_alert(null_report, critical_bands={4})
assert not null_decision["alert"], "False positive: null data raised a drift alert"

# Injection test — shift band 4 by 2 standard deviations.
sd4 = reference["bands"]["4"]["std"]
inj_report = synthetic_report(reference, shift={4: 2 * sd4})
inj_decision = evaluate_alert(inj_report, critical_bands={4})
assert inj_decision["alert"], "False negative: injected 2-sigma shift not detected"
assert 4 in inj_decision["drifted_bands"], "Wrong band flagged"
print("Drift monitor verified: null quiet, injection detected")

Evidently provides a higher-level report that cross-checks your hand-rolled PSI and adds a browsable HTML artefact for the monitoring dashboard. It is worth running in parallel as an independent confirmation.

import pandas as pd
from evidently.report import Report
from evidently.metric_preset import DataDriftPreset

# Column-per-band frames of sampled pixels (reference vs current batch).
ref_df = pd.DataFrame({f"band_{b}": s for b, s in ref_samples.items()})
cur_df = pd.DataFrame({f"band_{b}": np.concatenate(
    [masked_band_values(rasterio.open(t), int(b)) for t in tiles])[: len(ref_df)]
    for b in reference["bands"]})

drift_report = Report(metrics=[DataDriftPreset(stattest="psi", stattest_threshold=0.25)])
drift_report.run(reference_data=ref_df, current_data=cur_df)
drift_report.save_html("drift_reports/evidently_2026-07-10.html")
summary = drift_report.as_dict()
print("Dataset drift detected:",
      summary["metrics"][0]["result"]["dataset_drift"])

Troubleshooting & Common Errors

NaN or Inf PSI because a bin is empty

Cause: A live batch populates zero pixels in a reference bin, so live_i / ref_i is 0 and ln(0) is -inf, or a reference bin is empty and the ratio explodes.

Fix: Floor both proportion vectors with an epsilon before the logarithm, exactly as population_stability_index does with np.clip(..., _EPS, None). Never compute PSI on raw counts.

Nodata pixels manufacturing fake drift

Cause: Cloud-masked or fill pixels carry a sentinel value (0, -9999, or NaN). Left in the histogram they create a giant spike that swamps the real distribution and reports catastrophic drift on every band.

Fix: Always route reads through masked_band_values, which drops src.nodata and non-finite pixels. If clouds are flagged in a separate mask band, apply that mask too before pooling.

Comparing distributions across mismatched CRS or resolution

Cause: Live tiles delivered in a different projection or pixel size than the training scene change the pixel population — resampling smooths tails and reprojection re-weights land-cover proportions — which reads as drift that is really a preprocessing mismatch.

Fix: Enforce that every tile shares the reference CRS and resolution before scoring; reproject in the ingest step if not. Assert src.crs == reference_crs and src.res == reference_res at read time, and standardise projection handling per CRS alignment and projection handling.

Seasonal signal mistaken for drift

Cause: A single global reference compared against winter tiles flags every vegetation-sensitive band because phenology, not sensor or model failure, moved the distribution.

Fix: Maintain a reference per season (or per phenological phase) and compare each tile against the matching baseline. Prefer spectral indices that normalise illumination and canopy state over raw bands, and treat drift as real only when it persists across the seasonal cycle.

KS test always significant on large tiles

Cause: With millions of pixels the KS test has enormous power and returns p < 0.001 for trivially small, operationally meaningless differences.

Fix: Do not gate on the p-value alone. Use PSI (an effect size) as the primary trigger and the KS p-value only as a secondary confirmation, or subsample each band to a fixed size (for example 20,000 pixels) before the KS test so its sensitivity is bounded.

Per-tile scores swing wildly between batches

Cause: A single tile is a spatially autocorrelated patch of one land-cover type, so its distribution is not representative and PSI jumps around.

Fix: Pool pixels across all tiles in a batch before scoring, as run_drift_batch does, and require a minimum pooled pixel count before you trust a report. This is the same spatial dependence covered in handling spatial autocorrelation.

Performance Optimisation

Read windowed blocks instead of whole tiles

Full-scene src.read() loads every band into memory at once. For large mosaics, iterate over internal tiling with rasterio’s block windows and accumulate histogram counts incrementally, so memory stays flat regardless of scene size.

import numpy as np
import rasterio


def streaming_counts(path: str, band: int, edges: np.ndarray,
                     nodata) -> np.ndarray:
    """Accumulate histogram counts over block windows without loading the full band."""
    counts = np.zeros(len(edges) - 1, dtype=np.int64)
    with rasterio.open(path) as src:
        for _, window in src.block_windows(band):
            arr = src.read(band, window=window).astype("float64")
            mask = np.isfinite(arr)
            if nodata is not None:
                mask &= arr != nodata
            block_counts, _ = np.histogram(arr[mask], bins=edges)
            counts += block_counts
    return counts

Accumulate histograms, not raw values

PSI needs only bin counts, so there is no reason to hold every pixel. Accumulating counts per band across a batch turns an O(pixels) memory footprint into O(bins), letting a single worker monitor continental mosaics. Keep the raw-value path only for the subsampled KS test.

Parallelise across tiles

Tiles are independent until the pooling step, so histogram accumulation maps cleanly across a process pool. Distribute tiles to workers, return per-tile count vectors, and sum them in the parent — this scales the batch job linearly with cores and is the cheapest speed-up for a large catalogue.

Cache the reference sample

The KS reference sample is deterministic for a fixed seed and training scene. Persist it next to drift_reference.json so every batch reuses it instead of re-reading the training raster, removing the largest per-run I/O cost.

FAQ

What PSI threshold indicates real drift in spectral bands?

The conventional bands are PSI below 0.1 for a stable distribution, 0.1 to 0.25 for a moderate shift worth investigating, and above 0.25 for significant drift that warrants retraining. For satellite bands, calibrate these against a seasonal baseline because vegetation phenology alone can push a healthy band to 0.1 to 0.2. Track the PSI trend over successive batches rather than reacting to a single value — a steadily climbing PSI is a more reliable retraining signal than one noisy spike.

Should I use PSI or the KS test for drift detection?

Use both, for different jobs. PSI is a binned, symmetric divergence that is robust to sample size and trivial to threshold, which makes it ideal for routine dashboards and automated policies. The Kolmogorov-Smirnov statistic is distribution-free and sensitive to shifts anywhere in the cumulative distribution, and its p-value gives a statistical confirmation when PSI flags a band. Gate on PSI, confirm with KS, and never gate on a KS p-value alone at raster scale where it is almost always significant.

How do I stop seasonal change from triggering false drift alerts?

Build a per-season reference distribution rather than a single global baseline, compare each incoming tile against the matching season, and favour spectral indices that are more stable than raw bands. Persistent drift across multiple seasons is real degradation; a spike that reverts on the next seasonal cycle is phenology. If normalization statistics are recomputed per batch instead of frozen with the model, stale scaling can also look like seasonal drift — pin them as described in feature scaling for geospatial inputs.

Does input drift always mean the model is degrading?

No. Covariate shift raises the probability of degraded predictions but does not confirm it — a model can be robust to shifts in bands it barely uses. Confirm concept drift by monitoring the prediction distribution and, where ground truth arrives later, by tracking rolling accuracy against those labels. Treat input drift as a cheap, immediate early warning that gates a more expensive label-based evaluation, not as proof of failure on its own.

Detecting Covariate Shift in Satellite Input Features — the spectral-feature deep dive on quantifying input drift
ONNX Export for Geospatial Model Inference — package the monitored model into a portable inference runtime
Spatial Dataset Versioning with DVC and lakeFS — version the reference distribution alongside training data
Containerizing Geospatial Inference Pipelines — ship inference and retraining through one reproducible image
Feature Scaling for Geospatial Inputs — freeze normalization statistics so stale scaling is not mistaken for drift

Part of: Geospatial MLOps and Model Deployment