Regression Kriging: Combining Trend & Residual Interpolation

Regression kriging (RK) is the workhorse of environmental and soil mapping because it lets you exploit cheap, wall-to-wall covariates — elevation, temperature, distance-to-coast, spectral indices — to explain the large-scale pattern of a variable, while still honouring the fine-scale spatial autocorrelation that pure regression ignores. It decomposes the field into a deterministic trend plus a stochastic residual, models each with the tool best suited to it, and sums them back together. This page develops the method within the wider Kriging, Interpolation & Surface Generation Techniques workflow, with annotated scikit-learn and pykrige code, residual diagnostics, and production scaling notes.

Concept Summary #

Regression kriging models an observed variable $Z(s)$ at location $s$ as the sum of a deterministic trend and a spatially correlated residual, then predicts each part separately and adds them. It is the practical answer to a recurring problem: ordinary kriging captures spatial autocorrelation but wastes the explanatory power of covariates, while ordinary regression captures covariate effects but throws away the spatially structured residual. Regression kriging keeps both. It is a close relative of the drift-based estimators covered in Ordinary & Universal Kriging, and the residual step depends entirely on a well-fitted variogram.

Prerequisites #

• Python 3.10+
• numpy>=1.24, scipy>=1.10, scikit-learn>=1.3, pykrige==1.7.2, geopandas>=0.14, rasterio>=1.3
• Optional: scikit-gstat>=1.0 for flexible residual variogram fitting
• Target observations as a GeoDataFrame in a projected (metric) CRS
• Covariates available as exhaustive rasters covering the full prediction extent (not just sample points), co-registered to a common grid and CRS

Mathematical Core #

Regression kriging rests on the decomposition of the target field into a mean function and a residual:

$$Z(s) = m(s) + \varepsilon(s)$$

where $m(s)$ is the deterministic trend — the conditional expectation of $Z$ given a vector of covariates $\mathbf{x}(s) = [x_1(s), \dots, x_p(s)]$ — and $\varepsilon(s)$ is a zero-mean, second-order stationary residual capturing the spatial autocorrelation the trend does not explain.

The trend is estimated by any regression that maps covariates to the target. For an ordinary least-squares (OLS) trend,

$$\hat{m}(s) = \beta_0 + \sum_{k=1}^{p} \beta_k \, x_k(s),$$

though $\hat{m}$ may equally be a Random Forest or gradient-boosting prediction. The residuals at the sample locations are

$$\hat{\varepsilon}(s_i) = Z(s_i) - \hat{m}(s_i), \qquad i = 1, \dots, n.$$

These residuals are interpolated with ordinary kriging, which requires their variogram. The empirical semivariance at lag $h$ is

$$\hat{\gamma}(h) = \frac{1}{2 \, |N(h)|} \sum_{(i,j) \in N(h)} \big[ \hat{\varepsilon}(s_i) - \hat{\varepsilon}(s_j) \big]^2,$$

where $N(h)$ is the set of point pairs separated by approximately $h$. A theoretical model (spherical, exponential, Gaussian) is fitted to $\hat{\gamma}(h)$, giving nugget, sill, and range. The kriged residual at an unsampled location $s_0$ is a weighted linear combination

$$\hat{\varepsilon}(s_0) = \sum_{i=1}^{n} \lambda_i \, \hat{\varepsilon}(s_i),$$

with weights $\lambda_i$ solved from the ordinary kriging system built on the residual variogram. The final regression kriging prediction is simply the sum of the two independently predicted components:

$$\hat{Z}(s_0) = \hat{m}(s_0) + \hat{\varepsilon}(s_0).$$

Relationship to Universal Kriging (KED) #

Universal kriging, also called kriging with external drift (KED), embeds the trend directly into the kriging system by adding unbiasedness constraints on the covariates, so the trend coefficients and the interpolation weights are solved together in one linear system. Regression kriging separates the two stages. When the trend is linear in the covariates and both approaches use the same residual variogram, the predictions are algebraically identical. The practical differences drive the choice:

The one subtlety in universal kriging is that its residual variogram should be estimated from the detrended residuals, not the raw data — a step that regression kriging makes explicit and hard to skip. That transparency is the main reason practitioners reach for regression kriging when covariates matter.

Annotated Implementation #

The implementation below is a complete regression kriging pipeline. It fits a trend, computes residuals, kriges them, and sums the two surfaces. Each function has an explicit signature so you can drop it into a larger workflow.

1. Sample Covariates at Observation Points #

The trend model needs one row per observation, with the target value and the covariate values sampled at that point. When covariates arrive as rasters, sample them with rasterio.

import numpy as np
import geopandas as gpd
import rasterio


