⚙️ How-To Spatial machine learning Raster & rasterio Spatial Machine Learning · step 8 of 14 Updated

How to Predict a Continuous Surface with Random Forest

Problem statement

Predicting a continuous variable across a landscape from point samples and covariate rasters is the workhorse of digital soil mapping, species distribution modelling and environmental prediction. A random forest is the usual first choice, and it has three properties that matter here and are easy to forget.

Measured on a real problem — predicting a vegetation index from terrain covariates, 3,000 training cells:

random 5-fold CV        +0.555
block cross-validation  +0.221
true error over the
whole unsampled map     +0.557

The model works, and which number you quote depends entirely on what it will be asked to do.

Three properties: it cannot predict outside the training target range, it treats coordinates as extremely powerful features, and it gives no calibrated uncertainty without extra work.

Quick answer 

import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GroupKFold
from sklearn.cluster import KMeans
from sklearn.metrics import r2_score

model = RandomForestRegressor(n_estimators=200, min_samples_leaf=2,
                              n_jobs=-1, random_state=0)

blocks = KMeans(10, random_state=0, n_init=10).fit_predict(coords)
predicted = np.empty(len(y))
for train, test in GroupKFold(5).split(X, y, blocks):
    predicted[test] = model.fit(X[train], y[train]).predict(X[test])
print(f"block CV R2 {r2_score(y, predicted):+.3f}")

model.fit(X, y)
surface = model.predict(X_grid)
Sample points, covariate rasters, a fitted forest, and a prediction surface with an applicability mask.
The mask is part of the product, not an optional extra. A forest cannot signal that it is extrapolating.

Step-by-step solution

1. Set min_samples_leaf above 1

The default of 1 lets every leaf hold a single training sample, which memorises the training set. With spatially correlated data that memorisation looks like skill under random cross-validation.

Two to five is a reasonable range for a few thousand samples. It also shrinks the model substantially, which matters when a forest of 200 unpruned trees over 10,000 samples is hundreds of megabytes.

2. Do not tune n_estimators for accuracy

More trees reduce variance and never overfit, so the score rises and then plateaus. Two hundred is enough for most problems; a thousand costs five times the memory and prediction time for a change in the third decimal.

Spend the tuning effort on min_samples_leaf and max_features, which genuinely change the fit.

3. Choose the validation scheme from the deployment

The three numbers at the top are all correct and they answer different questions. Random CV estimates gap-filling inside the sampled area; block CV estimates prediction in unsampled parts of the same region; neither estimates a new region.

4. Remember the forest cannot extrapolate

A forest predicts the mean of training targets in a leaf, so its output is bounded by the training target range. Applied where the covariates are unfamiliar, it reuses its most extreme splits and saturates at a constant.

The map then looks smooth, plausible and in range while carrying no information. That is why the applicability mask matters more for a forest than for a linear model, which at least fails visibly.

5. Get uncertainty from the ensemble spread, carefully 

per_tree = np.stack([tree.predict(X) for tree in model.estimators_])
spread = per_tree.std(axis=0)

The spread across trees measures disagreement, which correlates with uncertainty and is not a calibrated interval. It is useful as a relative map of confidence and should not be quoted as a standard error.

A random forest prediction saturating at the training target range outside the covariate space it saw.
Predictions are bounded by the training targets. Outside the training feature space the surface flattens rather than failing loudly.

Code examples

Example 1 — fit, validate and report in one function 

import numpy as np
from sklearn.cluster import KMeans
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_absolute_error, r2_score
from sklearn.model_selection import GroupKFold, KFold


def fit_spatial_forest(X, y, coords, feature_names, n_blocks=10,
                       n_estimators=200, min_samples_leaf=2, seed=0):
    """Fit, score both ways, and report the difference."""
    def build():
        return RandomForestRegressor(
            n_estimators=n_estimators, min_samples_leaf=min_samples_leaf,
            max_features="sqrt", n_jobs=-1, random_state=seed)

    def score(splitter):
        predicted = np.empty(len(y))
        for train, test in splitter:
            predicted[test] = build().fit(X[train], y[train]).predict(X[test])
        return r2_score(y, predicted), mean_absolute_error(y, predicted)

    random_r2, random_mae = score(
        KFold(5, shuffle=True, random_state=seed).split(X))
    blocks = KMeans(n_blocks, random_state=seed, n_init=10).fit_predict(coords)
    block_r2, block_mae = score(GroupKFold(5).split(X, y, blocks))

    print(f"  random CV  R2 {random_r2:+.3f}  MAE {random_mae:.4f}")
    print(f"  block  CV  R2 {block_r2:+.3f}  MAE {block_mae:.4f}")
    print(f"  gap {random_r2 - block_r2:+.3f} — the part that depends on "
          "spatial proximity")

    model = build().fit(X, y)
    order = np.argsort(-model.feature_importances_)
    for i in order:
        print(f"    {feature_names[i]:16} {model.feature_importances_[i]:.3f}")

    print(f"  target range {y.min():.3f}..{y.max():.3f} — predictions cannot "
          "leave it")
    return model, {"random_r2": random_r2, "block_r2": block_r2}

Printing the target range next to the model is a small thing that prevents a large misunderstanding. A forest asked to predict a warmer future or a taller canopy than it was trained on cannot do it, and the range makes that concrete.

Example 2 — prediction with ensemble spread 

import numpy as np


