Spatial Machine Learning

A tree ensemble is the default predictive tool almost everywhere else, and on spatial data it fails in two specific, repeatable ways. The first is structural: a random forest predicts a leaf mean, so its surface is piecewise-constant and cannot leave the range of the observed response no matter what the covariates do. The second is procedural and does more damage: these models are still routinely validated with random k-fold, which on autocorrelated data measures interpolation between near-duplicates and reports an accuracy that evaporates the moment the model is asked to predict somewhere nobody sampled. This page, part of Python Workflows for Spatial Modeling & Regression, names both failures, fixes them with feature engineering that supplies spatial information honestly, spatial cross-validation, permutation importance computed under those folds, and an area-of-applicability mask that says where the map should not be trusted — and it is explicit about the cases where kriging remains the better instrument.

Prerequisites #

• Python 3.10+
• scikit-learn>=1.4, geopandas>=1.0, libpysal>=4.9, scikit-gstat>=1.0, numpy>=1.24, pandas>=2.0, scipy>=1.11
• A projected (metric) CRS on both the samples and the prediction grid, so block sizes and neighbourhood radii are metres
• Point observations with a continuous response and a covariate stack sampled at those points, plus the same stack rasterised over the prediction area
• A working knowledge of resampling design — see Cross-Validation Strategies before choosing a fold scheme

Mathematical Core #

A forest prediction is a convex combination #

Let $\hat{f}$ be a random forest of $B$ trees fitted to training pairs $(\mathbf{x}_i, y_i)$, $i = 1, \dots, n$. Write $L_b(\mathbf{x})$ for the set of training indices that fall in the same leaf as $\mathbf{x}$ in tree $b$. The forest prediction is the average over trees of the leaf means:

$$\hat{f}(\mathbf{x}) = \frac{1}{B}\sum_{b=1}^{B} \frac{1}{|L_b(\mathbf{x})|} \sum_{i \in L_b(\mathbf{x})} y_i .$$

Collecting the coefficient on each $y_i$ gives the more revealing form

$$\hat{f}(\mathbf{x}) = \sum_{i=1}^{n} \alpha_i(\mathbf{x})\, y_i , \qquad \alpha_i(\mathbf{x}) = \frac{1}{B}\sum_{b=1}^{B} \frac{\mathbb{1}\{ i \in L_b(\mathbf{x}) \}}{|L_b(\mathbf{x})|} ,$$

where $\mathbb{1}\{\cdot\}$ is the indicator function. Every $\alpha_i(\mathbf{x}) \ge 0$ and $\sum_i \alpha_i(\mathbf{x}) = 1$, because each tree contributes weights summing to one. The prediction is therefore a convex combination of the observed responses, and a convex combination satisfies

$$\min_i y_i \;\le\; \hat{f}(\mathbf{x}) \;\le\; \max_i y_i \quad \text{for every } \mathbf{x}.$$

This is not a tuning problem. No depth, no number of trees, no feature set changes it. The weights $\alpha_i(\mathbf{x})$ are what Lin and Jeon called adaptive nearest-neighbour weights: the forest is a data-driven kernel smoother, and outside the range of the covariates it saw, the kernel simply stops moving. Gradient boosting inherits a weakened version of the same limitation, because each additive stage is itself a tree fitted to residuals within the sampled range.

What random k-fold actually measures #

Take $n$ observations in a study area of size $A$ and hold out a random tenth. Under complete spatial randomness the training points have intensity $\lambda = 0.9n/A$, and the mean distance from a held-out point to its nearest training point is

$$\mathbb{E}[d_{\text{NN}}] = \frac{1}{2\sqrt{\lambda}} .$$

Under blocking with square blocks of side $L$, a held-out point in the interior of a removed block has no training point nearer than the block edge. For a point uniform in the block, $\Pr(D > t) = (1 - 2t/L)^2$ for $0 \le t \le L/2$, so

$$\mathbb{E}[D] = \int_0^{L/2}\left(1 - \frac{2t}{L}\right)^{2} dt = \frac{L}{6},$$

and that is a floor, not the value: the nearest training point lies some distance beyond the edge. The two numbers are the whole argument. A validation scheme interrogates the model at a characteristic prediction distance, and reports an error appropriate to that distance. If the distance implied by your folds is far smaller than the distance implied by your prediction grid, the reported error is not the error you will get.

