Computing SHAP Values for a Land Cover Classifier

Compute and map SHAP values for a multiclass land-cover model: TreeExplainer output shapes, per-class attributions, reconstruction assertions and joining values back to geometry.

Use shap.TreeExplainer(model) on a boosted or random-forest land-cover classifier, normalise the output to a (n_samples, n_features, n_classes) array, assert that the attributions reconstruct the model margin, and keep the sample identifiers so every row can be joined back to its geometry. The multiclass case is where most implementations quietly go wrong, because the output shape has changed across SHAP releases and an incorrect axis assumption produces attributions that look plausible and describe the wrong class.

This page is the multiclass workflow. For the surrounding methodology — block permutation, collinearity, mapping — see Model Explainability for Spatial Predictions.

The multiclass output is a cube; you must choose how to reduce it The SHAP output for a multiclass model is drawn as a three dimensional array indexed by sample, feature and class. Two reductions are shown: slicing a single class for every sample, which answers what would push each sample towards that class, and gathering each sample's own predicted class, which answers why the model chose what it chose. shap_values shape: (samples, features, classes) samples × features one layer per class classes → sv[:, :, k] fixed class k for all samples “what pushes towards water?” gather on predicted class per-sample argmax layer “why did it choose this?” 2-D matrix samples × features joinable to geometry Indexing the wrong axis gives a plausible map of the wrong class.

Why This Fails in Geospatial ML Pipelines

The output-shape ambiguity is the first trap and it is entirely avoidable. Different SHAP versions have returned a list of per-class arrays and a stacked three-dimensional array; XGBoost, LightGBM and scikit-learn wrappers have each behaved slightly differently. Code written against one combination indexes the wrong axis under another, and because the resulting matrix still has the right shape and plausible magnitudes, nothing complains. The output is a confident explanation of a class the analyst did not intend.

The second trap is losing the join. SHAP takes a matrix and returns a matrix; geometry is not carried through. If the feature frame was filtered, re-sorted or de-duplicated between the model fit and the explanation, row i of the attribution matrix no longer corresponds to row i of the GeoDataFrame — and the resulting map is a spatially scrambled but visually smooth surface, which is worse than no map at all.

Third, memory. Attributions are n_samples × n_features × n_classes in float64. A million pixels, twenty features and eight classes is 1.3 GB before any copies. Teams discover this at the end of a long run, so the sensible move is to explain a stratified sample from the start rather than to retrofit sampling after the first MemoryError.

Core Principles

  • Normalise the output shape immediately and assert it.
  • Choose the reduction deliberately — fixed class or predicted class — and name the output for it.
  • Assert reconstruction. Base value plus attributions must equal the raw margin.
  • Carry identifiers, not just arrays. Explain a frame that keeps its index.
  • Stratify the sample by predicted class and spatial block.
  • Version the attribution matrix as an artifact next to the model.

Production-Ready Code

from __future__ import annotations

import logging
import numpy as np
import pandas as pd
import shap

logger = logging.getLogger(__name__)


def explain_multiclass(model, x: pd.DataFrame) -> tuple[np.ndarray, np.ndarray]:
    """Return SHAP values normalised to (n_samples, n_features, n_classes).

    Handles both the list-of-arrays and stacked-array conventions that different
    SHAP releases use, so downstream code can index one known layout.
    """
    explainer = shap.TreeExplainer(model)
    raw = explainer.shap_values(x)

    if isinstance(raw, list):                       # older convention
        sv = np.stack(raw, axis=-1)
    else:
        sv = np.asarray(raw)
        if sv.ndim == 2:                            # binary model
            sv = np.stack([-sv, sv], axis=-1)

    n_samples, n_features = x.shape
    if sv.shape[:2] != (n_samples, n_features):
        raise ValueError(f"unexpected SHAP shape {sv.shape} for input {x.shape}")

    base = np.atleast_1d(np.asarray(explainer.expected_value))
    if base.shape[0] != sv.shape[2]:
        base = np.repeat(base, sv.shape[2])[: sv.shape[2]]

    logger.info("SHAP cube %s, %d class base values", sv.shape, base.shape[0])
    return sv, base


