Spatial Machine Learning vs Kriging

TL;DR: Benchmark both on one fold assignment: GroupKFold over spatial blocks, pykrige.ok.OrdinaryKriging(...).execute("points", x, y) against RandomForestRegressor. With few covariates and strong spatial structure kriging wins; with many informative covariates the ensemble wins; regression kriging — the ensemble as trend, ordinary kriging on its out-of-bag residuals — beat both at every rung of the sweep below.

Why This Matters #

The question “should I krige this or throw a gradient booster at it?” is usually settled by whichever library the analyst already knows, and that is a poor way to pick an estimator. The two families fail in opposite directions, so the wrong choice does not degrade gracefully — it produces a surface with the wrong shape, not merely the wrong numbers. Kriging is an explicit stochastic model. It assumes the target is a realisation of a random field with a covariance structure you fit, and everything else follows: it reproduces the observations exactly, it moves smoothly between them, and it returns a prediction variance derived from the same model. Buy those properties and you also buy the obligations — a credible variogram, roughly stationary structure, and a trend that you have to specify yourself, as set out in ordinary and universal kriging.

A tree ensemble makes no distributional commitment at all. It will happily eat forty covariates, discover interactions and thresholds you never thought to write down, and ignore units. What it will not do is produce a continuous surface, respect the observations it was trained on, or extrapolate: every prediction is an average of training targets, so the output is piecewise constant and bounded by the observed range. Uncertainty is not part of the object; quantile forests bolt an approximation on afterwards. This page benchmarks the two families and their hybrid on identical folds, and it belongs beside the rest of spatial machine learning rather than replacing any of it.

Environment and Version Pinning #

The benchmark needs a field simulator to build a target with known structure, a kriging engine, an ensemble, and a quantile forest for the uncertainty comparison.

pip install "numpy==1.26.4" "scipy==1.13.1" "pandas==2.2.3" \
            "scikit-learn==1.5.2" "pykrige==1.7.2" "gstools==1.6.0" \
            "scikit-gstat==1.0.18" "quantile-forest==1.3.11"

import numpy as np
import pandas as pd
import gstools as gs
from pykrige.ok import OrdinaryKriging
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GroupKFold
from quantile_forest import RandomForestQuantileRegressor

Step-by-Step Implementation #

1. Build one dataset with a known split between covariate signal and spatial signal #

A fair comparison needs a target whose composition you control. The construction below mixes a covariate-driven component with an exponential random field in a ratio $\alpha$, so that $\alpha = 0$ is pure spatial structure and $\alpha = 0.85$ is almost entirely covariate-driven.

rng = np.random.default_rng(2026)
N, SIDE = 600, 20_000.0                      # 600 samples over a 20 km square

xs = rng.uniform(0, SIDE, N)
ys = rng.uniform(0, SIDE, N)

# The spatially structured component: one exponential field, range 3 km.
field = gs.SRF(gs.Exponential(dim=2, var=1.0, len_scale=3000.0),
               seed=20260807)((xs, ys))
field = (field - field.mean()) / field.std()

# Eight smooth covariates, each its own field at its own range.
LENS = (1200, 1800, 2500, 3200, 4000, 5000, 6000, 8000)
COV = np.column_stack([
    gs.SRF(gs.Gaussian(dim=2, var=1.0, len_scale=L), seed=100 + i)((xs, ys))
    for i, L in enumerate(LENS)
])

def build_target(k, alpha, noise=0.35, seed=0):
    """Mix k covariates (share alpha) with the spatial field (share 1-alpha)."""
    r = np.random.default_rng(seed)
    if k:
        g = np.tanh(1.2 * COV[:, :k]).sum(axis=1) + 0.8 * COV[:, 0] * COV[:, 1]
        g = (g - g.mean()) / g.std()
    else:
        g = np.zeros(N)
    return (np.sqrt(alpha) * g + np.sqrt(1.0 - alpha) * field
            + r.normal(0.0, noise, N))

z = build_target(k=2, alpha=0.25, seed=1)
print(f"n = {N}, var(z) = {z.var():.3f}")

n = 600, var(z) = 1.118

The tanh and the product term matter: they give the ensemble something non-linear to find that a linear trend surface could not, which is the case where machine learning is supposed to pay for itself.

2. Assign spatial blocks and freeze one fold split #

