Geostatistical Modeling Spatial Statistics with Python

Combining Trend Models with Kriging Residuals

Fit a trend with model.fit(X, y), take residuals = y - model.predict(X), krige them with OrdinaryKriging(x, y, residuals, ...), then add model.predict(X_grid) + kriged_residuals. That two-line idea — trend plus kriged residual — is the whole of regression kriging. The steps below make it runnable end to end and report a leave-one-out RMSE.

Why This Matters

Ordinary kriging interpolates a variable from its own spatial autocorrelation alone, ignoring the covariates — elevation, temperature, distance-to-coast — that often explain most of its variation. Regression kriging fixes that by letting a regression model absorb the covariate-driven trend and reserving kriging for the leftover spatial structure. This page is the hands-on companion to the Regression Kriging guide and sits inside the broader Kriging, Interpolation & Surface Generation Techniques workflow. Get the residual step right and you typically beat both plain regression and plain kriging.

Environment

bash
pip install numpy==1.26.4 scipy==1.13.0 scikit-learn==1.5.0 pykrige==1.7.2 matplotlib==3.9.0
python
import numpy as np                                  # 1.26.4
from sklearn.ensemble import RandomForestRegressor  # scikit-learn 1.5.0
from sklearn.linear_model import LinearRegression
from pykrige.ok import OrdinaryKriging              # pykrige 1.7.2
import matplotlib.pyplot as plt                     # 3.9.0

Python 3.10 or 3.11 is recommended. Pin numpy<2 if your pykrige build predates NumPy 2 support.

Step-by-Step Implementation

Step 1 — Build a Covariate-Driven Sample

We simulate a field whose mean depends on two covariates plus a spatially autocorrelated residual, so we know regression kriging should help. In practice, replace this block with your GeoDataFrame and covariate rasters sampled at the points.

python
rng = np.random.default_rng(42)
n = 200

# Sample locations in a projected CRS (metres)
x = rng.uniform(0, 10000, n)
y = rng.uniform(0, 10000, n)

# Two covariates available everywhere (e.g. elevation, temperature)
elev = 200 + 0.03 * y + rng.normal(0, 15, n)
temp = 18 - 0.0009 * x + rng.normal(0, 0.5, n)

# Spatially autocorrelated residual (smooth bump field)
cx, cy = 4000, 6000
resid_true = 3.0 * np.exp(-((x - cx) ** 2 + (y - cy) ** 2) / (2 * 2500 ** 2))

# Observed target: covariate trend + spatial residual + noise
z = 0.02 * elev - 0.4 * temp + resid_true + rng.normal(0, 0.2, n)

X = np.column_stack([elev, temp])   # design matrix for the trend

Step 2 — Fit the Trend Model

Fit a Random Forest for a non-linear trend, or swap in LinearRegression for the OLS (KED-equivalent) version. Keep the forest shallow so it does not soak up the spatial residual.

python
def fit_trend(X, z, kind="rf"):
    if kind == "rf":
        model = RandomForestRegressor(
            n_estimators=400, max_depth=5, min_samples_leaf=5, random_state=0
        )
    else:
        model = LinearRegression()
    model.fit(X, z)
    return model

trend = fit_trend(X, z, kind="rf")
m_hat = trend.predict(X)
print(f"Trend R^2: {trend.score(X, z):.3f}")

Step 3 — Compute and Check Residuals

The residuals are what kriging will interpolate. Confirm their mean is near zero — a biased mean signals a mis-specified trend.

python
residuals = z - m_hat
print(f"Residual mean: {residuals.mean():+.4f}")
print(f"Residual std : {residuals.std():.4f}")

Step 4 — Variogram and Ordinary-Krige the Residuals

Fit an ordinary kriging model to the residuals. PyKrige fits the residual variogram automatically; inspect the parameters to confirm real structure (a rising variogram, not pure nugget).

python
ok_res = OrdinaryKriging(
    x, y, residuals,
    variogram_model="spherical",
    nlags=12,
    weight=True,
    verbose=False,
)
psill, vrange, nugget = ok_res.variogram_model_parameters
print(f"Residual variogram -> psill={psill:.3f}, range={vrange:.0f} m, nugget={nugget:.3f}")

# Prediction grid in the same metric units
gridx = np.linspace(0, 10000, 120)
gridy = np.linspace(0, 10000, 120)
resid_grid, resid_var = ok_res.execute("grid", gridx, gridy)

Step 5 — Predict the Trend on the Grid and Sum

Evaluate the trend at every grid cell (from exhaustive covariate rasters in production), then add the kriged residual surface.

python
# Reconstruct the covariates on the grid. In production, read these from rasters.
GX, GY = np.meshgrid(gridx, gridy)
elev_grid = 200 + 0.03 * GY
temp_grid = 18 - 0.0009 * GX
X_grid = np.column_stack([elev_grid.ravel(), temp_grid.ravel()])

trend_grid = trend.predict(X_grid).reshape(GY.shape)
rk_surface = trend_grid + np.asarray(resid_grid)   # final regression kriging map

print(f"RK surface range: {rk_surface.min():.2f} to {rk_surface.max():.2f}")

Step 6 — Cross-Validate and Report RMSE

Leave-one-out cross-validation refits the trend and residual kriging on all-but-one point each fold, then predicts the held-out point. We compare regression kriging against ordinary kriging of the raw values on the same folds.

