Detecting Covariate Shift in Satellite Input Features

Detect covariate shift in satellite input features with PSI, the KS test, and a classifier two-sample test to catch sensor drift before it degrades predictions.

TL;DR Answer

Covariate shift is when the distribution of your model’s inputs, P(X), changes while the relationship it learned, P(y|X), may stay the same. On satellite feeds it shows up as sensor calibration drift, a new season, or a changed atmospheric correction pipeline. Detect it per-band with the Population Stability Index (PSI) and the two-sample Kolmogorov–Smirnov test, and jointly with a classifier two-sample test — train a classifier to tell reference tiles from live tiles; a ROC AUC near 0.5 means the two are indistinguishable and there is no shift. This is one concrete diagnostic inside the wider Model Drift Detection for Geospatial Inference workflow.

import numpy as np

def psi(reference: np.ndarray, current: np.ndarray, bins: int = 10) -> float:
    """Population Stability Index between a reference and current sample."""
    edges = np.quantile(reference, np.linspace(0, 1, bins + 1))
    edges[0], edges[-1] = -np.inf, np.inf
    ref_pct = np.histogram(reference, edges)[0] / reference.size
    cur_pct = np.histogram(current, edges)[0] / current.size
    eps = 1e-6
    ref_pct = np.clip(ref_pct, eps, None)
    cur_pct = np.clip(cur_pct, eps, None)
    return float(np.sum((cur_pct - ref_pct) * np.log(cur_pct / ref_pct)))

Part of: Model Drift Detection for Geospatial Inference

What Covariate Shift Means for a Satellite Model

A deployed land-cover or biophysical model was fit on a reference distribution of inputs — the surface reflectances, spectral indices, and terrain features present in its training tiles. Covariate shift occurs when new imagery arrives with inputs drawn from a different distribution while the underlying physics the model encodes has not changed. The decision boundary is still correct; the data has simply moved into regions the model saw rarely or never during training, so its predictions become extrapolations.

Three causes dominate on real satellite feeds:

  • Sensor calibration drift. Optical sensors degrade over their mission lifetime. A slow radiometric gain change nudges every reflectance value a few percent, shifting the whole feature distribution without any change on the ground.
  • A new season. The phenological cycle drives large, legitimate swings in vegetation indices. A summer NDVI histogram looks nothing like a January one over the same fields.
  • A changed atmospheric correction. Reprocessing a scene with a new Sen2Cor or LaSRC version, or switching from top-of-atmosphere to surface reflectance, re-centres every band at once.

The danger is that covariate shift is silent. Prediction volumes and class ratios can look plausible while accuracy quietly collapses, because you have no ground-truth labels at inference time to measure it directly. Watching P(X) is the earliest signal you get, which is why it belongs at the front of any monitoring stack. This page is the concrete detection recipe; the broader retraining triggers, label-drift checks, and alerting policy live in the parent workflow.

Per-Feature Detection: PSI and the KS Test

Start with cheap, interpretable, one-band-at-a-time tests. They tell you which feature moved, which is exactly what an on-call engineer needs first.

PSI bins the reference sample by its quantiles, then measures how much probability mass moved between bins in the current sample. It is a symmetrised relative-entropy score. The widely used bands are: below 0.1 no meaningful shift, 0.1–0.25 moderate shift worth a look, and above 0.25 a major shift. Because it works on binned marginals it is robust to modest sample-size differences and trivially cheap to compute per tile.

The two-sample KS test compares the empirical cumulative distributions of the reference and current samples and returns the maximum vertical gap between them plus a p-value. It is more sensitive than PSI to shifts in the tails and needs no binning choice, but with the millions of pixels in a satellite scene even a trivially small, operationally meaningless difference becomes statistically significant. Treat the KS statistic (the distance) as the effect size and lean on PSI for the go/no-go decision; use the p-value only as a secondary sanity check.

A critical detail for raster inputs: ignore nodata. Cloud masks, scene edges, and fill values inject sentinel numbers (0, -9999, NaN) that will dominate any histogram if left in. Strip them before scoring, or every test fires on cloud cover rather than real drift. Consistent handling of nodata and band scaling is a shared concern with Feature Scaling for Geospatial Inputs, and the same masking logic should back both the training and monitoring paths.

import numpy as np
from scipy.stats import ks_2samp

def clean_band(arr: np.ndarray, nodata: float | None = None) -> np.ndarray:
    """Flatten a raster band and drop nodata / non-finite pixels."""
    flat = np.asarray(arr, dtype=float).ravel()
    flat = flat[np.isfinite(flat)]
    if nodata is not None:
        flat = flat[flat != nodata]
    return flat