def predict_with_spread(model, X, batch=100_000):
    """Mean prediction and the standard deviation across trees."""
    n = len(X)
    mean = np.empty(n, dtype="float32")
    spread = np.empty(n, dtype="float32")

    for start in range(0, n, batch):
        chunk = X[start:start + batch]
        per_tree = np.stack([tree.predict(chunk) for tree in model.estimators_])
        mean[start:start + batch] = per_tree.mean(axis=0)
        spread[start:start + batch] = per_tree.std(axis=0)

    print(f"  prediction {mean.min():.3f}..{mean.max():.3f}")
    print(f"  spread     {spread.min():.3f}..{spread.max():.3f} "
          f"(median {np.median(spread):.3f})")
    print("  note: tree spread measures disagreement, not a calibrated interval")
    return mean, spread

Batching matters because the per-tree stack is n_trees × n_samples floats. For 200 trees and a million cells that is 800 MB in one array.

Example 3 — quantile regression forest for real intervals 

import numpy as np


def quantile_predictions(model, X_train, y_train, X_new, quantiles=(0.1, 0.9),
                         batch=20_000):
    """Empirical quantiles from the training targets in each leaf."""
    leaves_train = model.apply(X_train)              # (n_train, n_trees)
    n_new = len(X_new)
    out = {q: np.empty(n_new) for q in quantiles}

    for start in range(0, n_new, batch):
        chunk = X_new[start:start + batch]
        leaves_new = model.apply(chunk)
        for i in range(len(chunk)):
            matching = np.zeros(len(y_train), bool)
            for tree in range(leaves_new.shape[1]):
                matching |= leaves_train[:, tree] == leaves_new[i, tree]
            values = y_train[matching]
            for q in quantiles:
                out[q][start + i] = np.quantile(values, q) if values.size \
                    else np.nan

    width = out[quantiles[-1]] - out[quantiles[0]]
    print(f"  {int((quantiles[-1] - quantiles[0]) * 100)}% interval width: "
          f"median {np.nanmedian(width):.3f}, "
          f"p95 {np.nanpercentile(width, 95):.3f}")
    return out

A quantile regression forest keeps all the training targets in each leaf rather than only their mean, so it can report empirical quantiles. It is genuinely more informative than the tree spread, and much slower — the loop above is illustrative rather than production code.

For real work, use an implementation that stores the leaf memberships efficiently.

Explanation

Why the forest cannot leave the training range

Each tree predicts the mean of the training targets in a leaf. An average of values from a set is bounded by that set's minimum and maximum, and averaging across trees preserves the bound.

So a forest trained on vegetation indices from 0.2 to 0.9 will never predict 0.95, whatever the covariates say. For interpolation that is a useful conservatism; for any question about conditions outside the training range it is a hard limit.

Linear and gradient-based models extrapolate instead, which is worse in a different way: they produce confident values with no support at all.

Why min_samples_leaf matters more than tree count

With min_samples_leaf=1, a tree can isolate every training sample. Under random cross-validation on spatially correlated data, a held-out sample often has a near neighbour in the training set, so memorisation scores well.

Raising it to two or five forces each leaf to average several samples, which both regularises the model and shrinks it substantially.

Tree count, by contrast, only reduces the variance of the ensemble average. It cannot overfit, and its returns flatten quickly — which is why it is the parameter people tune and the one that matters least.

Why tree spread is not a confidence interval

The spread across trees measures how much the trees disagree, which depends on the bootstrap samples and the feature subsampling as much as on genuine uncertainty.

In particular, it is small where the trees agree — including where they all agree on a wrong answer because none of them saw relevant training data. In an extrapolation region, spread can be low precisely where confidence should be lowest.

Use it as a relative map, pair it with an applicability mask, and use quantile regression forests where a real interval is needed.

Why the two CV numbers are both worth reporting

The gap between random and block cross-validation — 0.555 against 0.221 here — measures how much the model relies on spatial proximity.

A small gap means the covariates are doing the work and the model has some chance of transferring. A large gap means much of the skill comes from being near a training point, which is fine for gap-filling and not for anything else.

Reporting both, with the deployment task stated, is more informative than either alone.

Random forest parameters ranked, with min_samples_leaf having the largest effect and n_estimators the smallest.
The default `min_samples_leaf=1` memorises the training set, which scores well under random CV.

Edge cases or notes

  • Predictions are bounded by the training target range. Report the range.
  • min_samples_leaf=1 memorises; use 2–5 for a few thousand samples.
  • n_estimators cannot overfit and plateaus quickly. 200 is usually enough.
  • max_features="sqrt" decorrelates the trees; the regression default of all features often overfits.
  • Tree spread is not a calibrated interval and can be low in extrapolation regions.
  • Batch the per-tree predictions; the stack is trees × samples.
  • Report both CV numbers and the gap between them.
  • Ship an applicability mask — a forest fails flat, not loudly.

FAQ

How do I predict a continuous surface with a random forest?

Sample covariates at the training points, fit a RandomForestRegressor, validate with a scheme matching the deployment, then predict over the covariate stack and write the result with an applicability mask.

How many trees should I use?

About 200. More trees cannot overfit but the returns flatten quickly, and memory and prediction time scale linearly.

What should min_samples_leaf be?

Two to five for a few thousand samples. The default of 1 memorises the training set, which scores well under random cross-validation on spatial data.

Can a random forest extrapolate?

No. Predictions are bounded by the training target range, and outside the training feature space the surface saturates at a constant.

Is the spread across trees an uncertainty estimate?

Only a relative one. It measures disagreement, and it can be low in extrapolation regions where all trees agree on an unsupported answer.

How do I get real prediction intervals?

Quantile regression forests, which keep the training targets in each leaf rather than only their mean.

Which cross-validation number should I report?

Both, with the deployment task stated. The gap between them measures how much the model depends on spatial proximity.