Random Forest with Spatial Features and Spatial CV

TL;DR: Fit RandomForestRegressor on coordinates, oblique coordinates, reference distances and libpysal neighbourhood means, then score it twice — KFold(n_splits=5, shuffle=True) and GroupKFold(n_splits=5) with block identifiers as groups. On the worked example the two return R² of 0.851 and 0.548. The gap of 0.303 is the optimism; the second number is what deploys.

Why This Matters #

A random forest handed the coordinates of its training points will always score well under random k-fold, because random k-fold hands it back a held-out point sitting a few hundred metres from a dozen training points. The forest does not need to have learned anything about the process; it only needs to find the leaf containing the neighbours and return their mean. That is interpolation between near-duplicates, and it is not the task you deploy the model for. This page is part of Spatial Machine Learning, and it is the practical companion to Spatial Machine Learning vs Kriging: the argument there is about which family of model to choose, the argument here is about how to find out whether the one you chose is any good.

The whole pipeline is short — perhaps sixty lines — but three of its parts are routinely got wrong. The features are usually just x and y, which forces the tree to approximate every oblique gradient with a staircase. The validation is usually random, which reports a number nobody will ever reproduce on new ground. And the predictions are usually made everywhere, including in feature-space regions the training set never covered, where a forest silently returns the average of whatever leaf it lands in rather than admitting it has nothing to say.

Environment and Version Pinning #

pip install "numpy==1.26.4" "pandas==2.2.2" "scipy==1.13.1" \
            "geopandas==1.0.1" "libpysal==4.12.1" "esda==2.6.0" \
            "scikit-learn==1.5.1"

import numpy as np
import pandas as pd
import geopandas as gpd
import libpysal
from libpysal.weights import lag_spatial
from esda.moran import Moran
from scipy.spatial.distance import pdist
from sklearn.base import clone
from sklearn.ensemble import RandomForestRegressor
from sklearn.inspection import permutation_importance
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import KFold, GroupKFold, cross_val_score, cross_val_predict
from sklearn.neighbors import NearestNeighbors
from sklearn.preprocessing import StandardScaler

Step-by-Step Implementation #

1. Build an autocorrelated sample with a known structure #

The synthetic field below has two environmental covariates and a smooth residual component that no covariate explains. That residual is what makes the validation strategy matter: it is predictable from location within a few kilometres and unpredictable beyond, which is exactly the situation in soil, air-quality and yield mapping.

rng = np.random.default_rng(2026)
N, SIDE = 1200, 20_000.0                       # 1,200 samples in a 20 km square

px = rng.uniform(0, SIDE, N)
py = rng.uniform(0, SIDE, N)


def smooth_field(x, y, wavelengths, seed):
    """Sum of randomly oriented sinusoids: a cheap, reproducible smooth surface."""
    g = np.random.default_rng(seed)
    out = np.zeros_like(x)
    for lam in wavelengths:
        theta = g.uniform(0, np.pi)
        phase = g.uniform(0, 2 * np.pi)
        k = 2 * np.pi / lam
        out += np.sin(k * (np.cos(theta) * x + np.sin(theta) * y) + phase)
    return out / np.sqrt(len(wavelengths))


elev = 180.0 + 45.0 * smooth_field(px, py, (18000., 11000., 7000.), 11)
ndvi = 0.55 + 0.12 * smooth_field(px, py, (9000., 5000., 3000.), 12)
resid_field = smooth_field(px, py, (12000., 6500.), 13)

z = (18.0
     + 0.050 * (elev - 180.0)          # covariate response
     + 10.5 * (ndvi - 0.55)
     + 3.90 * resid_field              # smooth, unexplained by any covariate
     + rng.normal(0, 0.55, N))         # nugget

pts = gpd.GeoDataFrame(
    {"elev": elev, "ndvi": ndvi, "z": z},
    geometry=gpd.points_from_xy(px, py),
    crs="EPSG:32630",                  # projected, metres
)
print(f"n = {N}   z: mean {z.mean():.2f}  sd {z.std(ddof=1):.2f}")

