Monitoring Prediction Distribution Shift in Production

Watch the output side of a geospatial model: per-tile class shares, PSI against a reference, seasonal baselines, spatial attribution of shift, and alert thresholds that do not cry wolf.

Log the per-tile class shares of every run, compare them against a seasonal reference with a population stability index, and attribute any shift back to the tiles that caused it. Output-side monitoring is cheap — the predictions already exist — and it catches an entire class of failures that input monitoring misses, including model rollbacks, preprocessing regressions and silent partial runs.

This page is the output-side method. For covariate shift on the input side, see Model Drift Detection for Geospatial Inference.

Predictions in, class shares out, compared against the same month last year Each production run emits per-tile class shares. Those are aggregated and compared against a seasonal reference built from the same calendar month in previous years, producing a population stability index. When the index exceeds the calibrated threshold, the contribution of each tile is computed so the shift can be located on a map. Monitor the output; it is already computed and it is cheap run output 6 000 tiles class shares per tile, per class seasonal reference same month, prior years PSI calibrated threshold attribute to tiles map, then investigate A shift concentrated in a few tiles is a data problem; a shift spread evenly is a model problem. Attribution is what turns an alert into an action.

Why This Fails in Geospatial ML Pipelines

Seasonality swamps everything. Land cover class shares in a temperate zone swing by tens of percent between March and July, so a monitor comparing this week against last week fires continuously in spring and autumn and is muted within a month. The comparison has to be against the same phase of the cycle — the same calendar month in previous years — or the signal is buried under phenology.

The second failure is aggregate-only monitoring. A national class share can be perfectly stable while one province flips from cropland to bare soil and another compensates. Monitoring only the total is monitoring the sum of two errors, and it will not fire until both go the same way. Per-tile shares, aggregated into a stability index and then attributed back, catch it.

Third, thresholds imported from finance. The conventional 0.1 / 0.25 population-stability-index bands were calibrated on credit scorecards, not on satellite mosaics whose share vectors move with cloud cover and sun angle. Adopting them unexamined produces either an alert every week or an alert never. The threshold has to come from your own history of runs known to be fine.

Core Principles

  • Compare to the same season, not to the previous run.
  • Store per-tile shares, not just national totals.
  • Calibrate the threshold on your own quiet history.
  • Attribute every alert to tiles before anyone is asked to act on it.
  • Log the model version with every share record. Half of all shifts are deployments.
  • Pair with input monitoring. Output shift alone cannot say why.

Production-Ready Code

from __future__ import annotations

import logging
from dataclasses import dataclass

import numpy as np
import pandas as pd
import rasterio

logger = logging.getLogger(__name__)
IGNORE = 255


def tile_class_shares(tile_path: str, n_classes: int) -> np.ndarray:
    """Fraction of valid pixels in each class, for one tile."""
    with rasterio.open(tile_path) as src:
        arr = src.read(1)
    valid = arr != IGNORE
    if not valid.any():
        return np.full(n_classes, np.nan)
    counts = np.bincount(arr[valid], minlength=n_classes)[:n_classes]
    return counts / counts.sum()


def collect_run_shares(tile_paths: list[str], n_classes: int, run_id: str,
                       model_version: str, scene_date: str) -> pd.DataFrame:
    """One row per tile per run — the record everything else is built on."""
    rows = []
    for path in tile_paths:
        shares = tile_class_shares(path, n_classes)
        rows.append({"run_id": run_id, "model_version": model_version,
                     "scene_date": scene_date,
                     "tile_id": path.rsplit("/", 1)[-1].removesuffix(".tif"),
                     **{f"class_{i}": float(s) for i, s in enumerate(shares)}})
    return pd.DataFrame(rows)


def psi(current: np.ndarray, reference: np.ndarray, eps: float = 1e-6) -> float:
    """Population stability index between two share vectors."""
    c = np.clip(np.asarray(current, dtype="float64"), eps, None)
    r = np.clip(np.asarray(reference, dtype="float64"), eps, None)
    c, r = c / c.sum(), r / r.sum()
    return float(np.sum((c - r) * np.log(c / r)))


def seasonal_reference(history: pd.DataFrame, month: int, n_classes: int,
                       min_runs: int = 3) -> np.ndarray:
    """Median class shares for this calendar month across previous years."""
    hist = history.copy()
    hist["month"] = pd.to_datetime(hist["scene_date"]).dt.month
    same = hist.loc[hist["month"] == month]
    if same["run_id"].nunique() < min_runs:
        raise ValueError(
            f"only {same['run_id'].nunique()} historical runs for month {month}; "
            "not enough to build a seasonal reference")
    cols = [f"class_{i}" for i in range(n_classes)]
    return same.groupby("run_id")[cols].mean().median().to_numpy()


def calibrate_threshold(history: pd.DataFrame, n_classes: int,
                        quantile: float = 0.99, margin: float = 1.5) -> float:
    """Threshold from quiet history: above anything normal operation has produced."""
    cols = [f"class_{i}" for i in range(n_classes)]
    per_run = history.groupby(["run_id", "scene_date"])[cols].mean().reset_index()
    per_run["month"] = pd.to_datetime(per_run["scene_date"]).dt.month

    values = []
    for month, grp in per_run.groupby("month"):
        if len(grp) < 3:
            continue
        ref = grp[cols].median().to_numpy()
        values += [psi(row, ref) for row in grp[cols].to_numpy()]
    threshold = float(np.quantile(values, quantile) * margin)
    logger.info("calibrated PSI threshold %.4f from %d historical runs",
                threshold, len(values))
    return threshold