Random k-fold cross-validation on autocorrelated data leaks neighbours across the split and flatters every method that can memorise location — which is exactly the ensemble. Blocked folds are the only honest referee here; the reasoning is set out in full under cross-validation strategies.

BLOCK = 4000.0                                # 4 km blocks -> 25 blocks
block = (np.floor(xs / BLOCK).astype(int) * 100
         + np.floor(ys / BLOCK).astype(int))

folds = list(GroupKFold(n_splits=5).split(COV, z, groups=block))
print(f"{pd.Series(block).nunique()} blocks -> test sizes "
      f"{[len(te) for _, te in folds]}")

25 blocks -> test sizes [128, 122, 120, 118, 112]

Every method below consumes this same folds list. Nothing else in the comparison is allowed to vary.

3. Ordinary kriging, refitted inside each fold #

The variogram is part of the model, so it must be refitted on the training rows of each fold. Fitting it once on all 600 points and reusing it leaks the test blocks into the model and inflates the kriging score by a few percent.

def krige_fold(tr, te, target):
    ok = OrdinaryKriging(
        xs[tr], ys[tr], target[tr],
        variogram_model="exponential",
        nlags=15,
        exact_values=True,
        coordinates_type="euclidean",
        enable_plotting=False,
    )
    pred, var = ok.execute("points", xs[te], ys[te])
    return np.asarray(pred), np.asarray(var)

def scores(obs, pred):
    e = np.asarray(pred) - np.asarray(obs)
    ss_tot = np.sum((obs - obs.mean()) ** 2)
    return {"rmse": float(np.sqrt(np.mean(e ** 2))),
            "mae": float(np.mean(np.abs(e))),
            "r2": float(1.0 - np.sum(e ** 2) / ss_tot)}

def run_folds(fn, target):
    obs, pred = [], []
    for tr, te in folds:
        obs.append(target[te])
        pred.append(fn(tr, te, target))
    return scores(np.concatenate(obs), np.concatenate(pred))

ok_scores = run_folds(lambda tr, te, t: krige_fold(tr, te, t)[0], z)

4. The tree ensemble on the same folds #

The forest gets the covariates plus the raw coordinates, which is the minimum spatial feature set; richer geometry is covered in random forest with spatial features and spatial CV.

FEAT = np.column_stack([COV, xs, ys])

def rf_fold(tr, te, target, n_cov=2):
    cols = list(range(n_cov)) + [COV.shape[1], COV.shape[1] + 1]
    rf = RandomForestRegressor(n_estimators=500, min_samples_leaf=3,
                               max_features=0.5, n_jobs=-1, random_state=0)
    rf.fit(FEAT[np.ix_(tr, cols)], target[tr])
    return rf.predict(FEAT[np.ix_(te, cols)])

rf_scores = run_folds(lambda tr, te, t: rf_fold(tr, te, t, n_cov=2), z)

5. Regression kriging: ensemble trend, kriged residuals #

The trend model is fitted on covariates only — no coordinates. If the forest can see the coordinates it will model the spatial field itself, the residuals will be white noise, and the kriging step will have nothing left to do. The residuals must also be out-of-bag: in-bag residuals of a 500-tree forest are far smaller than its true errors, which drags the fitted sill down and makes the residual kriging almost inert.

def rk_fold(tr, te, target, n_cov=2):
    cols = list(range(n_cov))
    trend = RandomForestRegressor(n_estimators=500, min_samples_leaf=3,
                                  max_features=0.5, bootstrap=True,
                                  oob_score=True, n_jobs=-1, random_state=0)
    trend.fit(COV[np.ix_(tr, cols)], target[tr])

    resid = target[tr] - trend.oob_prediction_       # honest residuals
    ok = OrdinaryKriging(xs[tr], ys[tr], resid,
                         variogram_model="exponential", nlags=15,
                         exact_values=False)          # residuals are noisy
    r_hat, _ = ok.execute("points", xs[te], ys[te])
    return trend.predict(COV[np.ix_(te, cols)]) + np.asarray(r_hat)

rk_scores = run_folds(lambda tr, te, t: rk_fold(tr, te, t, n_cov=2), z)

print(pd.DataFrame({"ordinary kriging": ok_scores,
                    "random forest": rf_scores,
                    "regression kriging": rk_scores}).T
        .to_string(float_format=lambda v: f"{v:.3f}"))

                     rmse    mae     r2
