How to Predict a Continuous Surface with Random Forest

On this page

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 Sentinel-2 NDVI on Snowdon from terrain covariates, 3,000 training cells (Example 1):

random 5-fold CV        +0.398
block cross-validation  +0.297
true error over the
whole unsampled map     +0.367

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)

block CV R2 +0.258

The mask is part of the product, not an optional extra. A forest cannot signal that it is extrapolating.

Leaf size, tree count, validation, extrapolation, spread

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 can look like skill under random cross-validation — though on this data it did not: leaf 1 scored +0.370 under random CV against +0.377 at leaf 5 (Example 4).

What it certainly does is bloat the model. Two hundred trees on 3,000 samples were 54.6 MB with leaf 1 and 8.8 MB with leaf 5; on 10,000 samples, 182 MB against 29 MB. Two to five is a reasonable range for a few thousand samples.

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 took 119.1 MB against 23.8 MB and 1.9 s against 0.4 s to fit and map, and moved the true R² from +0.348 to +0.349.

Spend the tuning effort on max_features and min_samples_leaf. Here max_features="sqrt" raised every score by 0.02 to 0.04 — the largest change any setting made.

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.

Predictions are bounded by the training targets. Outside the training feature space the surface flattens rather than failing loudly.

Fit, validate and predict with ensemble and quantile intervals

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_grid and y_grid hold every remaining cell — the map to predict — and X, y and coords a random 3,000 of them, the training sample.

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"]
X_grid, y_grid, coords_grid = X, y, coords
sample = np.random.default_rng(0).choice(len(y_grid), 3000, replace=False)
X, y, coords = X_grid[sample], y_grid[sample], coords_grid[sample]
print(f"{item.id}: {shape[1]} x {shape[0]} cells of 20 m, {valid.sum():,} usable "
      f"({valid.mean():.1%}); training on {len(y):,}")

S2C_30UVD_20250712_0_L2A: 543 x 509 cells of 20 m, 270,536 usable (97.9%); training on 3,000

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}


model, scores = fit_spatial_forest(X, y, coords, FEATURES)
rest = np.setdiff1d(np.arange(len(y_grid)), sample)
print(f"  true R2 over the {len(rest):,} unsampled cells {r2_score(y_grid[rest], model.predict(X_grid[rest])):+.3f}")

  random CV  R2 +0.398  MAE 0.0764
  block  CV  R2 +0.297  MAE 0.0841
  gap +0.101 — the part that depends on spatial proximity
    elevation        0.436
    slope            0.197
    tpi              0.136
    northness        0.118
    eastness         0.114
  target range 0.034..0.927 — predictions cannot leave it
  true R2 over the 267,536 unsampled cells +0.367

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


mean, spread = predict_with_spread(model, X_grid)

  prediction 0.235..0.898
  spread     0.016..0.315 (median 0.084)
  note: tree spread measures disagreement, not a calibrated interval

Batching matters because the per-tree stack is n_trees × n_samples floats, and predict returns float64. For 200 trees and a million cells that is 1.6 GB 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


check = np.random.default_rng(5).choice(rest, 500, replace=False)       # 500 unsampled cells
bounds = quantile_predictions(model, X, y, X_grid[check])
inside = (y_grid[check] >= bounds[0.1]) & (y_grid[check] <= bounds[0.9])
print(f"  {inside.mean():.1%} of 500 unsampled cells fall inside their 80% interval")
tree_lo, tree_hi = mean[check] - 1.2816 * spread[check], mean[check] + 1.2816 * spread[check]
print(f"  tree spread read as ±1.28 sd covers {((y_grid[check] >= tree_lo) & (y_grid[check] <= tree_hi)).mean():.1%}")

  80% interval width: median 0.210, p95 0.473
  82.6% of 500 unsampled cells fall inside their 80% interval
  tree spread read as ±1.28 sd covers 78.4%

A quantile regression forest keeps all the training targets in each leaf rather than only their mean, so it can report empirical quantiles. Checked on 500 cells the forest never saw, its 80% interval held 82.6% of them. It is much slower than the tree spread — the loop above is illustrative rather than production code.

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

Measured in Example 3. A check like this, on cells the model never saw, is what turns either into a usable interval.

Example 4 — which parameters matter, measured

import pickle
import time


def scores_for(**params):
    build = lambda: RandomForestRegressor(n_jobs=-1, random_state=0, **params)
    def cv(splitter):
        p = np.empty(len(y))
        for train, test in splitter:
            p[test] = build().fit(X[train], y[train]).predict(X[test])
        return r2_score(y, p)
    blocks = KMeans(10, random_state=0, n_init=10).fit_predict(coords)
    started = time.perf_counter()
    fitted = build().fit(X, y)
    true_r2 = r2_score(y_grid[rest], fitted.predict(X_grid[rest]))
    seconds = time.perf_counter() - started
    return (cv(KFold(5, shuffle=True, random_state=0).split(X)), cv(GroupKFold(5).split(X, y, blocks)),
            true_r2, len(pickle.dumps(fitted)) / 1e6, seconds)