n = 1200   z: mean 18.31  sd 3.34

2. Construct the spatial feature block #

Four families of feature go in. Raw coordinates give the forest absolute position. Oblique coordinates — the projection of each point onto axes rotated by 30°, 60°, 120° and 150° — let a single split cut diagonally, which raw x and y cannot do. Distances to a fixed set of reference points supply a smooth radial basis. Neighbourhood means of the covariates, computed from a k-nearest-neighbour weights object, give the forest local context.

def oblique_coords(x, y, angles_deg=(30, 60, 120, 150)):
    """Project coordinates onto rotated axes so trees can cut obliquely."""
    cols = {}
    for a in angles_deg:
        r = np.deg2rad(a)
        cols[f"oc_{a}"] = np.cos(r) * x + np.sin(r) * y
    return pd.DataFrame(cols)


REF = np.array([[0., 0.], [SIDE, 0.], [0., SIDE],
                [SIDE, SIDE], [SIDE / 2, SIDE / 2]])
xy = np.column_stack([px, py])
dist_ref = pd.DataFrame(
    np.linalg.norm(xy[:, None, :] - REF[None, :, :], axis=-1),
    columns=[f"d_ref{i}" for i in range(len(REF))],
)

w = libpysal.weights.KNN.from_dataframe(pts, k=8)
w.transform = "R"                                  # row-standardised
pts = pts.reindex(w.id_order)                      # align before lagging
lags = pd.DataFrame({
    "elev_lag": lag_spatial(w, pts["elev"].values),
    "ndvi_lag": lag_spatial(w, pts["ndvi"].values),
})

X = pd.concat([
    pd.DataFrame({"x": px, "y": py}),
    oblique_coords(px, py),
    dist_ref,
    pd.DataFrame({"elev": elev, "ndvi": ndvi}),
    lags,
], axis=1)

print(X.shape)
print(list(X.columns))
print(f"Moran's I of z (k=8): {Moran(z, w, permutations=999).I:.3f}")

(1200, 15)
['x', 'y', 'oc_30', 'oc_60', 'oc_120', 'oc_150', 'd_ref0', 'd_ref1',
 'd_ref2', 'd_ref3', 'd_ref4', 'elev', 'ndvi', 'elev_lag', 'ndvi_lag']
Moran's I of z (k=8): 0.713

A Moran’s I of 0.713 confirms the target is strongly autocorrelated, which is the precondition for everything that follows. If your own data returns something near zero, random k-fold and blocked CV will agree and none of the machinery below buys you anything. The weights object is the same construction described in spatial weight matrices, and the reindex(w.id_order) line matters: lag_spatial assumes the value array is in the weights object’s own ordering, not the GeoDataFrame’s.

3. Fit the forest and score it under random k-fold #

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

kf = KFold(n_splits=5, shuffle=True, random_state=0)
r2_rand = cross_val_score(rf, X, z, cv=kf, scoring="r2")
pred_rand = cross_val_predict(rf, X, z, cv=kf)

print("random 5-fold R2  :", np.round(r2_rand, 3), f" mean {r2_rand.mean():.3f}")
print("random 5-fold RMSE:", f"{mean_squared_error(z, pred_rand) ** 0.5:.3f}")

random 5-fold R2  : [0.862 0.849 0.836 0.858 0.85 ]  mean 0.851
random 5-fold RMSE: 1.290

An R² of 0.851 against a target standard deviation of 3.34 looks like a strong model. The five folds agree closely — a standard deviation of 0.009 across them — which reads as stability but is really a symptom: every fold is the same easy problem.

4. Score it again under spatial block cross-validation #

Assign each sample to a 4 km square block, then pass the block identifier as groups to GroupKFold. Every block lands wholly inside one fold, so no held-out point has a training neighbour from its own block.

BLOCK = 4000.0
block_id = (np.floor(px / BLOCK).astype(int) * 100
            + np.floor(py / BLOCK).astype(int))
print("blocks:", len(np.unique(block_id)))

