A tabular explanation of a spatial model answers the wrong question. “Elevation contributed 0.31” tells you nothing about where elevation mattered, and in geospatial work the where is usually the finding. A flood model may be driven by upslope area in the headwaters and by impervious surface in the city, with the two effects cancelling to an unremarkable global average. Mapping the attributions surfaces that structure immediately; a bar chart of mean absolute SHAP values hides it completely.
This topic is part of Training Geospatial Predictive Models in Python. It builds on the models trained in gradient boosting for raster data and depends absolutely on the evaluation discipline of spatial cross-validation strategies — explaining a model that was validated with random folds is explaining a memorisation artefact.
Problem Framing
Three properties of spatial data break the standard explainability recipe.
Autocorrelation makes permutation dishonest. Permutation importance measures the score drop when a column is shuffled. Shuffling elevation across a whole country produces rows where a 3,000 m elevation sits beside a coastal rainfall value — a combination the model never saw and the world cannot produce. The score collapses, and the feature is declared essential. Permuting within spatial blocks keeps the values geographically plausible and typically halves the apparent importance of smooth surfaces.
Coordinates are a trap. Give a boosted tree raw x and y and it will carve the target surface into rectangles. Under a random split that looks like the strongest feature in the model; under the blocked evaluation described in choosing block size for spatial block cross-validation it collapses to noise. Any explanation that ranks coordinates first is telling you about the validation design, not about the process.
Geospatial feature sets are heavily collinear by construction. Slope, curvature and topographic position all derive from the same elevation grid; NDVI and EVI share two of three bands. Collinear features split their credit, so each looks unimportant while the group is decisive. Grouped importance over a correlation clustering is the fix, and it changes conclusions far more often than choosing between SHAP variants does.
Prerequisites & Environment Setup
# Pinned requirements for spatial explainability
shap==0.45.1
xgboost==2.0.3
scikit-learn==1.5.0
geopandas==0.14.4
numpy==1.26.4
pandas==2.2.2
scipy==1.13.1
matplotlib==3.9.0
Install with:
pip install "shap==0.45.1" "xgboost==2.0.3" "scikit-learn==1.5.0" \
"geopandas==0.14.4" "numpy==1.26.4" "pandas==2.2.2" \
"scipy==1.13.1" "matplotlib==3.9.0"GDAL/PROJ system dependencies (Ubuntu/Debian):
sudo apt-get install -y gdal-bin libgdal-dev libproj-devYou need a fitted model, a held-out feature frame, and — critically — the geometry or coordinates for every held-out row. Explanations that cannot be joined back to a location are half an analysis.
Step-by-Step Implementation
Step 1 — Keep identifiers alongside the design matrix
import geopandas as gpd
import numpy as np
import pandas as pd
FEATURES = ["elevation", "slope_deg", "northness", "tpi_11",
"ndvi_mean", "impervious_pct", "rain_24h", "upslope_area"]
samples = gpd.read_file("samples.gpkg") # one row per observation
X = samples[FEATURES].astype("float32")
y = samples["flooded"].astype("int8")
coords = np.c_[samples.geometry.x.values, samples.geometry.y.values]
assert not X.isna().any().any(), "impute or drop missing values before explaining"
assert len(X) == len(coords) == len(y)Step 2 — Assign spatial blocks and fit under a blocked split
from xgboost import XGBClassifier
def block_id(coords: np.ndarray, block_m: float = 5000.0) -> np.ndarray:
"""Integer id of the square block each point falls in."""
ix = np.floor(coords[:, 0] / block_m).astype("int64")
iy = np.floor(coords[:, 1] / block_m).astype("int64")
return ix * 100_000 + iy
blocks = block_id(coords, block_m=5000.0)
rng = np.random.default_rng(7)
uniq = np.unique(blocks)
rng.shuffle(uniq)
holdout_blocks = set(uniq[: max(1, len(uniq) // 4)].tolist())
test_mask = np.isin(blocks, list(holdout_blocks))
model = XGBClassifier(
n_estimators=400, max_depth=5, learning_rate=0.05,
subsample=0.8, colsample_bytree=0.8, eval_metric="logloss",
random_state=0,
)
model.fit(X[~test_mask], y[~test_mask])
print("held-out AUC:", model.score(X[test_mask], y[test_mask]))Step 3 — Block permutation importance
from sklearn.metrics import roc_auc_score
def block_permutation_importance(model, X, y, blocks, n_repeats=10, seed=0):
"""Permute each feature WITHIN spatial blocks, so shuffled values stay plausible."""
rng = np.random.default_rng(seed)
base = roc_auc_score(y, model.predict_proba(X)[:, 1])
order = np.argsort(blocks, kind="stable")
starts = np.flatnonzero(np.r_[True, np.diff(blocks[order]) != 0])
groups = np.split(order, starts[1:])
rows = []
for col in X.columns:
drops = []
for _ in range(n_repeats):
Xp = X.copy()
values = Xp[col].to_numpy()
for g in groups:
values[g] = rng.permutation(values[g])
Xp[col] = values
drops.append(base - roc_auc_score(y, model.predict_proba(Xp)[:, 1]))
rows.append({"feature": col, "auc_drop": float(np.mean(drops)),
"std": float(np.std(drops))})
return pd.DataFrame(rows).sort_values("auc_drop", ascending=False)
imp = block_permutation_importance(model, X[test_mask], y[test_mask], blocks[test_mask])
print(imp.to_string(index=False))
assert imp["auc_drop"].max() > 0, "no feature affects the score — check the model is fitted"Comparing this table against sklearn.inspection.permutation_importance on the same model is a diagnostic in itself: features whose rank drops sharply when permutation is blocked were being credited for spatial smoothness rather than for signal.
Step 4 — SHAP values, joined back to geometry
import shap
explainer = shap.TreeExplainer(model)
sv = explainer.shap_values(X[test_mask])
if isinstance(sv, list): # binary classifiers return one array per class
sv = sv[1]
shap_df = pd.DataFrame(sv, columns=[f"shap_{c}" for c in FEATURES])
explained = gpd.GeoDataFrame(
pd.concat([samples.loc[test_mask, ["sample_id"]].reset_index(drop=True), shap_df], axis=1),
geometry=samples.loc[test_mask, "geometry"].reset_index(drop=True),
crs=samples.crs,
)
explained.to_file("shap_attributions.gpkg", driver="GPKG")
# Local attributions must reconstruct the model output, up to the base value.
recon = explainer.expected_value + sv.sum(axis=1)
raw = model.predict(X[test_mask], output_margin=True)
assert np.allclose(recon, raw, atol=1e-3), "SHAP values do not sum to the model margin"That assertion is the one check worth keeping permanently: if attributions do not reconstruct the margin, the explainer was pointed at a different model, a different feature order, or a transformed matrix.
Step 5 — Grouped importance for collinear features
from scipy.cluster.hierarchy import fcluster, linkage
from scipy.spatial.distance import squareform
def correlation_groups(X: pd.DataFrame, threshold: float = 0.3) -> dict:
"""Cluster features by 1 - |Spearman rho| and return group membership."""
rho = X.corr(method="spearman").abs().to_numpy()
dist = 1.0 - rho
np.fill_diagonal(dist, 0.0)
dist = (dist + dist.T) / 2.0
z = linkage(squareform(dist, checks=False), method="average")
labels = fcluster(z, t=threshold, criterion="distance")
groups = {}
for name, lab in zip(X.columns, labels):
groups.setdefault(int(lab), []).append(name)
return groups
groups = correlation_groups(X[test_mask])
for gid, members in groups.items():
total = float(np.abs(sv[:, [FEATURES.index(m) for m in members]]).sum(axis=1).mean())
print(f"group {gid}: {members} -> mean |SHAP| {total:.4f}")Verification & Testing
Explanations deserve tests as much as models do, because a silently mis-joined SHAP matrix produces confident, entirely wrong maps.
def test_shap_row_order_matches_samples():
"""Attribution i must belong to sample i — the join is positional and fragile."""
sub = X[test_mask].reset_index(drop=True)
idx = 17
one = explainer.shap_values(sub.iloc[[idx]])
one = one[1] if isinstance(one, list) else one
assert np.allclose(one[0], sv[idx], atol=1e-4)
def test_constant_feature_gets_zero_attribution():
"""A feature with no variance cannot move any prediction."""
Xc = X[test_mask].copy()
Xc["rain_24h"] = Xc["rain_24h"].mean()
m2 = XGBClassifier(n_estimators=50, max_depth=3, random_state=0).fit(Xc, y[test_mask])
s2 = shap.TreeExplainer(m2).shap_values(Xc)
s2 = s2[1] if isinstance(s2, list) else s2
assert np.abs(s2[:, list(Xc.columns).index("rain_24h")]).max() < 1e-6Beyond code, run the residual map check: render the held-out residuals and the attribution maps side by side. If a region has large residuals and no feature shows notable attribution there, the model is missing a driver in that region — the most actionable output explainability produces.
Interpreting What You Get Back
An attribution is a statement about the model, and translating it into a statement about the world takes a further step that is easy to skip. Three habits make that translation safer.
Read attributions against the base value, not in isolation. A SHAP value of +0.4 means the feature pushed this sample four tenths of a logit above the model’s average prediction for the dataset. If the base value already sits near a decision boundary, that push may flip the class; if it sits far from one, the same push changes nothing. Reporting attributions without the base value invites readers to treat a small nudge as a decisive cause.
Distinguish a driver from a proxy. Geospatial feature sets are dense with proxies: road density stands in for population, elevation stands in for temperature and land use together, and distance to a river stands in for both flood exposure and soil type. The model has no way to prefer the cause over the proxy, and neither does the attribution. When a proxy tops the ranking, the useful next step is usually to add the thing it is proxying for and see whether the proxy’s importance collapses — which is both a modelling improvement and a much stronger piece of evidence than any single explanation run.
Treat a regional disagreement as a hypothesis, not a conclusion. If a feature dominates in one region and vanishes in another, there are three ordinary explanations before any interesting one: the feature has almost no variance in the second region, the second region has too few samples for the estimate to be stable, or the two regions were processed with different upstream settings. Check all three before concluding that the process differs. The residual-versus-importance matrix earlier in this page is the fastest way to triage which explanation applies.
Finally, be explicit about what an explanation cannot do. It cannot tell you what would happen if a feature changed, because the model was never asked that question — attribution decomposes a prediction, it does not simulate an intervention. If the decision at hand is “what happens if we build this drainage channel”, the honest answer requires either a causal design or a physical model, and an attribution map is at best a way of generating candidate hypotheses to test with one. Saying so plainly, in the same document as the maps, is what keeps an explainability report useful rather than persuasive.
Troubleshooting & Common Errors
shap.TreeExplainer returns a list, and later code fails on .shape — binary classifiers return one array per class. Take index 1 (positive class) explicitly rather than relying on the shape, which changed across SHAP releases.
Attributions do not sum to the prediction — the explainer saw a different feature order than the model, or the model was fitted on a NumPy array while SHAP got a DataFrame with reordered columns. Always pass a DataFrame with the same column order in both places, and keep the reconstruction assertion from Step 4.
MemoryError computing SHAP on a million rows — SHAP output is n_samples × n_features in float64. Explain a stratified sample of 20,000 held-out rows; the maps look identical and the memory drops by 50×.
Importance ranking changes every run — permutation with too few repeats. Ten repeats with a fixed seed and a reported standard deviation is the minimum publishable configuration.
Every attribution map looks like the elevation map — a symptom of collinearity plus a smooth dominant covariate. Run the grouped importance from Step 5 before drawing conclusions, and consider the decorrelation approaches in dimensionality reduction for spatial data.
Explanations differ between the training environment and production — the served model is not the explained model. Explain the exact exported artifact, and verify it as part of the ONNX export for geospatial model inference checks.
Performance Optimisation
Sample, do not subsample randomly. Stratify the explanation sample by predicted class and by spatial block, so rare classes and remote regions are represented. A 20,000-row stratified sample gives stable global rankings and usable maps at a tiny fraction of the cost.
Prefer TreeExplainer for tree models. It is exact and polynomial-time, whereas KernelExplainer is a sampling approximation that costs thousands of model evaluations per row. Reserve kernel or permutation explainers for models with no tree structure.
Cache the SHAP matrix as a data artifact. It is expensive, deterministic given a model and a sample, and useful to many downstream consumers — dashboards, drift reports, model cards. Version it next to the model with spatial dataset versioning with DVC and lakeFS.
Parallelise permutation, not SHAP. Block permutation importance is n_features × n_repeats independent model evaluations and scales linearly across processes; TreeExplainer is already multi-threaded internally and gains little from an outer pool.
Finally, wire the attribution summary into monitoring. A shift in the ranking of features between the training sample and a production batch is an early, interpretable signal of covariate shift — often visible before the score itself moves, which is why it pairs naturally with model drift detection for geospatial inference.
FAQ
Why do coordinates always come out as the most important features?
Because raw x and y let a tree memorise the target surface, absorbing credit that belongs to physical drivers. Under a random split this looks like strong signal; under a blocked split it collapses. Either drop raw coordinates, or keep them and evaluate exclusively with spatial folds.
Is permutation importance valid for spatially autocorrelated data?
Not in its standard form. Global shuffling destroys spatial structure and manufactures impossible feature combinations, which inflates the apparent importance of smooth surfaces such as elevation. Permute within spatial blocks, as in Step 3, and compare the two rankings — the difference is itself informative.
What does a map of SHAP values actually show?
Where a feature pushed the prediction above or below the model’s base value. Coherent patches of positive attribution usually reflect a real regional process; salt-and-pepper attribution usually means the feature is behaving as noise or proxying for something unmodelled.
Do SHAP values prove causation?
No — they decompose the model, not the world. A proxy feature (road density standing in for population) will be credited as if it were the cause. Treat attributions as hypotheses about the model, then test them against domain knowledge.
Related
- Computing SHAP Values for a Land Cover Classifier — the multiclass workflow end to end
- Mapping Permutation Importance Across Space — per-region importance surfaces
- Spatial Cross-Validation Strategies — the split design every explanation depends on
- Gradient Boosting for Raster Data — the models most commonly explained here
- Dimensionality Reduction for Spatial Data — handling the collinearity that dilutes attributions