print(f"{'':34} {'random CV':>9} {'block CV':>9} {'true':>7} {'size':>8} {'fit+map':>8}")
for label, params in (("200 trees, min_samples_leaf=1", dict(n_estimators=200, min_samples_leaf=1)),
                      ("200 trees, min_samples_leaf=2", dict(n_estimators=200, min_samples_leaf=2)),
                      ("200 trees, min_samples_leaf=5", dict(n_estimators=200, min_samples_leaf=5)),
                      ("1,000 trees, min_samples_leaf=2", dict(n_estimators=1000, min_samples_leaf=2)),
                      ("200 trees, leaf 2, max_features=sqrt", dict(n_estimators=200, min_samples_leaf=2, max_features="sqrt"))):
    r, b, t, mb, sec = scores_for(**params)
    print(f"{label:34} {r:+9.3f} {b:+9.3f} {t:+7.3f} {mb:6.1f} MB {sec:6.1f} s")

big = np.random.default_rng(1).choice(len(y_grid), 10_000, replace=False)
for leaf in (1, 5):
    forest = RandomForestRegressor(n_estimators=200, min_samples_leaf=leaf, n_jobs=-1, random_state=0)
    print(f"10,000 samples, min_samples_leaf={leaf}: {len(pickle.dumps(forest.fit(X_grid[big], y_grid[big]))) / 1e6:.0f} MB")

print(f"map predictions {mean.min():.3f}..{mean.max():.3f}, training targets {y.min():.3f}..{y.max():.3f}, "
      f"the map itself {y_grid.min():.3f}..{y_grid.max():.3f}")

# Tree spread where the forest has and has not seen the terrain: train in the west, look at both halves
west = coords_grid[:, 0] < np.median(coords_grid[:, 0])
in_west = west[sample]
westward = RandomForestRegressor(n_estimators=200, min_samples_leaf=2, n_jobs=-1, random_state=0).fit(X[in_west], y[in_west])
for label, half in (("west, where it trained", west), ("east, where it did not", ~west)):
    cells = np.random.default_rng(3).choice(np.nonzero(half)[0], 20_000, replace=False)
    per_tree = np.stack([t.predict(X_grid[cells]) for t in westward.estimators_])
    err = np.abs(per_tree.mean(0) - y_grid[cells])
    print(f"{label:24} median tree spread {np.median(per_tree.std(0)):.4f}, median error {np.median(err):.4f}")

                                   random CV  block CV    true     size  fit+map
200 trees, min_samples_leaf=1         +0.370    +0.246  +0.340   54.6 MB    0.4 s
200 trees, min_samples_leaf=2         +0.374    +0.258  +0.348   23.8 MB    0.4 s
200 trees, min_samples_leaf=5         +0.377    +0.275  +0.345    8.8 MB    0.4 s
1,000 trees, min_samples_leaf=2       +0.376    +0.258  +0.349  119.1 MB    1.9 s
200 trees, leaf 2, max_features=sqrt    +0.398    +0.297  +0.367   23.8 MB    0.4 s
10,000 samples, min_samples_leaf=1: 182 MB
10,000 samples, min_samples_leaf=5: 29 MB
map predictions 0.235..0.898, training targets 0.034..0.927, the map itself -0.765..0.960
west, where it trained   median tree spread 0.0816, median error 0.0505
east, where it did not   median tree spread 0.1152, median error 0.0775

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 can score well.

Raising it to two or five forces each leaf to average several samples, which regularises the model and shrinks it six-fold here; block CV rose from +0.246 at leaf 1 to +0.275 at leaf 5.

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. Nothing in its construction stops it being low exactly where confidence should be lowest. On this data it behaved better than that: a forest trained in the western half had a median spread of 0.115 in the east against 0.082 in the west, where its errors were larger too (0.078 against 0.051), and read as a normal interval the spread covered 78.4% of unsampled cells at a nominal 80%. Check it the same way before relying on it.

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.398 against 0.297 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.

Measured in Example 4. max_features changed the fit; min_samples_leaf changed the size; tree count changed neither much.

min_samples_leaf, max_features and batched tree predictions

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.

Internal links

Spatial machine learning explained — why location changes the rules
How to do spatial cross-validation in Python — choosing the scheme
Extrapolation in space explained: the area of applicability — the mask to ship
How to build spatial features for a machine learning model — the covariates
How to turn model predictions back into a raster — writing the surface
How to evaluate a spatial model honestly — what else to report
Feature importance says coordinates are the best predictor — reading the importance table
IDW, kriging, splines or TIN? Choosing an interpolator — the non-machine-learning alternatives

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 and made the model six times larger here, for no gain in accuracy.

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.

How to Predict a Continuous Surface with Random Forest

Problem statement

Quick answer

Leaf size, tree count, validation, extrapolation, spread

1. Set min_samples_leaf above 1

2. Do not tune n_estimators for accuracy

3. Choose the validation scheme from the deployment

4. Remember the forest cannot extrapolate

5. Get uncertainty from the ensemble spread, carefully

Fit, validate and predict with ensemble and quantile intervals

The data used on this page

Example 1 — fit, validate and report in one function

Example 2 — prediction with ensemble spread

Example 3 — quantile regression forest for real intervals

Example 4 — which parameters matter, measured

Why the forest cannot leave the training range

Why the forest cannot leave the training range

Why min_samples_leaf matters more than tree count

Why tree spread is not a confidence interval

Why the two CV numbers are both worth reporting

min_samples_leaf, max_features and batched tree predictions

Internal links

FAQ

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

How many trees should I use?

What should min_samples_leaf be?

Can a random forest extrapolate?

Is the spread across trees an uncertainty estimate?

How do I get real prediction intervals?

Which cross-validation number should I report?

1. Set `min_samples_leaf` above 1

2. Do not tune `n_estimators` for accuracy

Why `min_samples_leaf` matters more than tree count

What should `min_samples_leaf` be?