gkf = GroupKFold(n_splits=5)
r2_block = cross_val_score(rf, X, z, cv=gkf, groups=block_id, scoring="r2")
pred_block = cross_val_predict(rf, X, z, cv=gkf, groups=block_id)

print("blocked R2  :", np.round(r2_block, 3), f" mean {r2_block.mean():.3f}")
print("blocked RMSE:", f"{mean_squared_error(z, pred_block) ** 0.5:.3f}")
print("optimism gap:", f"{r2_rand.mean() - r2_block.mean():.3f}")

blocks: 25
blocked R2  : [0.612 0.489 0.571 0.508 0.56 ]  mean 0.548
blocked RMSE: 2.246
optimism gap: 0.303

Nothing about the model changed. The same 500 trees, the same 15 features, the same 1,200 samples. Only the definition of “held out” changed, and R² fell by 0.303 while RMSE rose from 1.29 to 2.25 — from 39% of the target’s standard deviation to 67% of it. Note also that the blocked folds now disagree with each other, with a fold standard deviation of 0.044 against 0.009 for the random folds. That spread is real information about how much your answer depends on which part of the study area you happen to be predicting into. The choice of block size and the alternatives to a square grid are covered in Spatial Block Cross-Validation in Python.

5. Compute permutation importance inside the spatial folds #

Permutation importance measured on a randomly held-out set will over-credit the coordinate features, because shuffling x destroys precisely the near-duplicate lookup the model was relying on. Measure it on blocked held-out sets instead.

def cv_permutation_importance(model, X, y, splitter, groups=None,
                              n_repeats=10, seed=0):
    parts = []
    for tr, te in splitter.split(X, y, groups):
        m = clone(model).fit(X.iloc[tr], y[tr])
        r = permutation_importance(m, X.iloc[te], y[te], scoring="r2",
                                   n_repeats=n_repeats, random_state=seed,
                                   n_jobs=-1)
        parts.append(pd.Series(r.importances_mean, index=X.columns))
    return pd.concat(parts, axis=1).mean(axis=1)


imp = pd.DataFrame({
    "random_cv": cv_permutation_importance(rf, X, z, kf),
    "blocked_cv": cv_permutation_importance(rf, X, z, gkf, groups=block_id),
}).sort_values("random_cv", ascending=False)

print(imp.head(8).round(3).to_string())

          random_cv  blocked_cv
oc_60         0.181       0.048
y             0.166       0.036
x             0.152       0.031
oc_120        0.144       0.027
d_ref4        0.098       0.022
elev          0.041       0.129
elev_lag      0.022       0.071
ndvi          0.017       0.054

The ordering inverts. Under random folds the five most important features are all positional and elev is a distant sixth; under blocked folds elev is first with 0.129, its neighbourhood mean second with 0.071, and every coordinate feature drops below 0.05. Read left to right, the table says: the coordinates are worth 0.181 of R² when the answer is a few hundred metres away, and 0.048 when it is four kilometres away. Report the right-hand column.

6. Mask the area of applicability #

A forest returns a leaf mean for any input you give it, however unlike the training data that input is. The mask below, a simplified form of the Meyer and Pebesma area of applicability, measures each prediction location’s distance to the nearest training point in importance-weighted, standardised feature space, scales it by the mean pairwise training distance, and rejects anything above the 95th percentile of the training set’s own leave-self-out distances.

def aoa(X_train, X_new, weights, q=0.95):
    """Dissimilarity index and threshold in importance-weighted feature space."""
    w_pos = np.clip(np.asarray(weights, float), 0, None)
    w_pos = np.sqrt(w_pos / w_pos.sum())
    sc = StandardScaler().fit(X_train)
    A = sc.transform(X_train) * w_pos
    B = sc.transform(X_new) * w_pos

    scale = pdist(A).mean()                       # mean pairwise training distance
    nn = NearestNeighbors(n_neighbors=2).fit(A)
    d_train = nn.kneighbors(A)[0][:, 1] / scale   # column 0 is the point itself
    thr = np.quantile(d_train, q)
    d_new = nn.kneighbors(B, n_neighbors=1)[0][:, 0] / scale
    return d_new, thr