def sample_covariates(gdf: gpd.GeoDataFrame,
                      raster_paths: dict[str, str]) -> np.ndarray:
    """Sample each covariate raster at the point geometries of gdf.

    Returns an (n_points, n_covariates) design matrix aligned with gdf order.
    All rasters must share the CRS of gdf.
    """
    coords = [(geom.x, geom.y) for geom in gdf.geometry]
    columns = []
    for name, path in raster_paths.items():
        with rasterio.open(path) as src:
            sampled = np.array([v[0] for v in src.sample(coords)], dtype=float)
        columns.append(sampled)
    return np.column_stack(columns)


# Observations in a projected CRS (e.g. UTM zone 32N)
gdf = gpd.read_file("soil_samples.gpkg").to_crs("EPSG:32632")
raster_paths = {
    "elevation": "dem.tif",
    "temperature": "lst.tif",
    "dist_coast": "dist_coast.tif",
}

X = sample_covariates(gdf, raster_paths)      # trend design matrix
y = gdf["organic_carbon"].to_numpy(dtype=float)
x_coord = gdf.geometry.x.to_numpy()
y_coord = gdf.geometry.y.to_numpy()

2. Fit the Deterministic Trend #

Any scikit-learn regressor works. OLS is interpretable and gives the classic KED-equivalent model; a Random Forest captures non-linear covariate effects at the cost of harder residual behaviour.

from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor


def fit_trend(X: np.ndarray, y: np.ndarray, kind: str = "ols"):
    """Fit and return a trend regressor mapping covariates X to target y."""
    if kind == "ols":
        model = LinearRegression()
    elif kind == "rf":
        # Shallow forest: avoid overfitting so residuals keep spatial structure
        model = RandomForestRegressor(
            n_estimators=300, max_depth=6, min_samples_leaf=5, random_state=0
        )
    else:
        raise ValueError("kind must be 'ols' or 'rf'")
    model.fit(X, y)
    return model


trend = fit_trend(X, y, kind="ols")
m_hat = trend.predict(X)          # trend estimate at sample points
residuals = y - m_hat             # what kriging must interpolate

3. Diagnose Residual Stationarity #

Before kriging, confirm the residuals behave like a zero-mean stationary field. A non-zero mean or a visible spatial gradient means the trend is mis-specified.

print(f"Residual mean : {residuals.mean():.4f}  (target ~0)")
print(f"Residual std  : {residuals.std():.4f}")
print(f"Trend R^2     : {trend.score(X, y):.3f}")

# Quick spatial-drift check: correlate residuals with coordinates.
# Strong correlation implies leftover trend the variogram cannot absorb.
for name, c in (("easting", x_coord), ("northing", y_coord)):
    r = np.corrcoef(residuals, c)[0, 1]
    print(f"corr(residual, {name}) = {r:+.3f}")

A near-zero mean and weak coordinate correlations indicate the residuals are ready for kriging. If a gradient remains, add the offending covariate, use a more flexible trend, or move to the drift-based estimators in the ordinary and universal kriging guide.

4. Model the Residual Variogram and Krige #

Krige the residuals with pykrige. The variogram here describes the residuals, not the raw target — this is the defining feature of regression kriging.

from pykrige.ok import OrdinaryKriging

# Fit an ordinary kriging model to the RESIDUALS
ok_res = OrdinaryKriging(
    x_coord, y_coord, residuals,
    variogram_model="spherical",
    nlags=15,
    weight=True,
    verbose=False,
    enable_plotting=False,
)
print("Residual variogram [psill, range, nugget]:",
      ok_res.variogram_model_parameters)

If you prefer to fit the residual variogram explicitly (for anisotropy control or robust estimators), scikit-gstat exposes the full model and can hand parameters back to pykrige.

5. Predict the Trend on the Grid and Sum #

The trend must be evaluated at every prediction cell, which is why exhaustive covariate rasters are mandatory. Read each raster block, build the per-cell feature matrix, predict the trend, krige residuals on the same grid, and add.

def build_grid_features(raster_paths: dict[str, str]):
    """Read co-registered rasters into a (n_cells, n_covariates) matrix
    plus the grid coordinate axes. Assumes identical transforms."""
    stack, gx, gy = [], None, None
    for path in raster_paths.values():
        with rasterio.open(path) as src:
            arr = src.read(1).astype(float)
            if gx is None:
                rows, cols = arr.shape
                xs = src.xy(np.zeros(cols), np.arange(cols))[0]
                ys = src.xy(np.arange(rows), np.zeros(rows))[1]
                gx, gy, shape = np.array(xs), np.array(ys), arr.shape
        stack.append(arr.ravel())
    return np.column_stack(stack), gx, gy, shape


Xg, gridx, gridy, shape = build_grid_features(raster_paths)