Autocorrelation is what converts this geometry into an inflated score. If the response has semivariance $\gamma(h)$, the irreducible error of predicting at lag $h$ from a single neighbour is of order $\gamma(h)$. Under random folds $h$ is tiny, $\gamma(h)$ is close to the nugget, and the model looks excellent. Under blocked folds $h$ approaches or exceeds the range, $\gamma(h)$ approaches the sill, and the model reports what it will actually deliver. Both numbers are correct answers to different questions; only one of them is your question.

Making spatial information available honestly #

Trees split on one feature at a time, so raw coordinates $(x, y)$ buy only axis-aligned partitions of the map. Rotating them supplies oblique boundaries at negligible cost:

$$x_i^{(\theta)} = x_i \cos\theta + y_i \sin\theta , \qquad \theta \in \{30^\circ, 60^\circ, 120^\circ\}.$$

A more expressive alternative is a Euclidean distance field to a small fixed set of anchors $\mathbf{a}_1, \dots, \mathbf{a}_K$ — typically the corners and centroid of the study frame:

$$d_{ik} = \lVert \mathbf{s}_i - \mathbf{a}_k \rVert_2 , \qquad k = 1, \dots, K.$$

The published buffer-distance approach sets $K = n$, using the distance to every training point. It works, and it is exactly the construction that makes random k-fold look magnificent, because a held-out point’s own buffer-distance column is a near-perfect location fingerprint of its neighbours. It is also $O(n)$ features and $O(nG)$ memory over a $G$-cell grid, which is untenable at map scale. A handful of anchors captures the same broad positional information for a fixed cost.

Neighbourhood aggregates are the most defensible spatial features, because they encode context rather than identity. Given a row-standardised spatial weight matrix $W$ with entries $w_{ij}$, the spatial lag of a predictor $x$ is

$$(Wx)_i = \sum_{j} w_{ij} x_j , \qquad \sum_j w_{ij} = 1 .$$

Lag the predictors, never the response. The spatial lag of $y$ is available at training points and unavailable at prediction cells, and including it produces a model that scores beautifully under any cross-validation scheme and cannot be deployed at all.

Distance-to-feature covariates sit in a different category again, and they are the ones worth the most effort. The distance from a location to the nearest watercourse, road, coastline, mine adit or urban edge is a genuine physical driver, computable everywhere on the grid, and unrelated to where anyone happened to sample. Unlike coordinates it encodes a mechanism, so it survives the move to spatial folds. Build these with a spatial index rather than a nested loop — a cKDTree over densified vertices of the linear features, or geopandas.sjoin_nearest when the geometry is complex — and compute them once for both the samples and the prediction grid from the same source layer, in the same projected CRS. A distance field derived from a road network clipped to the survey extent but from an unclipped network for the grid is a silent, expensive bug: the two feature matrices then mean different things in the same column.

Annotated Implementation #

The worked example is a soil survey: 1,200 topsoil samples of organic carbon collected in 40 field clusters across a 40 km square, with a covariate stack of elevation, slope, mean NDVI, clay percentage and distance to the nearest watercourse. Clustered sampling is the norm in soil, ecology and air quality work, and it is precisely the design under which random k-fold is most misleading.

Build the sample frame and measure the sampling geometry #

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from scipy.spatial import cKDTree

rng = np.random.default_rng(11)
SIDE = 40_000.0            # 40 km study square, metres
N_CLUSTER, PER_CLUSTER = 40, 30

# Field clusters: 40 survey sites, 30 samples each within a 400 m radius.
centres = rng.uniform(1_500, SIDE - 1_500, size=(N_CLUSTER, 2))
theta = rng.uniform(0, 2 * np.pi, size=(N_CLUSTER, PER_CLUSTER))
radius = 400.0 * np.sqrt(rng.uniform(0, 1, size=(N_CLUSTER, PER_CLUSTER)))
px = (centres[:, 0:1] + radius * np.cos(theta)).ravel()
py = (centres[:, 1:2] + radius * np.sin(theta)).ravel()
cluster_id = np.repeat(np.arange(N_CLUSTER), PER_CLUSTER)

# Covariates: a smooth elevation field, its gradient, greenness and texture.
elev = 180 + 140 * np.sin(px / 9_000) + 90 * np.cos(py / 11_000) + rng.normal(0, 12, px.size)
slope = np.abs(2.2 + 0.012 * (elev - 200) + rng.normal(0, 1.1, px.size))
ndvi = np.clip(0.34 + 0.0009 * (elev - 200) - 0.02 * slope + rng.normal(0, 0.05, px.size), 0, 1)
clay = np.clip(22 + 0.03 * (elev - 200) + rng.normal(0, 4.5, px.size), 2, 60)
d_river = np.abs(py - 0.35 * px - 6_000) / np.hypot(1, 0.35)