wts = imp["blocked_cv"].reindex(X.columns).values
inside, thresholds = np.zeros(len(X), dtype=bool), []
for tr, te in gkf.split(X, z, block_id):
    di, thr = aoa(X.iloc[tr], X.iloc[te], wts)
    inside[te] = di <= thr
    thresholds.append(thr)

print(f"mean DI threshold : {np.mean(thresholds):.3f}")
print(f"inside the mask   : {inside.mean():.1%}")
print(f"RMSE inside  : {mean_squared_error(z[inside], pred_block[inside]) ** 0.5:.3f}")
print(f"RMSE outside : {mean_squared_error(z[~inside], pred_block[~inside]) ** 0.5:.3f}")

mean DI threshold : 0.472
inside the mask   : 91.0%
RMSE inside  : 1.980
RMSE outside : 4.052

That last pair is the test the mask has to pass. If the 9% of held-out points it rejects had the same error as the 91% it accepts, the mask would be decoration. Here the rejected points are more than twice as wrong — RMSE 4.05 against 1.98 — and the two combine back to the blocked RMSE of 2.246, so the mask is separating the predictions that can be trusted from the ones that cannot.

Interpreting the Output #

The headline is the pair (0.851, 0.548), not either number alone. The blocked figure is the estimate of skill at a location with no nearby sample, which is what a prediction map is made of. The random figure is the estimate of skill at a location surrounded by samples, which is a situation you have no reason to make a prediction in. The difference between them, 0.303, is a measure of how much of the model’s apparent performance comes from spatial redundancy in the sample rather than from any learned relationship.

What good looks like: a blocked R² comfortably above zero, a gap of under about 0.15, and a blocked permutation importance table led by covariates rather than by coordinates. What that combination says is that the model has found a covariate relationship that transfers, and is using position only to fill in the residual. Warning signs come in three flavours. A gap above 0.3, as here, means the model is largely a spatial interpolator wearing a machine-learning costume — which may be perfectly acceptable, but should be compared against ordinary kriging before you accept the extra complexity. A blocked R² near or below zero means the blocked folds do worse than predicting the global mean, and the model should not be deployed at all. And a blocked fold spread wider than the mean itself means the study area is not homogeneous enough for one model.

The prediction surface is worth looking at directly, because a random forest is piecewise constant by construction. Each tree returns a leaf mean, so the ensemble surface is a mosaic of flat patches with discontinuities at every split value, and the patches are rectangles in feature space. On a map this shows up as visible blockiness, particularly where coordinates dominate the splits. It also has a consequence that no amount of tuning removes: the forest cannot predict outside the range of the training targets. On this dataset the observed z spans [9.81, 27.36] while the blocked out-of-fold predictions span only [12.94, 24.58]. A trend that continues past your sampled range will be flattened, which is the specific failure that Spatial Machine Learning vs Kriging weighs against the forest’s ability to absorb many covariates without a parametric form.

Critical Best Practices #

Never lag the target across the whole sample #

Adding a neighbourhood mean of z to the feature block is the single most common way to destroy a spatial validation. Under any cross-validation scheme the lag of a training point contains the values of its held-out neighbours, so information leaks backwards across the split and both the random and the blocked score become meaningless. Lag the covariates, which are available at every prediction location anyway. If you want the target’s neighbourhood in the model, recompute the lag inside each fold from the training points only, and accept that at prediction time you will have to substitute the model’s own predictions.

Choose the block size from the variogram range, not from convenience #

A 4 km block on this dataset is not arbitrary: it exceeds the range over which the residual field stays correlated, which is what makes the held-out blocks genuinely independent. Blocks smaller than the autocorrelation range leave neighbours on both sides of the split and reproduce the random-k-fold optimism at reduced strength; blocks far larger than the range throw away training data for nothing and inflate the fold-to-fold variance. Estimate an empirical variogram of the target or of a naive model’s residuals first, and size the blocks at or just above its range. Spatial K-Fold Cross-Validation Setup covers the mechanics of getting that number into a splitter.

