🧭 Concept Spatial machine learning Spatial Machine Learning · step 12 Notebook By · Updated

Extrapolation in Space Explained: The Area of Applicability

On this page

Problem statement

A prediction map covers every pixel. Some of those pixels resemble the training data closely; others do not resemble it at all, and the model is guessing there with no way to say so.

Measured on a real problem — Snowdon NDVI from terrain — a model trained on the western half of a 10 km window and applied to the eastern half:

random cross-validation on the training half   +0.338
block cross-validation on the training half    +0.033
actual R² on the eastern half                  -0.409

Both internal estimates were optimistic, by 0.75 and 0.44 R², and the true result was worse than predicting the training mean.

No cross-validation scheme fixes this, because the training data contains no information about the new area. What you can do is say where the model is extrapolating, which is what an area of applicability is.

Quick answer

Measure distance in scaled feature space from each prediction location to the training data, and threshold it against the distances within the training set:

import numpy as np
from scipy.spatial import cKDTree
from sklearn.preprocessing import StandardScaler


def applicability(X_train, X_new, quantile=0.95, k=5):
    scaler = StandardScaler().fit(X_train)
    train = scaler.transform(X_train)
    new = scaler.transform(X_new)

    tree = cKDTree(train)
    within, _ = tree.query(train, k=k + 1)        # k+1 because self is nearest
    threshold = np.quantile(within[:, 1:].mean(axis=1), quantile)

    outside, _ = tree.query(new, k=k)
    distance = outside.mean(axis=1)
    return distance <= threshold, distance, threshold


inside, distance, threshold = applicability(X_train, X_new)
print(f"{inside.mean():.1%} of the eastern half is inside the area learnt in the west")
74.5% of the eastern half is inside the area learnt in the west

Ship the resulting mask with the prediction. It does not estimate the error outside; it says where the covariates are unlike anything the model saw — 25.5% of the eastern half here. It cannot say whether the relationship still holds inside: in this case it did not (Example 1).

Prediction locations inside and outside the training data feature-space envelope, with the outside region flagged rather than assigned an error. where the covariates resemble the training datafeature spacered = outside the envelopetwo covariates shownThe mask says where the error is unknown — it does not estimate it spatialworkflow.io
The mask is about feature space, not geography. A nearby pixel can be outside it and a distant one inside.

Feature space not map space, scaled, thresholded, shipped

1. Distinguish geographic distance from feature-space distance

A pixel a hundred kilometres away with the same slope, elevation and aspect is well inside the model's competence. A pixel next door on a different geology may be far outside it.

Applicability is about the covariates, not the coordinates. That is why the mask often looks nothing like a buffer around the training points.

2. Scale the features before measuring distance

Elevation in metres ranges over a thousand; northness ranges over two. Unscaled, elevation dominates the distance entirely.

StandardScaler puts each feature on a comparable footing. Weighting by feature importance is a refinement: features the model barely uses should contribute less to the distance.

3. Derive the threshold from the training data

The useful comparison is against how far apart the training points themselves are. If a new location is further from the training data than 95% of training points are from their own neighbours, it is outside the envelope the model learned.

That makes the threshold self-calibrating: a dense training set gives a tight envelope, a sparse one a generous one.

4. Check the univariate ranges too

Multivariate distance can hide a single feature far outside its range, because the other features compensate.

lo, hi = np.percentile(X_train[:, i], [1, 99])
outside = ((X_new[:, i] < lo) | (X_new[:, i] > hi)).mean()

A feature with more than about 10% of its new values outside the training band is extrapolation regardless of what the distance says. For the eastern half, elevation was 12.0% outside the band learnt in the west (Example 2).

5. Ship the mask, not a corrected prediction

There is no defensible correction for extrapolation. The honest output is two rasters: the prediction, and a mask of where it is supported.