def attribute_shift(current: pd.DataFrame, reference: np.ndarray,
                    n_classes: int, top: int = 20) -> pd.DataFrame:
    """Which tiles moved the aggregate? Per-tile PSI against the same reference."""
    cols = [f"class_{i}" for i in range(n_classes)]
    out = current.copy()
    out["tile_psi"] = [psi(row, reference) for row in out[cols].to_numpy()]
    return out.nlargest(top, "tile_psi")[["tile_id", "tile_psi", *cols]]

Step-by-Step Walkthrough

Step 1 — Record shares on every run

shares = collect_run_shares(written_tiles, n_classes=5, run_id=spec.run_id,
                            model_version=spec.model_version,
                            scene_date=spec.scene_date)
shares.to_parquet(f"s3://monitoring/shares/{spec.scene_date}/{spec.run_id}.parquet")

This runs in seconds — the tiles are already written — and it is the only data the rest of the monitor needs.

Step 2 — Build the seasonal reference and the threshold once

history = pd.read_parquet("s3://monitoring/shares/")
reference = seasonal_reference(history, month=8, n_classes=5)
threshold = calibrate_threshold(history, n_classes=5)
print("reference shares:", np.round(reference, 4))
print("alert above PSI:", round(threshold, 4))

Step 3 — Score the run and attribute

cols = [f"class_{i}" for i in range(5)]
current = shares[cols].mean().to_numpy()
score = psi(current, reference)
print(f"run PSI {score:.4f} (threshold {threshold:.4f})")

if score > threshold:
    worst = attribute_shift(shares, reference, n_classes=5, top=20)
    print(worst.to_string(index=False))

Read the attribution before writing the incident note. Twenty tiles carrying the whole shift points at a data problem — one province, one source scene, one corrupt input. A shift spread evenly across thousands of tiles points at the model or the preprocessing.

Step 4 — Check the deployment log first

recent = history.loc[history["scene_date"] >= "2026-07-25", ["scene_date", "model_version"]]
print(recent.drop_duplicates().sort_values("scene_date").to_string(index=False))

A version change on the same day as the shift explains most alerts, and it is a one-line check that saves hours of investigation.

Step 5 — Alert with the context attached

if score > threshold:
    notify(
        title=f"Prediction shift on {spec.scene_date}: PSI {score:.3f}",
        body=(f"threshold {threshold:.3f}, model {spec.model_version}\n"
              f"top tiles: {', '.join(worst['tile_id'].head(5))}\n"
              f"class shares now {np.round(current, 3).tolist()} "
              f"vs reference {np.round(reference, 3).tolist()}"),
    )
Week-on-week fires all spring; seasonal fires once, correctly Two stability index traces run across a year. The week over week comparison rises above the alert line repeatedly through spring and autumn as phenology changes. The seasonal comparison stays flat all year and rises above the line only once, at the point where a model version was deployed. Stability index over one year FebMayAugNov alert line v7 deployed week-on-week seasonal The green trace has one alert all year, and it is the one worth reading.

Verification

import numpy as np
import pandas as pd
import pytest


def test_psi_is_zero_for_identical_distributions():
    v = np.array([0.5, 0.3, 0.15, 0.05])
    assert psi(v, v) == pytest.approx(0.0, abs=1e-12)


def test_psi_grows_with_divergence():
    ref = np.array([0.5, 0.3, 0.2])
    near = np.array([0.48, 0.31, 0.21])
    far = np.array([0.2, 0.3, 0.5])
    assert psi(near, ref) < psi(far, ref)


def test_zero_share_does_not_produce_inf():
    assert np.isfinite(psi(np.array([1.0, 0.0]), np.array([0.5, 0.5])))


def test_seasonal_reference_requires_history():
    thin = pd.DataFrame({"run_id": ["a"], "scene_date": ["2026-08-04"],
                         "class_0": [0.5], "class_1": [0.5]})
    with pytest.raises(ValueError, match="not enough"):
        seasonal_reference(thin, month=8, n_classes=2)


def test_attribution_finds_the_injected_tile():
    """A single tile flipped to one class must top the attribution table."""
    df = base_shares.copy()
    df.loc[df.index[7], ["class_0", "class_1", "class_2"]] = [1.0, 0.0, 0.0]
    worst = attribute_shift(df, reference=np.array([0.4, 0.35, 0.25]), n_classes=3, top=3)
    assert df.loc[df.index[7], "tile_id"] in set(worst["tile_id"])

The attribution test is the one that keeps the monitor useful: an alert nobody can localise gets acknowledged and ignored, which is the same as having no monitor.

Reading the shape of the shift Three interpretations are shown. A shift concentrated in a handful of tiles points to a data problem in a specific area. A shift spread evenly across all tiles points to the model or the preprocessing. A shift that coincides with a model deployment is usually the deployment itself. Where is the shift, and when did it start? concentrated a data problem one province, one source scene diffuse model or preprocessing check the scaler and the version step change a deployment check the version log first All three are worth knowing; only the first two need a data investigation.

FAQ

Why compare against a seasonal baseline instead of last week?

Because predictions move with phenology. A week-on-week comparison flags every tile in spring. Comparing this August against previous Augusts removes the seasonal signal and leaves the change you want to see.

What PSI value should trigger an alert?

Calibrate it. The conventional 0.1 / 0.25 bands come from credit scoring; compute the index across a year of quiet operation and set the threshold above the maximum you have actually observed.

Should prediction drift alone stop a deployment?

No — it should trigger investigation. Pair it with input-side monitoring and a small labelled audit sample so the alert can be resolved, not just acknowledged.

Part of: Model Drift Detection for Geospatial Inference Part of: Geospatial MLOps and Model Deployment