Align the weights object before computing lags #

libpysal weights carry their own id_order, and lag_spatial assumes the value array you hand it is in that order. When the weights come from a GeoDataFrame with a plain RangeIndex the two usually coincide, so the bug is invisible until someone filters rows, reads a subset, or builds weights from a shapefile with a non-sequential identifier. Then the lag silently attaches each point’s feature to some other point’s neighbours, which is undetectable in the fitted model and shows up only as a mysteriously weak covariate. Call reindex(w.id_order) before lagging, every time.

Compute permutation importance on the same folds you report #

Permutation importance is a property of a model and an evaluation set. Measured on a random held-out set it credits x, y and the oblique coordinates with roughly four times the importance they carry on a blocked set, because shuffling them breaks a lookup that only random splitting made available. Feature selection driven by the random-fold column will therefore keep coordinates and discard covariates — precisely backwards. Wrap permutation_importance in the blocked splitter, average across folds, and use that as the basis for both interpretation and any dimension reduction.

Treat the area-of-applicability mask as a claim to be tested #

The mask is only useful if the predictions it rejects really are worse. Every time you build one, split the blocked out-of-fold predictions by the mask and compare RMSE on each side, as in step six. If the two are within a few per cent of each other the threshold quantile is wrong, the importance weights are dominated by a feature that does not vary spatially, or the study area simply has no extrapolation problem. Reporting a mask that has never rejected anything harmful is worse than reporting none, because it implies a check that was not performed.

Troubleshooting #

Next Steps #

Take the blocked score as the baseline and compare it directly against ordinary kriging on the same folds, following Spatial Machine Learning vs Kriging. If the block geometry itself is what you now want to tune — size, buffering, or leaving a dead zone around each held-out block — that is the subject of Spatial Block Cross-Validation in Python.

Frequently Asked Questions #

Why is the blocked R-squared the number I should report?

Random k-fold puts a held-out point a few hundred metres from several training points, so the forest only has to interpolate between near-duplicates. That is not the task at deployment, where predictions are made for locations with no nearby sample. Blocking removes the near-duplicates from the training side of every split, so the blocked score measures the skill that actually transfers. On the worked example the two disagree by 0.303 in R-squared, and the blocked figure of 0.548 is the one that matched held-out block performance.

Should I feed the spatially lagged target into the features?

No. A neighbourhood mean of the target computed on the full sample carries each held-out value into its neighbours’ features, so cross-validation scores a model that has already seen the answer. Lag the covariates instead, which are available everywhere at prediction time. If you genuinely want the target’s neighbourhood structure in the model, recompute the lag inside each fold using only the training points, or move to a hybrid that models the residual spatially rather than smuggling it in as a feature.

Do rotated coordinates actually help, or are raw x and y enough?

A regression tree splits on one feature at a time, so with only x and y available every decision boundary is a staircase of axis-aligned steps. A trend running diagonally therefore costs many splits and is still approximated badly near the edges. Adding oblique coordinates at 30, 60, 120 and 150 degrees lets a single split cut along those directions. The gain is modest when covariates carry most of the signal, but it removes an artefact that is purely a property of the splitting rule.

How does the area-of-applicability mask differ from a prediction interval?

A prediction interval from quantile regression forests describes spread among the trees, and it stays narrow wherever the trees agree, including in feature-space regions where none of them has seen data. The area-of-applicability mask asks a different question: how far is this location from anything the model was trained on, measured in importance-weighted feature space. It returns a hard boundary rather than a width. Use both. The interval quantifies uncertainty where the model applies; the mask says where it does not apply at all.

Related

• Spatial Machine Learning vs Kriging — when the extra features are worth the loss of a variance model
• Spatial Block Cross-Validation in Python — block size, buffering and the dead zone around each held-out block
• Spatial Weight Matrices — the neighbour structure behind the lagged covariate features

← Back to Spatial Machine Learning