A nearby pixel outside the feature-space envelope and a distant pixel inside it, showing that applicability is not a buffer. applicability follows the covariates, not the mapgeographic distancefeature-space distanceapplicable?adjacent quarry✓near✕far✕novalley 100 km away✕far✓near✓yesThe model does not see coordinates unless you gave it coordinates spatialworkflow.io
Applicability follows the covariates. Distance on the map is a poor proxy for it.

An importance-weighted mask and a per-feature report

The data used on this page

A 10 km window on Snowdon, 4.16–4.00° W, 53.03–53.12° N. The target is NDVI from the least cloudy Sentinel-2 L2A scene of June–September 2025 on Earth Search, on the scene's 20 m grid; the features are terrain from the Copernicus DEM GLO-30 resampled to the same grid — elevation, slope, northness, eastness and a topographic position index over 300 m. Cells the scene classifies as water, cloud, cloud shadow, snow or no data are dropped. X, y and coords hold every remaining cell. The deployment the page audits is a model trained on 3,000 random cells of the western half (X_train, y_train) and applied to every cell of the eastern half (X_new, y_new).

import numpy as np
import rasterio
import rasterio.merge
from pystac_client import Client
from rasterio.enums import Resampling
from rasterio.warp import reproject, transform_bounds
from scipy import ndimage

AOI = (-4.16, 53.03, -4.00, 53.12)
items = Client.open("https://earth-search.aws.element84.com/v1").search(
    collections=["sentinel-2-l2a"], bbox=list(AOI), datetime="2025-06-01/2025-10-01",
    query={"eo:cloud_cover": {"lt": 60}}).items()
item = min(items, key=lambda i: i.properties["eo:cloud_cover"])

with rasterio.open(item.assets["scl"].href) as ds:                 # the 20 m grid
    win = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", ds.crs, *AOI),
                                       transform=ds.transform).round_offsets().round_lengths()
    scl, transform, crs = ds.read(1, window=win), ds.window_transform(win), ds.crs
    bounds = rasterio.windows.bounds(win, ds.transform)
shape = scl.shape
band = {}
for name in ("red", "nir"):                                         # 10 m bands averaged onto it
    with rasterio.open(item.assets[name].href) as ds:
        band[name] = ds.read(1, window=rasterio.windows.from_bounds(*bounds, transform=ds.transform),
                             out_shape=shape, resampling=Resampling.average).astype("float32") / 1e4
ndvi = (band["nir"] - band["red"]) / (band["nir"] + band["red"])

tiles = [rasterio.open(f"https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{t}_00_DEM/"
                       f"Copernicus_DSM_COG_10_N53_00_{t}_00_DEM.tif") for t in ("W005", "W004")]
mosaic, mosaic_transform = rasterio.merge.merge(tiles, bounds=(-4.18, 53.02, -3.98, 53.13))
elev = np.full(shape, np.nan, dtype="float32")
reproject(mosaic[0], elev, src_transform=mosaic_transform, src_crs=tiles[0].crs,
          dst_transform=transform, dst_crs=crs, resampling=Resampling.bilinear)
dz_dy, dz_dx = np.gradient(elev, -20, 20)                           # rows run south
slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))
aspect = np.arctan2(-dz_dx, -dz_dy)                                 # direction the slope faces, from north
tpi = elev - ndimage.uniform_filter(elev, size=15)                  # against the 300 m neighbourhood

valid = np.isin(scl, (2, 4, 5, 7)) & np.isfinite(ndvi) & np.isfinite(elev)
rows, cols = np.nonzero(valid)
xs, ys = rasterio.transform.xy(transform, rows, cols)
coords = np.column_stack([xs, ys])
X = np.column_stack([elev[valid], slope[valid], np.cos(aspect[valid]), np.sin(aspect[valid]), tpi[valid]])
y = ndvi[valid]
FEATURES = ["elevation", "slope", "northness", "eastness", "tpi"]
west = coords[:, 0] < np.median(coords[:, 0])
pick = np.random.default_rng(0).choice(np.nonzero(west)[0], 3000, replace=False)
X_train, y_train, X_new, y_new = X[pick], y[pick], X[~west], y[~west]
feature_names = FEATURES
print(f"{item.id}: {shape[1]} x {shape[0]} cells of 20 m, {valid.sum():,} usable "
      f"({valid.mean():.1%}); {len(X_train):,} training cells in the west, {len(X_new):,} to predict in the east")