python
def loocv_rmse(x, y, z, X, kind="rf"):
    n = len(z)
    err_rk, err_ok = np.zeros(n), np.zeros(n)
    idx = np.arange(n)
    for i in range(n):
        tr = idx != i
        # --- Regression kriging ---
        model = fit_trend(X[tr], z[tr], kind=kind)
        res_tr = z[tr] - model.predict(X[tr])
        ok = OrdinaryKriging(x[tr], y[tr], res_tr,
                             variogram_model="spherical", nlags=12,
                             weight=True, verbose=False)
        res_pred, _ = ok.execute("points", x[i:i+1], y[i:i+1])
        rk_pred = model.predict(X[i:i+1])[0] + float(res_pred[0])
        err_rk[i] = rk_pred - z[i]
        # --- Ordinary kriging of raw values (baseline) ---
        ok_raw = OrdinaryKriging(x[tr], y[tr], z[tr],
                                 variogram_model="spherical", nlags=12,
                                 weight=True, verbose=False)
        ok_pred, _ = ok_raw.execute("points", x[i:i+1], y[i:i+1])
        err_ok[i] = float(ok_pred[0]) - z[i]
    return np.sqrt(np.mean(err_rk ** 2)), np.sqrt(np.mean(err_ok ** 2))

rmse_rk, rmse_ok = loocv_rmse(x, y, z, X, kind="rf")
print(f"Regression kriging LOOCV RMSE : {rmse_rk:.4f}")
print(f"Ordinary kriging   LOOCV RMSE : {rmse_ok:.4f}")
print(f"Improvement: {100 * (1 - rmse_rk / rmse_ok):.1f}%")

Interpreting the Output

  • Trend R2R^2: the share of target variance the covariates explain. A moderate-to-high R2R^2 is what makes regression kriging pay off; near zero means the covariates are uninformative and you should fall back to ordinary kriging.
  • Residual variogram parameters: a range on the order of the true residual correlation length and a non-trivial partial sill confirm the residuals still carry spatial structure worth kriging. A pure-nugget residual variogram means the trend already captured everything.
  • LOOCV RMSE: the headline number. Regression kriging should report a lower RMSE than ordinary kriging on this covariate-driven field; the printed improvement quantifies the gain. If regression kriging is not better, the covariates are not helping and the extra complexity is unjustified.
  • Residual kriging variance (resid_var): uncertainty of the residual component only — near zero at samples, larger in gaps.

Critical Best Practices

Do Not Overfit the Trend

A deep Random Forest can drive training residuals to near zero, collapsing the variogram to pure nugget and neutralising the kriging step. Cap max_depth, raise min_samples_leaf, and confirm the residual variogram still rises with lag. The trend should explain the smooth, covariate-driven pattern, not memorise every point.

Refit Everything Inside the CV Loop

Fitting the trend on all the data and then cross-validating only the kriging leaks information and flatters the RMSE. Refit both the trend and the residual variogram on the training subset of every fold, as the loop above does, so the reported error reflects genuine out-of-sample performance.

Keep Residuals Stationary

The residual kriging step assumes a zero-mean, second-order stationary residual field. Check the residual mean and scan for leftover coordinate trends; if a gradient remains, add a covariate or move toward the drift-based approach in Ordinary & Universal Kriging.

Use Exhaustive Covariate Rasters for the Grid

The trend must be evaluated at every prediction cell, so each covariate needs a raster covering the full extent, co-registered to one grid and CRS. Covariates that exist only at sample points cannot drive a regression kriging surface.

Troubleshooting

Symptom Likely cause Fix
RK RMSE not better than OK Covariates uninformative or trend overfitted Reassess covariates; regularise the trend; confirm residual variogram has structure
Residual variogram is flat (pure nugget) Trend absorbed the spatial signal Use a shallower/less flexible trend model
LinAlgError during kriging Duplicate coordinates in a fold Deduplicate points or add tiny coordinate jitter
Residual mean far from zero Trend mis-specified or wrong scale Add covariates, transform the target, or use OLS with an intercept
Grid surface has implausible extremes Trend extrapolating beyond covariate range Clip covariates to the training range or use a bounded regressor
CV is very slow Refitting kriging n times on many points Use fewer nlags, moving-window kriging, or k-fold instead of LOOCV

Next Steps

Return to the Regression Kriging guide for the theory and the KED contrast, then decide whether a one-stage drift model fits better via Ordinary vs Universal Kriging: Which to Use.

Frequently Asked Questions

Should I use OLS or Random Forest for the regression kriging trend?

Use OLS when covariate effects are roughly linear and interpretability matters; it is equivalent to universal kriging with external drift. Use a shallow Random Forest when the covariate-target relationship is non-linear or involves interactions. Keep the forest regularised (limited depth, minimum leaf size) so residuals retain spatial structure for the kriging step.

Why krige the residuals instead of the raw values?

The trend already explains the large-scale, covariate-driven pattern. Kriging the residuals interpolates only what the trend missed — the spatially autocorrelated fluctuation. This usually lowers prediction error versus ordinary kriging of the raw values whenever the covariates are informative, and it keeps the residual field closer to the stationarity kriging assumes.

How do I validate a regression kriging model?

Use leave-one-out or spatially blocked cross-validation. Refit the trend and residual variogram on the training subset each fold, predict the held-out point, and accumulate errors into an RMSE. Compare against ordinary kriging on the same folds; regression kriging should win when covariates carry signal.

Related

← Back to Regression Kriging