# Response: a genuine covariate signal, a smooth spatial term the covariates
# do NOT explain, and white noise. The middle term is what makes honest
# validation hard and random k-fold dishonest.
smooth = 6.0 * np.sin(px / 5_200) * np.cos(py / 4_400)
soc = (4.0 + 26.0 * ndvi + 0.09 * clay - 0.0012 * d_river
       + smooth + rng.normal(0, 3.0, px.size))

samples = gpd.GeoDataFrame(
    {"soc": soc, "elevation": elev, "slope": slope, "ndvi_mean": ndvi,
     "clay_pct": clay, "dist_to_river": d_river, "cluster": cluster_id},
    geometry=[Point(a, b) for a, b in zip(px, py)],
    crs="EPSG:32631",                      # projected — metres
)

# The two distances that decide which validation scheme is honest.
tree = cKDTree(np.c_[px, py])
nn = tree.query(np.c_[px, py], k=2)[0][:, 1]

gx, gy = np.meshgrid(np.arange(25, SIDE, 50), np.arange(25, SIDE, 50))
grid_nn = tree.query(np.c_[gx.ravel(), gy.ravel()], k=1)[0]

print(f"samples: {len(samples)}   clusters: {N_CLUSTER}   CRS: {samples.crs}")
print(f"SOC (g/kg): mean {soc.mean():.2f}  sd {soc.std(ddof=1):.2f}  "
      f"range {soc.min():.2f} - {soc.max():.2f}")
print(f"median sample-to-sample distance : {np.median(nn):>7.0f} m")
print(f"median grid-cell-to-sample distance: {np.median(grid_nn):>6.0f} m")

samples: 1200   clusters: 40   CRS: EPSG:32631
SOC (g/kg): mean 24.13  sd 8.61  range 5.02 - 52.87
median sample-to-sample distance :      61 m
median grid-cell-to-sample distance:   1874 m

Those last two lines settle the validation argument before a model has been fitted. Samples sit 61 m apart; the cells you intend to predict sit 1,874 m from the nearest observation. Any scheme that tests the model at 61 m is answering a question nobody asked.

Add the spatial features #

from libpysal.weights import KNN, lag_spatial

coords = np.c_[samples.geometry.x.values, samples.geometry.y.values]

# Raw coordinates: axis-aligned splits only.
samples["coord_x"], samples["coord_y"] = coords[:, 0], coords[:, 1]

# Euclidean distance fields to four fixed anchors (the frame corners).
# Fixed anchors, not per-sample buffer distances: O(K) features, and the
# same four columns can be computed for the prediction grid.
ANCHORS = np.array([[0.0, 0.0], [SIDE, 0.0], [0.0, SIDE], [SIDE, SIDE]])
for k, a in enumerate(ANCHORS):
    samples[f"edf_anchor_{k}"] = np.linalg.norm(coords - a, axis=1)

# Neighbourhood aggregates of the PREDICTORS from a k-nearest weights
# matrix. Row-standardised, so each lag is a local mean, not a local sum.
w = KNN.from_dataframe(samples, k=8)
w.transform = "R"
for col in ("ndvi_mean", "elevation"):
    samples[f"lag_{col}_k8"] = lag_spatial(w, samples[col].values)

FEATURES = ["elevation", "slope", "ndvi_mean", "clay_pct", "dist_to_river",
            "lag_ndvi_mean_k8", "lag_elevation_k8",
            "coord_x", "coord_y",
            "edf_anchor_0", "edf_anchor_1", "edf_anchor_2", "edf_anchor_3"]

X = samples[FEATURES].to_numpy(dtype=np.float32)
y = samples["soc"].to_numpy(dtype=np.float64)
print(f"feature matrix: {X.shape[0]} rows x {X.shape[1]} columns")

feature matrix: 1200 rows x 13 columns

The lag features are computed once, on the full sample set. That is a mild optimism — inside a fold, a training point’s lag may involve a neighbour that is currently held out — and it is worth knowing about, but it is second-order compared with the fold scheme itself. If you need it eliminated, rebuild w inside each fold from the training rows only and recompute the lags there; on 1,200 points that costs a few milliseconds per fold.