S2C_30UVD_20250712_0_L2A: 543 x 509 cells of 20 m, 270,536 usable (97.9%); 3,000 training cells in the west, 135,292 to predict in the east

Example 1 — an importance-weighted applicability mask

import numpy as np
from scipy.spatial import cKDTree
from sklearn.preprocessing import StandardScaler


def applicability_mask(X_train, X_new, importances=None, quantile=0.95, k=5):
    """Feature-space distance, weighted by how much the model uses each feature."""
    scaler = StandardScaler().fit(X_train)
    train = scaler.transform(X_train)
    new = scaler.transform(X_new)

    if importances is not None:
        weights = np.asarray(importances, dtype=float)
        weights = weights / weights.sum() * len(weights)
        train = train * np.sqrt(weights)
        new = new * np.sqrt(weights)

    tree = cKDTree(train)
    within, _ = tree.query(train, k=min(k + 1, len(train)))
    reference = within[:, 1:].mean(axis=1)
    threshold = float(np.quantile(reference, quantile))

    distances, _ = tree.query(new, k=min(k, len(train)))
    mean_distance = distances.mean(axis=1)
    inside = mean_distance <= threshold

    print(f"  training neighbour distance: median {np.median(reference):.3f}, "
          f"p{quantile * 100:.0f} {threshold:.3f}")
    print(f"  deployment distance: median {np.median(mean_distance):.3f}, "
          f"p95 {np.quantile(mean_distance, 0.95):.3f}")
    print(f"  inside the area of applicability: {inside.mean():.1%}")
    return inside, mean_distance, threshold


from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score

forest = RandomForestRegressor(n_estimators=200, min_samples_leaf=2, n_jobs=-1, random_state=0).fit(X_train, y_train)
predicted = forest.predict(X_new)
for label, weights in (("unweighted", None), ("importance-weighted", forest.feature_importances_)):
    print(label)
    mask, _, _ = applicability_mask(X_train, X_new, importances=weights)
    for side, part in (("inside", mask), ("outside", ~mask)):
        print(f"    {side:7}: median error {np.median(np.abs(y_new[part] - predicted[part])):.3f}, "
              f"R2 {r2_score(y_new[part], predicted[part]):+.3f}")
unweighted
  training neighbour distance: median 0.404, p95 0.946
  deployment distance: median 0.657, p95 1.603
  inside the area of applicability: 74.5%
    inside : median error 0.067, R2 -0.891
    outside: median error 0.117, R2 -0.342
importance-weighted
  training neighbour distance: median 0.400, p95 0.925
  deployment distance: median 0.656, p95 1.603
  inside the area of applicability: 73.3%
    inside : median error 0.066, R2 -0.887
    outside: median error 0.118, R2 -0.363

Weighting by importance is the principled refinement: without it, a feature the model barely uses contributes as much to the distance as the one it depends on. Here it made almost no difference — 73.3% inside against 74.5% — because elevation, the feature the forest leant on, was also the one that shifted.

The more important result is on either side of the line. Inside the mask the errors were smaller, a median of 0.067 against 0.117 outside, so the mask did find the worse predictions. But R² was negative on both sides (−0.891 inside, −0.342 outside, each against its own variance): the model was worse than the mean everywhere in the east. The mask marks unfamiliar covariates; a relationship that changes between regions while the covariates stay familiar is invisible to it.