trend_grid = trend.predict(Xg).reshape(shape)          # m(s) everywhere
resid_grid, resid_var = ok_res.execute("grid", gridx, gridy)  # kriged residuals

rk_prediction = trend_grid + np.asarray(resid_grid)    # final RK surface

The kriging variance resid_var describes the uncertainty of the residual component; combining it with the trend model’s own prediction uncertainty is covered under Uncertainty & Variance Mapping.

Output Interpretation #

• Trend $R^2$: how much of the target’s variance the covariates explain. A high $R^2$ (say > 0.6) means the trend carries most of the signal and the residual surface is a fine correction. A low $R^2$ means the covariates are weak and regression kriging degenerates toward ordinary kriging of the raw values.
• Residual variogram: a residual variogram that is nearly pure nugget (flat, no rise with lag) means the trend already removed the spatial structure — kriging the residuals will add little, which is fine and expected when covariates are strong. A clear sill-and-range structure means genuine autocorrelation remains for kriging to exploit.
• Final surface range: predictions should sit within a plausible physical range. Values far outside the observed target range usually trace to an extrapolating trend model (especially OLS on covariates outside their training range), not the kriging step.
• Residual kriging variance: near zero at sample points, rising in data-sparse zones. It reflects only the residual component and understates total uncertainty when the trend itself is uncertain.

Production Considerations #

Kriging cost. Ordinary kriging of residuals inverts an $n \times n$ covariance matrix, so cost scales as $\mathcal{O}(n^3)$. Beyond a few thousand points, use moving-window kriging (n_closest_points with backend="loop"), or fit the variogram with scikit-gstat and krige with a neighbourhood search.

Raster alignment. The single most common production failure is mis-registered covariate rasters. Resample every covariate to one reference grid, CRS, and nodata mask before sampling or gridding. A one-pixel offset silently corrupts both the trend fit and the grid prediction.

Trend overfitting. A deep Random Forest can drive training residuals to near zero, collapsing the variogram to pure nugget and defeating the purpose of kriging. Regularise the trend (shallow depth, minimum leaf size) and always validate on held-out points, ideally with spatially blocked folds rather than random splits.

Memory on large grids. Predicting the trend on a national-scale raster produces huge arrays. Process the grid in windowed blocks with rasterio’s windowed reads and write predictions incrementally to a GTiff.

Troubleshooting #

Next Steps #

Work through the full runnable pipeline — Random Forest trend, residual variogram, and a reported RMSE — in Combining Trend Models with Kriging Residuals. To decide whether a trend even warrants the two-stage approach, compare the estimators in Ordinary vs Universal Kriging: Which to Use, and turn the residual variance into calibrated intervals with Uncertainty & Variance Mapping.

Frequently Asked Questions #

What is the difference between regression kriging and universal kriging?

Regression kriging fits the trend and the residual variogram in two separate, explicit stages: a regression model estimates the deterministic mean from covariates, then the residuals are kriged. Universal kriging (kriging with external drift, KED) solves the trend and the spatial interpolation simultaneously inside one kriging system. Mathematically the two are equivalent when the trend is linear in the covariates and the same variogram is used, but regression kriging lets you use any regressor (Random Forest, gradient boosting) for the trend and inspect the residuals directly.

Can I use a Random Forest as the trend model in regression kriging?

Yes. Random Forest regression kriging (RFRK) is a common variant: the Random Forest captures non-linear covariate effects for the trend, and the residuals are kriged as usual. The key requirement is that the residuals must remain approximately second-order stationary, so avoid overfitting the trend (which drives residuals toward zero at sample points and destroys the spatial structure the variogram needs).

Why must regression kriging residuals be stationary?

The residual kriging step assumes the residuals form a zero-mean, second-order stationary random field so a single variogram describes their covariance everywhere. If the trend model leaves systematic drift in the residuals, the variogram is biased and the kriged residual surface will be unreliable. Always inspect residual maps and residual variograms before kriging, and test for stationarity if a spatial pattern remains.

How do I get covariate values on the full prediction grid?

Regression kriging requires the covariates to be available as exhaustive rasters covering the whole prediction area, because the trend must be evaluated at every grid cell, not just at sample points. Read the covariate rasters with rasterio, stack them into a per-cell feature matrix, apply the fitted regressor, then add the kriged residuals. Covariates that exist only at sample points cannot be used for regression kriging.

Related

• Combining Trend Models with Kriging Residuals — the full step-by-step RK pipeline with a Random Forest trend and cross-validated RMSE
• Ordinary & Universal Kriging — the drift-based estimator that regression kriging generalises
• Uncertainty & Variance Mapping — combining residual kriging variance with trend uncertainty
• Variogram Modeling & Semivariance Analysis — fitting the residual variogram that drives the kriging step

← Back to Kriging, Interpolation & Surface Generation Techniques