Choose the block size from the residual variogram #

Block width is the one parameter of spatial cross-validation that matters, and it should not be guessed. Fit a quick non-spatial model, take its residuals, and estimate the range at which their autocorrelation dies out. Blocks narrower than that range still leak.

import skgstat as skg
from sklearn.ensemble import RandomForestRegressor

# Covariate-only model: no coordinates, no distance fields. Its residuals
# carry the spatial structure the covariates fail to explain, which is
# exactly the structure that leaks across fold boundaries.
plain = ["elevation", "slope", "ndvi_mean", "clay_pct", "dist_to_river"]
rf0 = RandomForestRegressor(n_estimators=300, min_samples_leaf=5,
                            random_state=0, n_jobs=-1)
rf0.fit(samples[plain], y)
resid = y - rf0.predict(samples[plain])

V = skg.Variogram(coords, resid, model="spherical", n_lags=15,
                  maxlag=15_000, normalize=False)
rng_m, psill, nugget = V.parameters
print(f"residual variogram (spherical): range {rng_m:.0f} m  "
      f"partial sill {psill:.1f}  nugget {nugget:.1f}")
print(f"nugget-to-total ratio: {nugget / (psill + nugget):.2f}")

BLOCK = 5_000.0
print(f"recommended block side >= {rng_m:.0f} m  ->  using {BLOCK:.0f} m")

residual variogram (spherical): range 4183 m  partial sill 41.2  nugget 12.6
nugget-to-total ratio: 0.23
recommended block side >= 4183 m  ->  using 5000 m

A nugget-to-total ratio of 0.23 says that three quarters of the residual variance is spatially structured — the covariates are not carrying the signal on their own, and the model has a great deal of autocorrelation available to memorise. Blocks of 5 km comfortably exceed the 4,183 m range.

Run both validation schemes and compare #

from sklearn.model_selection import KFold, GroupKFold, cross_val_predict
from sklearn.metrics import r2_score, root_mean_squared_error

rf = RandomForestRegressor(n_estimators=500, min_samples_leaf=5,
                           max_features=0.4, random_state=0, n_jobs=-1)

# Block label: integer cell index of a 5 km grid, used as a group.
bx = np.floor(coords[:, 0] / BLOCK).astype(int)
by = np.floor(coords[:, 1] / BLOCK).astype(int)
blocks = bx * 100 + by

pred_random = cross_val_predict(rf, X, y, cv=KFold(10, shuffle=True,
                                                   random_state=0))
pred_block = cross_val_predict(rf, X, y, cv=GroupKFold(n_splits=10),
                               groups=blocks)

def fold_distance(cv_iter):
    """Median distance from each held-out point to its nearest training point."""
    out = []
    for tr, te in cv_iter:
        d = cKDTree(coords[tr]).query(coords[te], k=1)[0]
        out.append(d)
    return np.median(np.concatenate(out))

d_rand = fold_distance(KFold(10, shuffle=True, random_state=0).split(X))
d_block = fold_distance(GroupKFold(n_splits=10).split(X, y, groups=blocks))

for name, p, d in (("random 10-fold", pred_random, d_rand),
                   ("blocked 5 km", pred_block, d_block)):
    print(f"{name:<16} R2 = {r2_score(y, p):.3f}   "
          f"RMSE = {root_mean_squared_error(y, p):.2f} g/kg   "
          f"median prediction distance = {d:.0f} m")

random 10-fold   R2 = 0.870   RMSE = 3.10 g/kg   median prediction distance = 68 m
blocked 5 km     R2 = 0.340   RMSE = 6.99 g/kg   median prediction distance = 3163 m

The same model, the same features, the same data. Reported skill falls from an $R^2$ of 0.870 to 0.340 and error more than doubles, because the second number is measured at a prediction distance two orders of magnitude larger. Neither figure is a bug; the first simply describes a task — filling 68 m gaps inside a survey cluster — that nobody is asking the model to perform.

When the design is clustered, grouping on the cluster identifier is often more faithful than a regular grid, because the cluster is the unit of dependence: GroupKFold(n_splits=10) with groups=samples["cluster"] leaves whole survey sites out and reproduces the deployment task directly. The mechanics of both variants, including how to keep fold sizes balanced when blocks hold very different sample counts, are worked through in Spatial Block Cross-Validation in Python.