The area of applicability of a model trained in the west of the Snowdon window, unweighted and weighted by feature importance, with the errors inside and outside each. east half, trained in the west: two masksmaskinsideerror insideerror outsideunweighted74.5%0.067 (R² −0.891)0.117 (R² −0.342)importance-weighted73.3%0.066 (R² −0.887)0.118 (R² −0.363)errors are median absolute NDVI; R² is against each side’s own varianceWeighting barely moved the line; neither side of it was a working model spatialworkflow.io
Measured in Example 1. The mask separates larger errors from smaller ones; it cannot certify the smaller ones.

Example 2 — a per-feature extrapolation report

import numpy as np


def extrapolation_report(X_train, X_new, feature_names, band=(1, 99)):
    """Which features are being extrapolated, and by how much."""
    rows = []
    for i, name in enumerate(feature_names):
        train, new = X_train[:, i], X_new[:, i]
        lo, hi = np.percentile(train, band)
        below = float((new < lo).mean())
        above = float((new > hi).mean())
        span = hi - lo
        reach = float(max((lo - new.min()) / span if new.min() < lo else 0,
                          (new.max() - hi) / span if new.max() > hi else 0))

        rows.append({"feature": name, "below": below, "above": above,
                     "reach": reach})
        flag = "  <- extrapolating" if below + above > 0.1 else ""
        print(f"  {name:14} below {below:6.1%}  above {above:6.1%}  "
              f"reach {reach:5.2f} training ranges{flag}")

    worst = max(rows, key=lambda r: r["below"] + r["above"])
    print(f"  worst: {worst['feature']} "
          f"({worst['below'] + worst['above']:.1%} outside)")
    return rows


extrapolation_report(X_train, X_new, feature_names)
  elevation      below   5.5%  above   6.5%  reach  0.34 training ranges  <- extrapolating
  slope          below   0.7%  above   2.6%  reach  0.48 training ranges
  northness      below   1.1%  above   1.2%  reach  0.00 training ranges
  eastness       below   1.0%  above   0.7%  reach  0.00 training ranges
  tpi            below   3.2%  above   2.1%  reach  0.83 training ranges
  worst: elevation (12.0% outside)

The reach column — how far beyond the training range the new data goes, in units of the training range — says how far outside the model is being pushed. A reach of 0.1 is a mild extension; a reach of 2 means the model has never seen anything like it. Here the worst was TPI at 0.83, though only 5.3% of the eastern cells were outside its band; elevation had the most cells outside, 12.0%, at a reach of 0.34.

Example 3 — writing the prediction and the mask together

import numpy as np
import rasterio


def write_prediction(prediction, inside, distance, profile, out_path):
    """Three bands: the prediction, the mask, and the distance behind it."""
    profile = dict(profile) | {
        "count": 3, "dtype": "float32", "nodata": np.nan,
        "compress": "deflate", "tiled": True,
    }
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(prediction.astype("float32"), 1)
        dst.write(inside.astype("float32"), 2)
        dst.write(distance.astype("float32"), 3)
        dst.set_band_description(1, "prediction")
        dst.set_band_description(2, "1 = inside the area of applicability")
        dst.set_band_description(3, "feature-space distance to training data")
        dst.update_tags(
            applicability_note="Band 1 is only supported where band 2 is 1. "
                               "Outside that, the model is extrapolating and "
                               "the error is unknown.")

    print(f"  {np.nanmean(inside):.1%} of pixels are inside the area of "
          f"applicability")                    # nanmean: nodata cells are NaN
    print(f"  wrote {out_path} with the mask as band 2")
    return out_path


grid = {name: np.full(shape, np.nan, dtype="float32") for name in ("prediction", "inside", "distance")}
cell_inside, cell_distance, _ = applicability(X_train, X)
grid["prediction"][valid] = forest.predict(X)
grid["inside"][valid] = cell_inside
grid["distance"][valid] = cell_distance
write_prediction(grid["prediction"], grid["inside"], grid["distance"],
                 {"driver": "GTiff", "width": shape[1], "height": shape[0], "crs": crs, "transform": transform},
                 "prediction.tif")
  84.8% of pixels are inside the area of applicability
  wrote prediction.tif with the mask as band 2