def per_feature_report(
    reference: np.ndarray,
    current: np.ndarray,
    nodata: float | None = None,
) -> dict[str, float]:
    """PSI plus KS statistic and p-value for a single band."""
    ref = clean_band(reference, nodata)
    cur = clean_band(current, nodata)
    if ref.size < 100 or cur.size < 100:
        raise ValueError("Need >= 100 valid pixels per band to score drift.")
    ks_stat, ks_p = ks_2samp(ref, cur)
    return {
        "psi": psi(ref, cur),
        "ks_stat": float(ks_stat),
        "ks_pvalue": float(ks_p),
        "n_ref": ref.size,
        "n_cur": cur.size,
    }

Joint Detection: The Classifier Two-Sample Test

Per-band tests miss shifts that live only in the joint distribution. If the correlation between the red and near-infrared bands changes — say a new atmospheric correction alters their ratio — each band’s marginal histogram can look untouched while the pair has clearly moved. NDVI and similar ratio features depend on exactly that joint structure.

The classifier two-sample test catches these. Label every reference sample 0 and every current sample 1, train a classifier to tell them apart under stratified cross-validation, and read the ROC AUC. If reference and live tiles are drawn from the same distribution, no classifier can separate them and the AUC sits at 0.5. An AUC meaningfully above 0.5 is direct evidence of multivariate covariate shift, and the classifier’s feature importances rank which bands drive it.

import numpy as np
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import cross_val_score, StratifiedKFold


def classifier_two_sample_test(
    reference: np.ndarray,
    current: np.ndarray,
    n_splits: int = 5,
    random_state: int = 42,
) -> dict[str, float]:
    """Detect joint covariate shift by trying to classify reference vs current.

    Parameters
    ----------
    reference, current : np.ndarray of shape (n_samples, n_features)
        Per-pixel feature matrices, nodata already removed.
    n_splits : int
        Stratified CV folds used to estimate held-out AUC.
    random_state : int
        Seed for the classifier and the fold shuffling.

    Returns
    -------
    dict with mean AUC, its std across folds, and an ``is_shift`` flag.
    """
    if reference.ndim != 2 or current.ndim != 2:
        raise ValueError("reference and current must be 2-D (n_samples, n_features).")
    if reference.shape[1] != current.shape[1]:
        raise ValueError("reference and current must have the same number of features.")

    X = np.vstack([reference, current])
    y = np.concatenate([
        np.zeros(reference.shape[0], dtype=int),
        np.ones(current.shape[0], dtype=int),
    ])

    clf = HistGradientBoostingClassifier(
        max_depth=4, max_iter=200, random_state=random_state
    )
    cv = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=random_state)
    auc = cross_val_score(clf, X, y, cv=cv, scoring="roc_auc")

    mean_auc = float(auc.mean())
    return {
        "auc_mean": mean_auc,
        "auc_std": float(auc.std()),
        "is_shift": mean_auc > 0.55,   # 0.5 == indistinguishable
    }

Two guardrails matter. Subsample balanced pixel sets — a full scene has millions of pixels, and you want the two classes roughly equal in size so the AUC is not an artefact of imbalance. And draw samples spatially spread across the scene; contiguous pixels are strongly correlated, so a classifier can exploit local texture as a shortcut rather than distributional difference. The reasoning is the same one that motivates block sampling in Handling Spatial Autocorrelation: treat spatially adjacent pixels as near-duplicates, not independent draws.

How the Three Tests Fit Together

Covariate shift detection flow for satellite features Reference tiles and current tiles are cleaned of nodata, then scored two ways: per-band PSI and KS tests on marginals, and a classifier two-sample test on the joint distribution. Both feed a single shift decision. Reference tiles training-era P(X) Current tiles live inference feed Drop nodata mask + finite Per-feature marginals PSI > 0.25 = major KS distance per band Joint distribution classifier two-sample AUC ≈ 0.5 = no shift Shift decision alert + name the drifting bands, then investigate cause Calibration drift sensor gain change Seasonal cycle expected, not drift Atmos. correction reprocessing change

Run both branches on every batch. PSI and the KS test name the offending band and are cheap enough to run continuously; the classifier test is the backstop for joint shifts the marginals hide. When either fires, the cause triage in the bottom row decides your response — a seasonal swing needs no action, while calibration drift or a correction change usually means retraining or a rollback.

Walkthrough: Scoring Reference vs Current Tiles

Read matched bands from a reference composite and a fresh acquisition, then run both layers of detection.

import numpy as np
import rasterio