Rank predictors under the spatial folds #

Permutation importance for feature $j$ is the average degradation in a loss $\mathcal{L}$ when that column is shuffled:

$$\mathrm{PI}_j = \frac{1}{R}\sum_{r=1}^{R}\Big[\, \mathcal{L}\big(y_{\text{test}}, \hat{f}(X^{(j,r)}_{\text{test}})\big) - \mathcal{L}\big(y_{\text{test}}, \hat{f}(X_{\text{test}})\big) \Big],$$

where $X^{(j,r)}$ is the test matrix with column $j$ permuted at repeat $r$. The quantity is only as meaningful as the test set it is measured on. Computed on randomly held-out rows it rewards features that identify where a point is; computed on spatially held-out blocks it rewards features that transfer.

from sklearn.inspection import permutation_importance

def importance_under(cv_iter, n_repeats=10):
    acc = np.zeros(len(FEATURES))
    n_folds = 0
    for tr, te in cv_iter:
        m = RandomForestRegressor(n_estimators=500, min_samples_leaf=5,
                                  max_features=0.4, random_state=0, n_jobs=-1)
        m.fit(X[tr], y[tr])
        # neg_root_mean_squared_error: importances_mean is the INCREASE in
        # RMSE (g/kg) caused by permuting the column.
        r = permutation_importance(m, X[te], y[te], n_repeats=n_repeats,
                                   scoring="neg_root_mean_squared_error",
                                   random_state=0, n_jobs=-1)
        acc += r.importances_mean
        n_folds += 1
    return acc / n_folds

imp_rand = importance_under(KFold(10, shuffle=True, random_state=0).split(X))
imp_block = importance_under(GroupKFold(n_splits=10).split(X, y, groups=blocks))

table = (pd.DataFrame({"feature": FEATURES,
                       "random_10fold": imp_rand,
                       "blocked_5km": imp_block})
           .sort_values("random_10fold", ascending=False)
           .head(10))
print(table.to_string(index=False,
                      float_format=lambda v: f"{v:.2f}"))

         feature  random_10fold  blocked_5km
    edf_anchor_3           2.41         0.12
         coord_y           1.98         0.09
         coord_x           1.77         0.11
    edf_anchor_1           1.52         0.07
       ndvi_mean           0.86         1.94
       elevation           0.61         1.55
lag_ndvi_mean_k8           0.44         0.87
   dist_to_river           0.39         0.72
        clay_pct           0.18         0.44
           slope           0.12         0.31

The ordering inverts completely. Under random folds the four geometry features are the top four, together worth 7.68 g/kg of RMSE; under blocked folds they are worth 0.39 g/kg in total and the covariates take over. That is the memorisation being made visible: coordinates and distance fields are useful for locating a held-out point among its own cluster-mates, and useless for predicting ground the model has never seen. Keep them in the model — they still absorb residual trend — but read their importance as a diagnostic, not as science. A model whose geometry features remain dominant under spatial folds has covariates that are not doing their job.

Configuring the Area of Applicability #

Cross-validation gives one error figure for the whole map. It does not say that the error is 5 g/kg in the surveyed valleys and 10 g/kg on the ridge nobody visited. The area of applicability, in the form proposed by Meyer and Pebesma, answers that by asking a feature-space question rather than a geographic one: is this prediction cell similar to something in the training set, in the predictors that actually matter?

Standardise each predictor to zero mean and unit variance using the training statistics, weight dimension $j$ by its normalised permutation importance $w_j$ with $\sum_j w_j = 1$, and compute for a new point $\mathbf{x}^*$ the minimum weighted distance to any training row:

$$d(\mathbf{x}^*) = \min_{i \in \text{train}} \Bigl[ \sum_{j=1}^{p} w_j \bigl( (x^*_j - x_{ij}) / \sigma_j \bigr)^{2} \Bigr]^{1/2} .$$

Divide by $\bar{d}$, the mean of all pairwise weighted distances among training rows, to obtain a scale-free dissimilarity index:

$$DI(\mathbf{x}^*) = \frac{d(\mathbf{x}^*)}{\bar{d}} .$$

The threshold comes from the training data itself. For each training row, compute its $DI$ against training rows from other folds — the same separation the cross-validation used — and set

$$t = Q_{0.75}\big(DI_{\text{train}}\big) + 1.5 \cdot \mathrm{IQR}\big(DI_{\text{train}}\big),$$