Writing the note into the file's tags is worth the line. The mask travels with the prediction, but the interpretation of the mask usually does not, and a raster tag survives every copy.

Why cross-validation cannot see extrapolation

Why cross-validation cannot estimate extrapolation error

Cross-validation resamples the training data. Every fold is drawn from the same distribution, so the estimate describes generalisation within that distribution.

A new region is a different distribution. The training data contains no observations of it, so no resampling of that data can say how the model behaves there.

The measured case is unambiguous: block CV — the scheme specifically recommended for spatial data — promised +0.033 where the actual transfer gave −0.409. It was less optimistic than random CV and still wrong by 0.44 R².

Why feature space, not geographic space

The model does not see coordinates unless you gave it coordinates. It sees slope, elevation, reflectance, distance to water.

So its competence extends to combinations of those it has seen, wherever they occur. A distant valley with familiar terrain is inside; an adjacent quarry with an unfamiliar combination is outside.

That is why an applicability mask often looks patchy rather than like a buffer, and why it is more useful than a distance-from-training-points map.

Why the threshold comes from the training data

An absolute distance threshold means nothing, because feature-space distances depend on the features, their scaling and their number.

Comparing against the within-training neighbour distances makes it relative and self-calibrating. If the training set is dense in feature space, the envelope is tight. If it is sparse, the envelope is generous, correctly reflecting that the model was fitted on scattered evidence anyway.

The 95th percentile is a convention. Reporting the mask at two or three thresholds is more honest than picking one.

Why tree models make this more important

A random forest predicts the mean of training targets in a leaf, so it cannot produce a value outside the training target range. Beyond the training feature space it reuses its most extreme splits and saturates.

The result is a smooth, in-range, plausible map that contains no information. It does not look wrong, which is exactly the problem.

Measured, the same model transferred to the eastern half with coordinates included scored −0.081 — close to predicting the mean — because it had been flattened towards a constant. Without an applicability mask, that map would be indistinguishable from a good one.

An applicability threshold derived from within-training neighbour distances rather than from an absolute value. the threshold is self-calibratingscale the featuresor the widest dominatesweight byimportanceignored features shouldnot countmeasurewithin-training…k nearest, per rowthreshold at the95th percentiledense sample → tightenvelopeAn absolute distance means nothing — it depends on the features and their number spatialworkflow.io
Comparing against the training set itself is what makes the envelope adapt to how densely it was sampled.

Weights, univariate ranges and thresholds

  • Scale the features before measuring distance, or the largest-range feature dominates.
  • Weight by importance; features the model ignores should not drive the mask.
  • Check univariate ranges too — multivariate distance can hide one extreme feature.
  • Derive the threshold from within-training distances, not from an absolute number.
  • Report the mask at several thresholds rather than one.
  • The mask is not an error estimate. It says where the error is unknown.
  • Tree models saturate outside their range, producing plausible flat maps.
  • Write the interpretation into the file tags, not just the documentation.

FAQ

What is the area of applicability?

The part of a prediction map where the covariates resemble the training data closely enough that the model is interpolating rather than guessing.

How is it different from a buffer around the training points?

It is measured in feature space, not geographic space. A distant location with familiar terrain is inside; an adjacent location with an unfamiliar combination is outside.

Can I estimate the error outside the area?

No. The training data contains no information about it. The mask says where the error is unknown, which is the honest available statement.

Why not just use cross-validation?

Because it resamples the training distribution. Block CV promised +0.033 where a real transfer test gave −0.409.

What threshold should I use?

The 95th percentile of within-training neighbour distances is a reasonable convention. Reporting two or three thresholds is more informative than one.

Should I weight by feature importance?

Yes, as a matter of principle — otherwise a feature the model barely uses contributes as much to the distance as the one it depends on. Here it moved the share inside by about one point.

What do I ship with the prediction?

The prediction, the mask, and the feature-space distance, with the interpretation written into the file's tags.