A geospatial model that scores 0.94 on a held-out fold is not a deployed system. Between that validation notebook and a service that returns reliable predictions over live satellite tiles lies an entire operational discipline: exporting the model into a portable runtime, packaging the geospatial dependency stack so it behaves identically on every host, versioning the terabytes of raster and vector data the model was trained on, and monitoring the input distribution so you know — before your users do — when the ground has shifted beneath the model. This is the operational half of the workflow, and it is where most geospatial ML projects quietly fail after a promising start.
The failure surface is distinctly geospatial. A conventional tabular model can drift when a marketing campaign changes user behaviour; a land-cover model drifts when a new Sentinel-2 processing baseline subtly shifts reflectance, when a wildfire re-paints a region overnight, or when your inference pipeline reprojects tiles through a different PROJ version than training used. The consequences hide well: predictions stay in a plausible range, no exception is raised, and the error only surfaces when someone overlays the output on ground truth. This page covers the four operational areas that contain that risk: model drift detection for geospatial inference, ONNX export for geospatial model inference, spatial dataset versioning with DVC and LakeFS, and containerizing geospatial inference pipelines.
Everything here assumes a model already trained well. If your predictions are wrong because of leakage or a poor evaluation protocol, no amount of deployment engineering will save them — start with training geospatial predictive models in Python and the spatial feature engineering that feeds it. The premise of this page is that the model is sound and the remaining work is making it reproducible, portable, observable, and durable in production.
Foundational Prerequisites
Geospatial MLOps compounds two dependency stacks that are each fragile on their own: the geospatial native libraries (GDAL, PROJ, GEOS) and the ML runtime stack (ONNX, model frameworks, monitoring). A version skew in either one produces predictions that are wrong without being obviously broken. Pin both, in one lockfile, and build every environment — training, CI, and the inference container — from the same pins.
# requirements.txt — tested stack for geospatial MLOps
onnxruntime==1.20.1
onnx==1.17.0
skl2onnx==1.18.0
dvc[s3]==3.59.1
mlflow==2.20.2
evidently==0.6.7
rasterio==1.4.3
geopandas==1.0.1
pyproj==3.7.1
numpy==2.2.4
xgboost==2.1.4The geospatial native libraries must be pinned at the system level too, because a container that ships a different GDAL than training used will reproject and resample tiles differently. On Debian slim base images:
apt-get install -y --no-install-recommends \
gdal-bin=3.8.* \
libgdal-dev=3.8.* \
libproj-dev=9.4.* \
proj-dataVerify the runtime matches the training environment before serving a single request. This check belongs in the container entrypoint, not in a notebook you ran once:
import onnxruntime as ort
import rasterio, pyproj
print(ort.__version__) # 1.20.1
print(rasterio.gdal_version()) # 3.8.x — MUST match training host
print(pyproj.proj_version_str) # 9.4.x
sess = ort.InferenceSession("model.onnx", providers=["CPUExecutionProvider"])
print([i.name for i in sess.get_inputs()]) # confirm input signatureA mismatch between the GDAL that trained the model and the GDAL that serves it is the single most common source of silent geospatial deployment error. Catch it at startup with an explicit assertion rather than discovering it as a slow drift in prediction quality weeks later.
Core Techniques
Model Drift Detection for Geospatial Inference
A deployed geospatial model sees a moving target. Sensors are recalibrated, atmospheric correction baselines are reprocessed, seasons turn, and land cover physically changes. Model drift detection for geospatial inference is the practice of continuously comparing the live input distribution against the distribution the model was trained on, and raising an alert before degraded predictions reach users. The baseline is not arbitrary: it should be the feature distribution from your spatial cross-validation folds, so the reference already reflects geographic variability rather than a single region.
The workhorse tests are distribution-distance metrics applied per feature. The Kolmogorov-Smirnov statistic and the Wasserstein (earth-mover) distance both quantify how far the current window of inputs has moved from the reference, without assuming normality:
import numpy as np
from scipy.stats import ks_2samp, wasserstein_distance
def feature_drift(reference: np.ndarray, current: np.ndarray) -> dict:
"""Per-feature drift score between a reference and a live window."""
ref = reference[np.isfinite(reference)]
cur = current[np.isfinite(current)]
ks_stat, p_value = ks_2samp(ref, cur)
return {
"ks": float(ks_stat),
"p_value": float(p_value),
"wasserstein": float(wasserstein_distance(ref, cur)),
"drifted": bool(p_value < 0.01),
}Frame drift spatially, not just globally: a model can be stable in aggregate while degrading badly over one biome or one acquisition date. Aggregate drift scores per tile and map them, so a hotspot of covariate shift is visible geographically. See detecting covariate shift in satellite input features for the full monitoring workflow with evidently.
ONNX Export for Geospatial Model Inference
Serving a model as a pickled scikit-learn or a .pt PyTorch file couples your inference host to the exact training framework version, and often to a Python interpreter you would rather not ship. ONNX export for geospatial model inference decouples the two: you export the trained graph once, into a portable format that onnxruntime executes with no framework dependency, identical numerics across hosts, and materially faster CPU inference — which matters when you are scoring millions of pixels per scene.
from skl2onnx import to_onnx
from skl2onnx.common.data_types import FloatTensorType
import numpy as np
# n_features must match the training feature matrix exactly
initial_type = [("input", FloatTensorType([None, n_features]))]
onnx_model = to_onnx(
trained_pipeline,
initial_types=initial_type,
target_opset=18,
)
with open("landcover.onnx", "wb") as f:
f.write(onnx_model.SerializeToString())The export contract is the input signature. A raster inference loop feeds pixels as a 2-D array of shape (n_pixels, n_features), and the feature order must match training exactly — ONNX will not catch a reordered column, it will simply produce wrong classes. Bake the feature order into a manifest that both training and serving read. See serving an ONNX land-cover model with onnxruntime for a complete tile-to-prediction service.
Spatial Dataset Versioning with DVC and LakeFS
Model reproducibility is impossible without data reproducibility, and geospatial datasets are uniquely hard to version: a single training set can be terabytes of Cloud-Optimized GeoTIFFs, satellite archives are silently reprocessed months after acquisition, and administrative boundaries are reissued yearly. Spatial dataset versioning with DVC and LakeFS tracks these large, mutable inputs with Git-like semantics, so a model trained today can be retrained from the identical bytes next year.
DVC stores a small hash pointer in Git while the raster bytes live in object storage, giving you a commit that pins exactly which imagery produced a model:
# Track a directory of cloud-optimized GeoTIFFs
dvc add data/sentinel2_2026_q2/
git add data/sentinel2_2026_q2.dvc .gitignore
git commit -m "Training imagery: Sentinel-2 L2A, 2026 Q2 composite"
# Push the bytes to remote object storage
dvc remote add -d s3remote s3://geo-ml-artifacts/dvc
dvc pushPair data versioning with a model registry — MLflow logging the ONNX artifact, its metrics, and the DVC data hash it was trained on — so every deployed model links back to an exact, retrievable input snapshot. This is what makes a drift-triggered retrain trustworthy rather than a leap of faith. See versioning cloud-optimized GeoTIFFs with DVC for large-raster patterns that avoid re-hashing terabytes on every commit.
Containerizing Geospatial Inference Pipelines
A geospatial inference service is only reproducible if its entire native stack — GDAL, PROJ, GEOS, and their data files — is frozen into an image. Containerizing geospatial inference pipelines covers building images that are both correct and small, since a naive GDAL install can balloon a container past 2 GB and slow every cold start on autoscaled inference nodes.
FROM python:3.12-slim AS base
RUN apt-get update && apt-get install -y --no-install-recommends \
gdal-bin libgdal-dev libproj-dev proj-data \
&& rm -rf /var/lib/apt/lists/*
ENV PROJ_LIB=/usr/share/proj
COPY requirements.txt .
RUN pip install --no-cache-dir -r requirements.txt
COPY model.onnx feature_manifest.json infer.py ./
ENTRYPOINT ["python", "infer.py"]Two failure modes dominate here: a missing or mismatched PROJ_LIB path silently disables datum transformations, and a multi-stage build that omits proj-data produces containers that reproject incorrectly without erroring. Multi-stage builds and slim base images shrink the image, but the PROJ data directory must survive the final stage. See building a slim GDAL Docker image for inference for a multi-stage build that reaches production-grade size without dropping datum grids.
MLOps and Pipeline Integration
The four techniques above are stages of one continuous loop, not independent tools. A production geospatial MLOps system wires them into a pipeline where a drift signal on the monitoring side can trigger a versioned retrain, a new ONNX artifact, and a rolling container deploy — all with full provenance from prediction back to source imagery.
Continuous Training Triggers
Retraining should be event-driven, not calendar-driven. When the drift monitor’s per-feature scores breach a threshold over a sustained window, it emits a retrain event. That event references the current DVC data hash, so the retrain job pulls a defined snapshot rather than “whatever is in the bucket today.” The retrained model is re-evaluated with the same spatial cross-validation strategies used originally — a fresh model that has not cleared spatial validation must never auto-promote, because random splits would report inflated gains from spatial autocorrelation.
Artifact Registry and Lineage
Every promoted model carries a lineage record: the ONNX file hash, the training data version, the fitted preprocessing artifacts (scalers, encoders, spatial weights), the library pins, and the validation metrics. MLflow’s model registry is the practical home for this. The rule is that no artifact travels without its transformer state — a model deployed with a scaler fitted on different data will produce systematically shifted features, exactly the class of silent error that feature scaling for geospatial inputs warns about at training time.
Reproducible Inference Contracts
The contract between training and serving is a versioned manifest: the feature order, the target CRS, the tile size and overlap buffer, the nodata value, and the expected input ranges per band. Both the training pipeline and the inference container read this manifest. When they diverge — a new band is added upstream, a CRS is changed — the container fails fast at startup rather than emitting plausible-but-wrong rasters for days.
Performance and Scalability
Tiled, Parallel Inference
Continental-scale prediction is embarrassingly parallel over tiles, but only if tiling is done correctly. Split incoming extents into fixed tiles with an overlap buffer sized to the largest neighbourhood any feature needs, run each tile through onnxruntime, then merge on the interior and discard the buffer.
import numpy as np
import onnxruntime as ort
import rasterio
from rasterio.windows import Window
sess = ort.InferenceSession("landcover.onnx", providers=["CPUExecutionProvider"])
input_name = sess.get_inputs()[0].name
def predict_window(src, window: Window) -> np.ndarray:
stack = src.read(window=window).astype(np.float32) # (bands, h, w)
bands, h, w = stack.shape
X = stack.reshape(bands, h * w).T # (pixels, bands)
logits = sess.run(None, {input_name: X})[0]
return logits.argmax(axis=1).reshape(h, w).astype(np.uint8)onnxruntime on CPU with batched pixel arrays typically outpaces a native scikit-learn predict loop by a wide margin, and it parallelises cleanly across tiles because each session call is independent and thread-safe for inference.
Cloud-Native I/O and Cold Starts
Read tiles directly from Cloud-Optimized GeoTIFFs over range requests rather than downloading whole scenes, and write predictions back as COGs so downstream consumers can do the same. On autoscaled inference, container cold-start time is dominated by image size, which is why the slim-image work in the containerization section pays off directly: a 400 MB image starts in a fraction of the time of a 2 GB one, and cold starts are on the critical path for bursty request patterns.
Memory Budgets
Budget roughly 3× a tile’s raw byte size as working memory per worker for the reshape, cast, and prediction buffers. A 512 × 512 × 12-band float32 tile is about 12 MB raw, so allow ~36 MB per inference thread. Size your tile and worker count so that peak concurrent tiles stay within the node’s memory, and prefer smaller tiles with more workers over large tiles that risk OOM-killing the container mid-scene.
Common Failure Modes
1. CRS drift between training and inference. The model was trained on features computed in EPSG:32633 (metric UTM), but the inference container reprojects incoming tiles through a subtly different pipeline — a different PROJ version, or a missing datum grid — so distances, slopes, and areas are computed on a slightly different grid. Predictions stay in range but degrade near tile edges and in high-relief terrain. Fix: pin the target CRS in the inference manifest, pin PROJ at the system level, and assert pyproj.proj_version_str matches the training host at container startup.
2. Tile-boundary artefacts. Neighbourhood features (moving-window texture, spatial lag) computed on an edge pixel see missing neighbours outside the tile, so predictions on tile borders differ systematically from interiors and produce a visible grid of seams in the merged output. Fix: process overlapping tiles with a buffer equal to the largest neighbourhood radius, predict on the buffered tile, and merge only the interior region.
3. Stale normalization statistics. The inference service loads a scaler or per-band normalization fitted on last year’s imagery while the incoming distribution has shifted, or worse, refits normalization per tile so every tile is scaled to its own statistics. Both produce features that no longer match training-time distributions. Fix: serialise the fitted transformer with the model, version it in the registry, and never refit at inference time.
4. Silent covariate shift. A new sensor processing baseline nudges reflectance, or a wildfire re-paints a region; the input distribution moves away from the training distribution while the model keeps emitting confident, in-range predictions. Nothing errors. Fix: run continuous distribution-distance monitoring (KS / Wasserstein) against the training baseline, aggregated per tile so the shift is visible geographically, and gate a retrain on sustained breaches.
5. GDAL version skew in containers. The image ships a different GDAL or PROJ than training used, so resampling, warping, and nodata handling differ. The model runs, tiles decode, predictions merge — and they are subtly wrong. Fix: pin gdal-bin and libgdal-dev to exact versions in the Dockerfile, assert rasterio.gdal_version() at startup, and build training, CI, and serving from the same base image.
6. ONNX feature-order drift. The exported graph expects features in training order, but an upstream change reorders or inserts a band and the serving code feeds the new order. ONNX runs happily and returns wrong classes with no exception. Fix: pin feature order in a manifest read by both sides, and validate the input signature (sess.get_inputs()) against it before serving.
Best Practices Checklist
FAQ
Why export a geospatial model to ONNX instead of serving the native pickle?
ONNX decouples the serving host from the training framework. A pickled scikit-learn or PyTorch model demands the exact framework version (and often the exact Python) at inference time, which couples your container to your training environment and makes upgrades risky. An ONNX graph runs under onnxruntime with no framework dependency, produces identical numerics across hosts, and delivers materially faster CPU inference — which matters when a single Sentinel-2 scene is tens of millions of pixels. The trade-off is the export contract: the input signature and feature order are frozen at export, so you must pin them in a manifest both sides read.
How is drift detection different for geospatial models than for tabular models?
The mechanisms of drift are geospatial and the monitoring must be too. Inputs drift because sensors are recalibrated, atmospheric-correction baselines are reprocessed, seasons turn, and land cover physically changes — not because of a UI change or a marketing campaign. Critically, a geospatial model can be stable in global aggregate while degrading badly over one biome or one acquisition date, so drift scores must be aggregated and mapped per tile, not computed only as a single global number. The reference baseline should come from spatial cross-validation folds so it already reflects geographic variability rather than one region.
What is the most common silent failure when deploying a geospatial model?
Version skew in the native geospatial stack — GDAL and PROJ — between training and serving. The container reprojects, resamples, or handles nodata slightly differently than the training environment did, so features are computed on a subtly different grid. The model runs without error, tiles decode, predictions merge, and everything looks plausible; the error only surfaces when someone overlays output on ground truth. Pin GDAL and PROJ to exact versions and assert them at container startup.
When should a deployed geospatial model be retrained?
Retrain on evidence, not on a schedule. When continuous drift monitoring shows per-feature distribution distance breaching a threshold over a sustained window — not a single noisy tile — emit a retrain event that references the current versioned data snapshot. Re-evaluate the retrained model with the same spatial cross-validation used originally, and only promote it if it clears that bar. A calendar-based retrain wastes compute when nothing has shifted and lags reality when something shifts fast, such as after a fire or flood.
How do I keep a GDAL inference container small enough for fast autoscaling?
Use a slim base image and a multi-stage build, install only the runtime GDAL/PROJ packages with --no-install-recommends, and clear the apt cache in the same layer. The one thing you must not strip is the PROJ data directory — dropping proj-data silently disables datum transformations. A carefully built image lands in the few-hundred-megabyte range instead of over 2 GB, and since cold-start time on autoscaled inference is dominated by image size, that difference is directly on the critical path for bursty workloads.