the standard outlier fence. The area of applicability is $\{\mathbf{x}^* : DI(\mathbf{x}^*) \le t\}$.

from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import NearestNeighbors
from scipy.spatial.distance import pdist

wj = np.clip(imp_block, 0, None)
wj = wj / wj.sum()                     # importance weights, sum to 1
scaler = StandardScaler().fit(X)
Zt = scaler.transform(X) * np.sqrt(wj)  # sqrt(w) so Euclidean d^2 = sum w_j d_j^2

dbar = pdist(Zt).mean()

# Threshold: each training row's distance to the nearest training row in a
# DIFFERENT fold, so the fence reflects the separation CV actually used.
di_train = np.empty(len(Zt))
for tr, te in GroupKFold(n_splits=10).split(X, y, groups=blocks):
    nn_tr = NearestNeighbors(n_neighbors=1).fit(Zt[tr])
    di_train[te] = nn_tr.kneighbors(Zt[te])[0].ravel() / dbar

q1, q3 = np.quantile(di_train, [0.25, 0.75])
thresh = q3 + 1.5 * (q3 - q1)

# Apply to the 800 x 800 prediction grid, chunked to bound peak memory.
nn_all = NearestNeighbors(n_neighbors=1).fit(Zt)
Zg = scaler.transform(grid_features) * np.sqrt(wj)   # grid_features: (640000, 13)
di_grid = np.concatenate([
    nn_all.kneighbors(Zg[s:s + 50_000])[0].ravel() / dbar
    for s in range(0, len(Zg), 50_000)
])
inside = di_grid <= thresh

print(f"mean pairwise training distance : {dbar:.3f}")
print(f"DI threshold (Q75 + 1.5 IQR)    : {thresh:.4f}")
print(f"grid cells inside AOA           : {inside.sum():,} of {inside.size:,} "
      f"({100 * inside.mean():.1f}%)")

# di_train was already computed fold-wise above, so it doubles as the
# per-observation DI under exactly the folds pred_block came from.
m_in = di_train <= thresh
print(f"blocked-CV RMSE inside AOA      : "
      f"{root_mean_squared_error(y[m_in], pred_block[m_in]):.2f} g/kg")
print(f"blocked-CV RMSE outside AOA     : "
      f"{root_mean_squared_error(y[~m_in], pred_block[~m_in]):.2f} g/kg")

mean pairwise training distance : 3.214
DI threshold (Q75 + 1.5 IQR)    : 0.4413
grid cells inside AOA           : 440,960 of 640,000 (68.9%)
blocked-CV RMSE inside AOA      : 5.42 g/kg
blocked-CV RMSE outside AOA     : 9.60 g/kg

The mask is doing real work: error inside the area of applicability is 5.42 g/kg, outside it is 9.60 g/kg, and the two combine to the 6.99 g/kg headline figure. Nearly a third of the map is a region where the model is being asked about covariate combinations it never saw, and where a piecewise-constant predictor will simply return the nearest leaf mean with no signal that it is guessing.

Output Interpretation #

Read four things off a spatial machine learning run, in this order.

The ratio between random and spatial scores. A gap of a few points is a well-specified model with informative covariates. A gap like 0.870 against 0.340 says the model was living on autocorrelation, and the honest number is the smaller one. If the two scores are almost equal, check that your blocks are actually wider than the residual range before congratulating yourself — equal scores are more often a sign of blocks that leak than of a model that transfers.

Where the geometry features sit in the importance ranking. Under spatial folds, coordinates and distance fields should fall to the bottom. If they stay at the top, the covariates are not explaining the response and you are fitting a very expensive, very badly regularised interpolator. That is the point at which to consider kriging, or a combination.

The shape of the predicted-versus-observed scatter at the extremes. Regression to the mean at both tails is the convex-combination bound made visible: high observations are systematically under-predicted, low ones over-predicted, and the slope of the fit against observed values is well below one. A slope of 0.6–0.8 under spatial folds is typical and tolerable; a slope near 0.4 means the model is mostly returning the global mean and the map is decoration.

The fraction of the grid outside the area of applicability. Below roughly 10% is comfortable. Above 30%, as here, the sampling design did not cover the covariate space and no amount of modelling will fix it — the answer is more samples in the unrepresented conditions, chosen by sampling the high-$DI$ region rather than by convenience. A figure near zero is worth checking rather than celebrating: it usually means the importance weights collapsed onto one feature, making the distance nearly one-dimensional.