BANDS = ["blue", "green", "red", "nir"]
NODATA = 0.0

def load_pixels(path: str, sample: int = 20_000, seed: int = 0) -> np.ndarray:
    """Return an (n_sample, n_band) matrix of valid, spatially spread pixels."""
    with rasterio.open(path) as src:
        arr = src.read().astype(float)              # (bands, rows, cols)
    stacked = arr.reshape(arr.shape[0], -1).T       # (pixels, bands)
    valid = np.all(np.isfinite(stacked), axis=1) & np.all(stacked != NODATA, axis=1)
    stacked = stacked[valid]
    rng = np.random.default_rng(seed)
    idx = rng.choice(stacked.shape[0], size=min(sample, stacked.shape[0]), replace=False)
    return stacked[idx]

ref = load_pixels("reference_2024_summer.tif")
cur = load_pixels("current_batch.tif")

# Per-band marginal tests
for i, band in enumerate(BANDS):
    report = per_feature_report(ref[:, i], cur[:, i], nodata=None)
    flag = "SHIFT" if report["psi"] > 0.25 else "ok"
    print(f"{band:>5}: PSI={report['psi']:.3f}  KS={report['ks_stat']:.3f}  [{flag}]")

# Joint test across all bands
joint = classifier_two_sample_test(ref, cur)
print(f"joint AUC={joint['auc_mean']:.3f} ± {joint['auc_std']:.3f}  shift={joint['is_shift']}")

To keep the comparison honest, the reference must come from the same phenological season as the current batch — here a prior summer, not an annual mean. Pin the exact reference composite so the baseline is reproducible; versioning those reference tiles alongside the model is what Spatial Dataset Versioning with DVC and lakeFS exists to guarantee, so a drift alert always points at a fixed, retrievable baseline.

Verification: Inject a Shift, Confirm No False Alarm

A drift detector you have not tested will either cry wolf or stay silent when it matters. Prove both directions before trusting it in production.

import numpy as np

rng = np.random.default_rng(7)
base = rng.normal(0.25, 0.05, size=(20_000, 4))   # stand-in for 4 surface-reflectance bands

# 1. Identical distribution -> no alarm
copy = rng.normal(0.25, 0.05, size=(20_000, 4))
assert psi(base[:, 0], copy[:, 0]) < 0.1
assert not classifier_two_sample_test(base, copy)["is_shift"]

# 2. Inject a calibration-style gain shift on the NIR band -> detected
shifted = copy.copy()
shifted[:, 3] = shifted[:, 3] * 1.15 + 0.03        # +15% gain, small offset
assert psi(base[:, 3], shifted[:, 3]) > 0.25        # marginal test catches it
assert classifier_two_sample_test(base, shifted)["is_shift"]  # joint test agrees

print("drift detector verified: silent on identical data, fires on injected shift")

The first block confirms that two independent draws from the same distribution produce a PSI under 0.1 and a classifier AUC near 0.5 — no false alarm. The second injects a +15% radiometric gain, exactly the fingerprint of sensor calibration drift, and confirms both the per-band and joint tests fire. Wire these two assertions into CI so a change to your masking or binning logic can never silently break detection. From here, connect the alert to the retraining and rollback policy in the Model Drift Detection for Geospatial Inference workflow.

FAQ

What PSI threshold should trigger an alert for satellite bands?

The convention treats PSI below 0.1 as no meaningful shift, 0.1–0.25 as moderate shift worth investigating, and above 0.25 as a major shift that warrants retraining or rollback. For noisy spectral bands, calibrate against your own historical tile-to-tile variation instead of adopting the defaults blindly: compute PSI between two reference periods you know are stable, and set the alert level above that natural baseline so ordinary noise does not page you at 3 a.m.

How do I tell a normal seasonal cycle from true covariate shift?

Compare each incoming window against a reference drawn from the same phenological season in prior years, not against an annual average. A March acquisition should be scored against previous Marches. Vegetation indices swing widely between summer and winter, so an annual reference makes every winter tile look shifted. If a same-season comparison still flags a shift, the cause is more likely sensor calibration, an atmospheric correction change, or genuine land-cover change than phenology.

Why use a classifier two-sample test instead of per-feature tests?

Per-feature tests such as PSI and the KS test examine one band at a time and miss shifts that appear only in the joint distribution — for example when the red/near-infrared correlation changes while each band’s marginal stays fixed. A classifier two-sample test trains a model to distinguish reference from live samples: if it cannot separate them the AUC sits near 0.5, and an AUC well above 0.5 signals multivariate shift while ranking which features drive it.

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