ordinary kriging    0.482  0.381  0.792
random forest       0.611  0.483  0.666
regression kriging  0.455  0.359  0.815

With two covariates carrying a quarter of the signal, kriging is comfortably ahead of the forest and the hybrid is ahead of both.

6. Sweep the covariate strength and find the crossover #

One dataset proves nothing about which family to prefer. Re-run the whole benchmark while the covariates take over more of the target.

RUNGS = [(0, 0.00), (2, 0.25), (4, 0.50), (6, 0.70), (8, 0.85)]
rows = []
for k, alpha in RUNGS:
    zk = build_target(k=k, alpha=alpha, seed=1)
    rows.append({
        "k": k, "alpha": alpha,
        "OK":  run_folds(lambda tr, te, t: krige_fold(tr, te, t)[0], zk)["rmse"],
        "RF":  run_folds(lambda tr, te, t: rf_fold(tr, te, t, n_cov=max(k, 1)), zk)["rmse"],
        "RK":  run_folds(lambda tr, te, t: rk_fold(tr, te, t, n_cov=max(k, 1)), zk)["rmse"],
    })

print(pd.DataFrame(rows).to_string(index=False,
                                   float_format=lambda v: f"{v:.3f}"))

 k  alpha     OK     RF     RK
 0  0.000  0.470  0.905  0.468
 2  0.250  0.482  0.611  0.455
 4  0.500  0.598  0.548  0.489
 6  0.700  0.742  0.531  0.502
 8  0.850  0.914  0.523  0.498

Interpreting the Output #

Read the sweep as three separate stories. Ordinary kriging degrades monotonically from 0.470 to 0.914 because it is being asked to interpolate a target that is increasingly not a smooth field: as $\alpha$ rises the residual spatial structure thins out, the fitted sill is progressively dominated by nugget, and kriging degenerates towards predicting the mean. The forest moves the other way, from a hopeless 0.905 with nothing but coordinates to 0.523 with eight covariates and an interaction. At $k = 0$ the forest is not a spatial model in any meaningful sense — it can only carve the coordinate plane into rectangles and average within them, which is a crude and biased interpolator.

Regression kriging is flat, between 0.455 and 0.502 across the whole range, and it is the lowest curve at every rung including the two where kriging alone is strongest. That is the practical headline: the hybrid does not have to be right about which regime it is in. When covariates are weak the trend model contributes almost nothing and the kriging step carries the prediction; when the field is weak the residual variogram flattens and the trend carries it. Details of the fitting sequence are in regression kriging.

Two warning signs are worth naming. If the forest beats kriging by a wide margin under random k-fold but the gap closes or reverses under blocked folds, the forest was memorising location, not learning the process — a near-universal artefact when coordinates are among the features. And if regression kriging scores no better than its own trend model, inspect the residual variogram before concluding anything: a residual sill that is nearly all nugget means there is no structure left to krige, and the extra machinery is buying nothing.

The uncertainty comparison is a separate test, and both methods fail it in different degrees.

tr, te = folds[0]
_, kvar = krige_fold(tr, te, z)
half = 1.645 * np.sqrt(np.maximum(kvar, 0.0))
kpred, _ = krige_fold(tr, te, z)
cov_k = np.mean((z[te] >= kpred - half) & (z[te] <= kpred + half))

cols = [0, 1, COV.shape[1], COV.shape[1] + 1]
qrf = RandomForestQuantileRegressor(n_estimators=500, min_samples_leaf=5,
                                    random_state=0)
qrf.fit(FEAT[np.ix_(tr, cols)], z[tr])
q = qrf.predict(FEAT[np.ix_(te, cols)], quantiles=[0.05, 0.95])
cov_q = np.mean((z[te] >= q[:, 0]) & (z[te] <= q[:, 1]))

print(f"kriging 90% interval : coverage {cov_k:.3f}, mean width {2*half.mean():.2f}")
print(f"quantile forest 90%  : coverage {cov_q:.3f}, mean width {(q[:,1]-q[:,0]).mean():.2f}")

kriging 90% interval : coverage 0.862, mean width 1.55
quantile forest 90%  : coverage 0.741, mean width 1.21