When Kriging Is Still the Better Tool #

Being even-handed here matters, because the tree ensemble is not always the right answer and often is not the best one on its own. Fitting all three approaches under the same blocked folds settles it empirically:

from pykrige.ok import OrdinaryKriging

# Ordinary kriging under the same GroupKFold splits, for a like-for-like
# comparison: no covariates at all, only the coordinates and the variogram.
ok_pred = np.empty(len(y))
rk_pred = np.empty(len(y))
for tr, te in GroupKFold(n_splits=10).split(X, y, groups=blocks):
    ok = OrdinaryKriging(coords[tr, 0], coords[tr, 1], y[tr],
                         variogram_model="spherical",
                         variogram_parameters={"sill": 53.8, "range": 4183.0,
                                               "nugget": 12.6},
                         enable_plotting=False)
    ok_pred[te], _ = ok.execute("points", coords[te, 0], coords[te, 1])

    # Regression kriging: the same forest, with its training residuals
    # kriged and added back at the held-out locations.
    m = RandomForestRegressor(n_estimators=500, min_samples_leaf=5,
                              max_features=0.4, random_state=0, n_jobs=-1)
    m.fit(X[tr], y[tr])
    rk = OrdinaryKriging(coords[tr, 0], coords[tr, 1], y[tr] - m.predict(X[tr]),
                         variogram_model="spherical",
                         variogram_parameters={"sill": 53.8, "range": 4183.0,
                                               "nugget": 12.6},
                         enable_plotting=False)
    resid_te, _ = rk.execute("points", coords[te, 0], coords[te, 1])
    rk_pred[te] = m.predict(X[te]) + resid_te

for name, p in (("ordinary kriging", ok_pred),
                ("random forest + spatial features", pred_block),
                ("regression kriging (RF + OK resid)", rk_pred)):
    print(f"{name:<34} RMSE = {root_mean_squared_error(y, p):.2f}  "
          f"R2 = {r2_score(y, p):.3f}")

ordinary kriging                   RMSE = 7.62  R2 = 0.217
random forest + spatial features   RMSE = 6.99  R2 = 0.340
regression kriging (RF + OK resid) RMSE = 6.31  R2 = 0.463

Three lessons sit in those numbers. The forest beats plain kriging here because the covariates carry genuine signal, and kriging alone ignores them. But the combination beats both, because the forest models the covariate response while the kriged residuals recover the smooth spatial term the covariates cannot express — the smooth component built into the synthetic data. That hybrid is Regression Kriging, and on covariate-rich spatial problems it is usually the strongest single choice.

Reach for kriging rather than an ensemble when sampling is dense relative to the variogram range, so that interpolation genuinely is the task; when covariates are weak, missing, or available only at coarse support; when the target is smooth and predictions outside the observed range are physically meaningful, since kriging’s weights are not constrained to be non-negative and it will extrapolate a trend; and above all when you need a calibrated prediction variance. Kriging returns a variance surface from the same model that produced the estimate. A forest returns a spread of leaf means, which is a measure of tree disagreement, not of predictive uncertainty, and is systematically too narrow in exactly the extrapolation regions where you most need it to be wide. The full side-by-side, including quantile regression forests as a partial answer to the uncertainty problem, is in Spatial Machine Learning vs Kriging.

Reach for the ensemble when you have a deep covariate stack with non-linear and interacting effects, when the relationship between covariates and response is non-stationary across the study area, when the response is bounded and the convex-combination property is a feature rather than a bug, and when the covariate rasters are exhaustive so that prediction is cheap everywhere.

Production Considerations #

Training cost. Fitting a forest is $O(B \cdot m \cdot n \log n)$ where $B$ is the number of trees, $m$ the features considered per split (max_features), and $n$ the training rows. Ten-fold spatial cross-validation multiplies this by $k$, and permutation importance with $R$ repeats over $p$ features adds $k \cdot p \cdot R$ full prediction passes on top. With $p = 13$, $R = 10$ and $k = 10$ that is 1,300 extra prediction passes — usually the dominant cost of the whole workflow. Reduce $R$ to 5 for routine runs and reserve 10 or 20 for the reported figure.