def assert_reconstructs(model, x: pd.DataFrame, sv: np.ndarray, base: np.ndarray,
                        atol: float = 1e-2) -> None:
    """Base value plus contributions must equal the model's raw margin, per class."""
    margins = model.predict(x, output_margin=True)
    margins = np.atleast_2d(margins)
    if margins.shape[0] == 1 and sv.shape[0] != 1:
        margins = margins.T
    recon = base[None, :] + sv.sum(axis=1)
    if not np.allclose(recon, margins, atol=atol):
        worst = float(np.max(np.abs(recon - margins)))
        raise ValueError(
            f"SHAP values do not reconstruct the margin (max error {worst:.4f}). "
            "The explainer and the model likely disagree on feature order.")


def attributions_for_prediction(sv: np.ndarray, predicted: np.ndarray) -> np.ndarray:
    """Gather each sample's own predicted-class layer -> (n_samples, n_features)."""
    idx = np.asarray(predicted, dtype="int64")
    if idx.shape[0] != sv.shape[0]:
        raise ValueError("predicted class vector does not match the sample count")
    return np.take_along_axis(sv, idx[:, None, None], axis=2)[:, :, 0]


def stratified_sample(frame: pd.DataFrame, predicted: np.ndarray, blocks: np.ndarray,
                      n: int = 20_000, seed: int = 0) -> np.ndarray:
    """Row positions stratified by predicted class and spatial block."""
    rng = np.random.default_rng(seed)
    key = pd.Series([f"{c}_{b}" for c, b in zip(predicted, blocks)])
    per_group = max(1, n // key.nunique())
    picks = []
    for _, positions in key.groupby(key).groups.items():
        pos = np.asarray(positions)
        picks.append(rng.choice(pos, size=min(per_group, len(pos)), replace=False))
    return np.concatenate(picks)

Step-by-Step Walkthrough

Step 1 — Take a stratified sample

import numpy as np

predicted = model.predict(X_test)
blocks = block_id(coords_test, block_m=5000.0)
rows = stratified_sample(X_test, predicted, blocks, n=20_000, seed=0)

X_sample = X_test.iloc[rows]
pred_sample = predicted[rows]
print(f"explaining {len(rows):,} rows over {len(np.unique(blocks[rows]))} blocks")

Step 2 — Explain and verify in the same breath

sv, base = explain_multiclass(model, X_sample)
assert_reconstructs(model, X_sample, sv, base)
print("cube:", sv.shape, "base values:", np.round(base, 3))

If assert_reconstructs raises, stop. Every downstream map would be wrong, and the usual cause is a feature-order mismatch between the frame passed to fit and the one passed to the explainer.

Step 3 — Reduce to a per-sample matrix

attrib = attributions_for_prediction(sv, pred_sample)     # (n, n_features)

CLASSES = ["background", "cropland", "forest", "water", "built"]
per_class = {
    name: np.abs(sv[:, :, k]).mean(axis=0)
    for k, name in enumerate(CLASSES)
}
import pandas as pd
print(pd.DataFrame(per_class, index=X_sample.columns).round(4))

The per-class table is often the most informative single output: it shows that elevation drives forest, that shortwave infrared drives water, and that the built class rests on a feature nobody expected — which is usually where the modelling insight is.

Step 4 — Join back to geometry and write a mappable artifact

import geopandas as gpd

explained = gpd.GeoDataFrame(
    pd.DataFrame(attrib, columns=[f"shap_{c}" for c in X_sample.columns],
                 index=X_sample.index),
    geometry=samples.loc[X_sample.index, "geometry"],
    crs=samples.crs,
)
explained["predicted"] = [CLASSES[i] for i in pred_sample]
explained.to_parquet("artifacts/shap_attributions.parquet")
assert len(explained) == len(X_sample)

Using the frame’s own index for the join — rather than positional alignment — is what makes this safe against any upstream filtering.

A per-class table says more than one global ranking A table of mean absolute SHAP values with features as rows and land cover classes as columns. Shortwave infrared dominates the water class, elevation dominates forest, and near infrared dominates cropland, so no single global ranking describes the model. Mean |SHAP| by feature and class featurecroplandforest waterbuilt nir0.410.220.180.09 swir20.120.140.530.16 elevation0.080.370.050.11 tpi_110.060.090.210.04 impervious0.050.030.070.44 Averaging these columns into one ranking would hide every one of these relationships.

Verification

import numpy as np
import pandas as pd
import pytest


def test_output_is_normalised_to_a_cube():
    sv, base = explain_multiclass(model, X_sample.head(50))
    assert sv.ndim == 3
    assert sv.shape[:2] == (50, X_sample.shape[1])
    assert base.shape[0] == sv.shape[2]


def test_reconstruction_holds():
    sv, base = explain_multiclass(model, X_sample.head(200))
    assert_reconstructs(model, X_sample.head(200), sv, base)


def test_feature_order_mismatch_is_detected():
    shuffled = X_sample.head(100)[list(reversed(X_sample.columns))]
    sv, base = explain_multiclass(model, shuffled)
    with pytest.raises(ValueError, match="reconstruct"):
        assert_reconstructs(model, shuffled, sv, base)


def test_gather_picks_the_predicted_layer():
    sv = np.arange(2 * 3 * 4, dtype="float64").reshape(2, 3, 4)
    predicted = np.array([1, 3])
    out = attributions_for_prediction(sv, predicted)
    assert np.array_equal(out[0], sv[0, :, 1])
    assert np.array_equal(out[1], sv[1, :, 3])


def test_join_survives_a_filtered_frame():
    """Filtering the frame before explaining must not scramble the geometry join."""
    sub = X_test[X_test["elevation"] > 500]
    sv, base = explain_multiclass(model, sub.head(100))
    attrib = attributions_for_prediction(sv, model.predict(sub.head(100)))
    frame = pd.DataFrame(attrib, index=sub.head(100).index)
    assert frame.index.equals(sub.head(100).index)

The feature-order test is worth writing once and keeping forever: it is the failure that produces the most convincing wrong answers.

Join on the index, never on position A feature frame is filtered before explanation, producing an attribution matrix whose row order matches the filtered frame rather than the original. Joining positionally to the original geometry scrambles the map; joining on the retained index produces the correct assignment. SHAP returns a bare matrix — the index is yours to preserve samples (10 000 rows) index 0 … 9 999 geometry attached filter filtered (3 200 rows) index 7, 12, 30, … order preserved shap matrix (3 200 rows) positions 0 … 3 199 no index of its own Positional join back to the 10 000-row frame → every attribution lands on the wrong parcel. Rebuild the DataFrame with index=filtered.index and the join is correct by construction.

FAQ

What shape does TreeExplainer return for a multiclass model?

It varies by release and model wrapper — a stacked (samples, features, classes) array in recent versions, a list of per-class arrays in older ones. Normalise to one layout immediately, as explain_multiclass does, and assert the shape rather than indexing on assumption.

Which class should I explain?

The predicted class per sample, if the question is “why did the model choose this”. Explaining a fixed class answers “what would push each pixel towards that class” — a legitimate but different question. Name the output after whichever you chose.

How large a sample do I need for stable global rankings?

A few tens of thousands of rows, stratified by predicted class and spatial block. Beyond that the ranking barely moves while memory doubles, because the cube is samples × features × classes in float64.

Part of: Model Explainability for Spatial Predictions Part of: Training Geospatial Predictive Models in Python