Neither reaches nominal coverage on held-out blocks, which is expected — the kriging variance conditions on a variogram treated as known, and the quantile forest cannot widen its intervals in regions it has never seen. But the kriging interval is close enough to be useful with a modest inflation factor, while the forest interval is short by fifteen coverage points. If a decision depends on the interval rather than the point prediction, that difference usually settles the choice on its own.

Critical Best Practices #

Refit the variogram inside every fold #

A variogram fitted on all the data and then reused across folds has seen the test blocks. The leak is small in absolute terms and reliably positive in direction, so it makes kriging look better than it is against an ensemble that was refitted honestly. Put OrdinaryKriging(...) inside the fold loop, not outside it.

Never give the regression-kriging trend the coordinates #

This is the failure that produces a hybrid indistinguishable from its own trend model. A forest with x and y among its features will split on them repeatedly and absorb the spatial field, leaving residuals with a pure-nugget variogram and a kriging step that adds noise. Give the trend covariates only, and keep the coordinates for the standalone ensemble comparison.

Use out-of-bag residuals, not in-bag ones #

A 500-tree forest fits its training data very closely, so z[tr] - trend.predict(COV[tr]) is far smaller than the model’s real error. A variogram fitted to those residuals has a sill perhaps a third of the correct value, and the kriged correction is scaled down to match. oob_score=True and trend.oob_prediction_ cost nothing and give residuals with roughly the right magnitude.

Match the block size to the autocorrelation range, not to a round number #

Blocks smaller than the variogram range still leak: a test point 500 m from a training point in a field with a 3 km range is not independent. Estimate the range first, then set the block edge at or above it. Here the range is 3 km and the blocks are 4 km. Blocks far larger than the range are safe but waste data by forcing every method to extrapolate further than it will in production.

Score the interval, not only the point prediction #

RMSE ranks estimators on their conditional mean and says nothing about whether their uncertainty is usable. Compute empirical coverage of the nominal interval on the same blocked folds. A method that wins on RMSE while under-covering by fifteen points is the wrong choice for any workflow where the output feeds a threshold, a permit decision or a risk map.

Troubleshooting #

Next Steps #

Build the hybrid properly with the full fitting sequence in regression kriging, and harden the referee that decided this comparison using cross-validation strategies before you report any of these numbers to someone else.

Frequently Asked Questions #

When does kriging beat a random forest?

Kriging wins when most of the target variance is spatially structured and the covariates you hold are weak or few. In the benchmark on this page, with two covariates carrying a quarter of the signal, ordinary kriging scored a blocked-CV RMSE of 0.482 against the forest’s 0.611. Kriging is also the better choice when you need an exact interpolator, a defensible prediction variance, or predictions on a grid where no covariate raster exists at the required resolution.

Why can a tree ensemble not extrapolate?

Every prediction from a forest is an average of training target values held in the leaves an observation falls into. An average of numbers drawn from the training set can never lie outside the training range, so the model is structurally incapable of predicting a value above the highest or below the lowest observation it was trained on. Kriging has no such bound: it forms a weighted combination with weights that can be negative, so it can overshoot local data, and away from the samples it decays smoothly towards the fitted mean.

Can I trust quantile forest intervals the way I trust the kriging variance?

Not by default. The kriging variance is derived from the fitted covariance model and the sampling geometry, so it grows in data-poor areas even though it ignores the uncertainty in the variogram itself. Quantile forest intervals are empirical spreads of training targets and are systematically too narrow where the model must extrapolate. In the benchmark here the nominal 90 percent kriging interval achieved 0.862 coverage on held-out blocks, while the quantile forest interval achieved only 0.741. Calibrate any ensemble interval on blocked folds before publishing it.

Does regression kriging always win?

It won at every rung of the sweep on this page, but it is not automatic. Regression kriging only helps when the trend model leaves genuine spatial structure in its residuals; if the residual variogram is a pure nugget, the kriging step adds fitting cost and nothing else. It also fails quietly when the trend is fitted with coordinates among the features, because the ensemble then absorbs the spatial signal itself and the residuals are white by construction.

Related

• Random Forest with Spatial Features and Spatial CV — the ensemble side of this comparison, built properly
• Regression Kriging — the hybrid that won every rung of the sweep
• Ordinary & Universal Kriging — the kriging side, including when a universal trend replaces the ensemble

← Back to Spatial Machine Learning