Prediction cost and grid size. Predicting one cell is $O(B \log n)$; a 640,000-cell grid with 500 trees is a few seconds with n_jobs=-1. The constraint is memory in the feature stack, not compute: 640,000 cells by 13 features in float64 is 66 MB, entirely manageable. Buffer-distance features are what break this. With $K = n = 1{,}200$ distance columns the same grid needs $640{,}000 \times 1{,}200 \times 8$ bytes, or 6.1 GB in float64 and 3.1 GB in float32. Four fixed anchors need 20 MB. That factor of 300 is the practical argument against per-sample buffer distances at map scale, quite apart from the validation problem they create.

Area-of-applicability cost. The nearest-neighbour search is $O(G \cdot n \cdot p)$ in the worst case for $G$ grid cells. KD-trees degrade towards brute force above roughly ten dimensions, so assume brute force and budget accordingly: chunk the grid as in the code above so peak memory is bounded by the chunk size rather than $G$, and cache Zt, dbar and thresh alongside the fitted model, because they are part of the model artefact. A model shipped without its area of applicability cannot be used responsibly.

Caching and reproducibility. Sample the covariate stack at the observation points once and persist the result; re-extracting from rasters on every experiment is the slowest step in most real pipelines and the easiest to get subtly wrong when CRS or resampling settings drift. Store the grid feature stack as a memory-mapped float32 array so that repeated prediction runs do not re-read it. Persist the block labels with the sample table: a fold scheme that changes between the tuning run and the reported run silently invalidates both.

Parallelisation. Use n_jobs=-1 on the forest and leave the outer cross-validation loop serial; nesting joblib pools oversubscribes cores and often runs slower than the serial version. For very large surveys, dask-geopandas handles the covariate extraction and the block assignment, but the forest itself should be fitted on a single machine unless $n$ exceeds a few million rows, at which point subsampling within blocks costs less accuracy than distributed training costs in complexity.

Troubleshooting #

Next Steps #

Two child guides carry this material further. Random Forest with Spatial Features and Spatial CV is the end-to-end build: feature construction, hyperparameter selection under blocked folds, and the reporting template that puts both scores side by side. Spatial Machine Learning vs Kriging runs the comparison properly across sampling densities and covariate strengths, and covers prediction uncertainty from both sides. For the fold machinery itself, work through Cross-Validation Strategies and its practical companion Spatial Block Cross-Validation in Python.

Frequently Asked Questions #

Why can a random forest never extrapolate beyond the observed range?

Every prediction is an average of leaf means, and every leaf mean is an average of training responses. Written out, the prediction is a weighted sum of the observed values with non-negative weights that sum to one — a convex combination. A convex combination of numbers cannot be smaller than the smallest or larger than the largest, so the fitted surface is bounded by the observed range whatever the covariates say. Gradient boosting escapes this only partially, because its additive updates are still fitted to residuals inside the sampled range.

Is random k-fold ever acceptable on spatial data?

Only when the model will be deployed at the same spatial density it was trained at, and the observations are not spatially clustered. If you are filling gaps between sensors that will remain in place, random k-fold answers the right question. If you are producing a map over ground where no sample was taken, it does not: it measures interpolation between near-duplicates. The test is to compare the median distance from a held-out point to its nearest training point against the same distance for the cells you intend to predict.

Should I put raw coordinates into the feature matrix?

Yes, but only alongside spatial cross-validation, and expect them to be ranked low once validation is honest. Raw coordinates give the trees axis-aligned splits, which carve the study area into rectangles and let the model memorise where samples are rather than what conditions produce the response. Rotated coordinates or distance fields to a few fixed anchors give oblique boundaries at the same cost. If coordinates dominate importance under spatial folds, the covariates are not carrying the signal.

When is kriging still the better tool than a tree ensemble?

When sampling is dense relative to the autocorrelation range, when covariates are weak or unavailable, when the target is smooth and you need to extrapolate past the observed range, and whenever a calibrated prediction variance matters more than a point estimate. Kriging gives an interpolator that honours the data and a variance surface derived from the same model. A forest gives neither. The usual best answer is to combine them: fit the trend with the ensemble and krige the residuals.

Related

• Random Forest with Spatial Features and Spatial CV — the end-to-end build with hyperparameter selection under blocked folds
• Spatial Machine Learning vs Kriging — the comparison run properly, across sampling densities and covariate strengths
• Spatial Block Cross-Validation in Python — block construction, fold balancing and buffered variants
• Cross-Validation Strategies — choosing a resampling design that matches the deployment task
• Regression Kriging — fitting the trend with a model and kriging what it leaves behind

← Back to Python Workflows for Spatial Modeling & Regression