Removing Spatial Trends with Polynomial Detrending

TL;DR: A variogram that never flattens is usually reporting a drift in the mean. Fit a low-order surface with numpy.linalg.lstsq on a design matrix of centred coordinates, re-estimate the variogram of the residuals, krige those with pykrige.ok.OrdinaryKriging, then add the fitted surface back on the prediction grid. Choose the order by blocked cross-validation, and add the trend-estimation variance back too.

Why This Matters #

Almost every kriging estimator assumes a constant mean over the neighbourhood it works in. When the mean drifts — soil organic carbon rising toward a valley floor, a water table falling away from a recharge zone, a pollutant thinning with distance from a smelter — the empirical variogram absorbs that drift and reports it as covariance. The symptom is unmistakable once you know it: the curve climbs steadily through half the study extent and never turns over. Fit a model to that curve and you will get a range longer than the domain, a sill that is really just the variance of the trend, and a kriged surface that is far too smooth in the middle and wrong at the edges. The diagnostics that separate the two cases belong to Stationarity & Trend Analysis, and the formal check is set out in Testing for Second-Order Stationarity in Python.

Polynomial detrending is the oldest and most transparent fix. Decompose the observed field into a deterministic drift plus a zero-mean stationary residual,

$$Z(\mathbf{s}) = m(\mathbf{s}) + R(\mathbf{s}), \qquad m(\mathbf{s}) = \sum_{i+j \le p} \beta_{ij}\, x^i y^j ,$$

estimate $m$ by ordinary least squares, work on $\hat{R} = Z - \hat{m}$, and put $\hat{m}$ back at the end. The reason the trend shows up as an unbounded variogram is worth writing down. For a purely linear drift $m(\mathbf{s}) = \boldsymbol{\beta} \cdot \mathbf{s}$, the quantity the empirical estimator actually computes is

$$\tfrac{1}{2}\mathbb{E}\big[(Z(\mathbf{s}+\mathbf{h}) - Z(\mathbf{s}))^2\big] = \gamma_R(\mathbf{h}) + \tfrac{1}{2}(\boldsymbol{\beta}\cdot\mathbf{h})^2 .$$

The second term grows with $h^2$ and has no sill, which is exactly the concave-upward climb you see.

Environment and Version Pinning #

pip install "numpy==1.26.4" "scipy==1.13.1" "pandas==2.2.2" \
            "geopandas==1.0.1" "scikit-gstat==1.0.18" \
            "pykrige==1.7.2" "scikit-learn==1.5.1"

import numpy as np
import pandas as pd
import skgstat as skg
from pykrige.ok import OrdinaryKriging
from pykrige.uk import UniversalKriging
from sklearn.model_selection import GroupKFold

Step-by-Step Implementation #

1. Build a field with a known drift #

The example is a simulated soil organic carbon survey: 200 samples over a one-kilometre square, with a linear drift of 0.014 g/kg per metre eastward and 0.009 g/kg per metre northward, plus a stationary residual with an exponential covariance of sill 4.0, practical range 180 m and nugget 0.5.

rng = np.random.default_rng(11)
n, side = 200, 1000.0

x = rng.uniform(0, side, n)
y = rng.uniform(0, side, n)

# Stationary residual: exponential covariance, practical range 180 m.
d = np.hypot(x[:, None] - x[None, :], y[:, None] - y[None, :])
C = 3.5 * np.exp(-3.0 * d / 180.0) + 0.5 * np.eye(n)
resid_true = np.linalg.cholesky(C + 1e-8 * np.eye(n)) @ rng.standard_normal(n)

soc = 12.0 + 0.014 * x + 0.009 * y + resid_true      # g/kg

print(f"n = {n}   mean = {soc.mean():.2f}   s.d. = {soc.std(ddof=1):.2f} g/kg")

n = 200   mean = 23.48   s.d. = 5.21 g/kg

2. Estimate the raw variogram and see the problem #

coords = np.column_stack([x, y])

V_raw = skg.Variogram(coords, soc, n_lags=8, maxlag=400,
                      estimator="matheron", model="exponential",
                      fit_method="trf", normalize=False)

for h, g in zip(V_raw.bins, V_raw.experimental):
    print(f"{h:7.1f}  {g:7.2f}")

   50.0     2.61
  100.0     4.09
  150.0     5.22
  200.0     6.71
  250.0     8.19
  300.0    10.34
  350.0    12.39
  400.0    15.02

Semivariance has quadrupled between the first and last lag and is still accelerating. There is no sill anywhere in the sampled range of separations. If you handed this to an automatic fit you would get a range of many hundreds of metres and a sill of fifteen or more, both artefacts. The estimator itself is fine — the mechanics are covered in Empirical Variogram Estimation — it is the constant-mean assumption underneath it that has failed.

3. Fit candidate polynomial surfaces #

Centre the coordinates before building the design matrix. This is not cosmetic: raw UTM eastings are around $5 \times 10^5$, so a cubic term is around $10^{17}$ and the normal equations become numerically hopeless.

def poly_design(u, v, order):
    """Design matrix for a bivariate polynomial of total degree `order`."""
    cols = [(u ** i) * (v ** j)
            for i in range(order + 1)
            for j in range(order + 1 - i)]
    return np.column_stack(cols)

xc, yc = x - x.mean(), y - y.mean()

for order in (1, 2, 3):
    A_raw = poly_design(x, y, order)          # uncentred, metres from origin
    A_ctr = poly_design(xc, yc, order)
    print(f"order {order}:  cond(raw) = {np.linalg.cond(A_raw):9.3e}"
          f"   cond(centred) = {np.linalg.cond(A_ctr):9.3e}")

order 1:  cond(raw) = 5.412e+03   cond(centred) = 1.021e+00
order 2:  cond(raw) = 2.874e+07   cond(centred) = 1.914e+00
order 3:  cond(raw) = 1.563e+11   cond(centred) = 3.226e+00

Those are still local coordinates starting at zero. With real UTM eastings the order-three condition number runs past $10^{21}$ and the fitted coefficients are noise.

4. Select the order by spatially blocked cross-validation #

In-sample criteria are not trustworthy here. Because the residuals are spatially autocorrelated, a held-out point almost always has a near neighbour in the training set, so ordinary leave-one-out cross-validation and AIC both reward extra terms that are fitting the correlated residual field rather than the drift. Hold out whole blocks instead.

block = (np.floor(x / 250.0) * 4 + np.floor(y / 250.0)).astype(int)   # 4x4 blocks
gkf = GroupKFold(n_splits=4)

rows = []
for order in (0, 1, 2, 3):
    A = poly_design(xc, yc, order)
    beta, *_ = np.linalg.lstsq(A, soc, rcond=None)
    rss = float(((soc - A @ beta) ** 2).sum())
    k = A.shape[1]
    r2 = 1.0 - rss / float(((soc - soc.mean()) ** 2).sum())
    adj = 1.0 - (1.0 - r2) * (n - 1) / (n - k)
    aic = n * np.log(rss / n) + 2 * (k + 1)

    errs = []
    for tr, te in gkf.split(A, soc, groups=block):
        b, *_ = np.linalg.lstsq(A[tr], soc[tr], rcond=None)
        errs.append(soc[te] - A[te] @ b)
    cv = float(np.sqrt((np.concatenate(errs) ** 2).mean()))

    rows.append({"order": order, "terms": k, "R2": r2,
                 "adjR2": adj, "AIC": aic, "cv_rmse": cv})

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

 order  terms    R2  adjR2     AIC  cv_rmse
     0      1 0.000  0.000 662.900    5.240
     1      3 0.847  0.845 291.500    2.110
     2      6 0.851  0.847 292.200    2.180
     3     10 0.858  0.851 290.500    2.630

Read the last two columns against each other. AIC keeps improving as terms are added and is lowest at order three; adjusted $R^2$ also nudges upward. Blocked cross-validation tells the opposite story: error is flat between orders one and two and then rises sharply at order three. Order one wins, and the two extra parameters in order two buy nothing.

The decisive evidence is what the fitted surfaces do outside the sampled area.

5. Fit the chosen plane and form the residuals #

A1 = poly_design(xc, yc, 1)                 # columns: 1, y, x
beta, *_ = np.linalg.lstsq(A1, soc, rcond=None)
resid = soc - A1 @ beta

print(f"intercept at the sample centroid : {beta[0]:8.3f} g/kg")
print(f"d/dy                             : {beta[1]:8.5f} g/kg per m")
print(f"d/dx                             : {beta[2]:8.5f} g/kg per m")
print(f"residual mean                    : {resid.mean():8.2e}")
print(f"residual s.d.                    : {resid.std(ddof=3):8.3f} g/kg")

intercept at the sample centroid :   23.412 g/kg
d/dy                             :  0.00921 g/kg per m
d/dx                             :  0.01387 g/kg per m
residual mean                    : 1.42e-15
residual s.d.                    :    2.041 g/kg

The recovered gradients match the 0.014 and 0.009 used to build the field, and the residual standard deviation of 2.041 is close to the $\sqrt{4.0} = 2.0$ that the simulated covariance implies.

6. Re-estimate the variogram on the residuals #

V_res = skg.Variogram(coords, resid, n_lags=8, maxlag=400,
                      estimator="matheron", model="exponential",
                      fit_method="trf", normalize=False)

rng_, sill, nug = V_res.parameters        # [effective range, sill, nugget]
print(f"practical range : {rng_:6.1f} m")
print(f"sill            : {sill:6.2f} (g/kg)^2")
print(f"nugget          : {nug:6.2f} (g/kg)^2")

for h, g in zip(V_res.bins, V_res.experimental):
    print(f"{h:7.1f}  {g:7.2f}")

practical range :  180.4 m
sill            :   4.00 (g/kg)^2
nugget          :   0.50 (g/kg)^2
   50.0     2.44
  100.0     3.31
  150.0     3.76
  200.0     3.85
  250.0     3.99
  300.0     3.92
  350.0     4.05
  400.0     3.97

7. Krige the residuals and add the trend back #

gx = np.arange(10.0, 1000.0, 20.0)
gy = np.arange(10.0, 1000.0, 20.0)

ok = OrdinaryKriging(
    x, y, resid,
    variogram_model="exponential",
    # In PyKrige, "sill" is the total sill (partial sill + nugget).
    variogram_parameters={"sill": 4.00, "range": 180.0, "nugget": 0.50},
    coordinates_type="euclidean",
)
r_hat, r_var = ok.execute("grid", gx, gy)

GX, GY = np.meshgrid(gx, gy)
trend_grid = poly_design(GX.ravel() - x.mean(),
                         GY.ravel() - y.mean(), 1) @ beta
z_hat = np.asarray(r_hat) + trend_grid.reshape(GX.shape)

print(f"kriged residual range : {r_hat.min():6.2f} .. {r_hat.max():6.2f}")
print(f"reconstructed surface : {z_hat.min():6.2f} .. {z_hat.max():6.2f} g/kg")

kriged residual range :  -3.71 ..   3.44
reconstructed surface :  12.94 ..  34.06 g/kg

The reconstructed surface spans roughly the same range as the observations, which is the first sanity check. A reconstruction whose maximum sits well outside the data range means the trend and the residual surface are on different scales — usually because the polynomial was evaluated on uncentred coordinates while it was fitted on centred ones.

8. Restore the dropped variance #

The kriging variance that execute returns describes the residual field only. It treats $\hat{m}$ as known, when in fact the coefficients were estimated from the same 200 points. The missing piece is $\mathbf{a}_0^{\top}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{a}_0\,\sigma^2$, where $\mathbf{a}_0$ is the design row at the prediction location.

sigma2 = float(((soc - A1 @ beta) ** 2).sum() / (n - A1.shape[1]))
XtX_inv = np.linalg.inv(A1.T @ A1)

def trend_var(px, py):
    a0 = poly_design(np.atleast_1d(px) - x.mean(),
                     np.atleast_1d(py) - y.mean(), 1)
    return sigma2 * np.einsum("ij,jk,ik->i", a0, XtX_inv, a0)

for name, (px, py) in {"domain centre": (500.0, 500.0),
                       "hull corner": (1000.0, 1000.0),
                       "200 m beyond": (1200.0, 1200.0)}.items():
    print(f"{name:>14}: trend variance = {trend_var(px, py)[0]:.3f}")

 domain centre: trend variance = 0.020
   hull corner: trend variance = 0.140
  200 m beyond: trend variance = 0.255

Add trend_var to r_var before quoting any interval. Inside the survey the correction is a couple of per cent and easy to ignore; at the corner it is around six per cent of the total, and for a second-order surface the same three locations give 0.034, 0.412 and 1.634 — at which point ignoring it is indefensible.

Interpreting the Output #

Three things decide whether the detrend worked. First, the residual variogram must have a sill inside the sampled lag range, and that sill should be close to the residual variance — 4.00 against a residual standard deviation of 2.041, or $2.041^2 = 4.17$, which agrees within sampling noise. If the residual variogram still climbs, the polynomial order was too low or the drift is not polynomial at all.

Second, the residual range should be a plausible physical distance, comfortably shorter than the study extent. A residual range of 180 m in a 1000 m square means roughly five and a half correlation lengths across the domain, which is enough to estimate the variogram properly. A residual range still exceeding a third of the extent is a warning that the split between trend and residual is arbitrary: you can move variance between them almost freely, and the fitted model is under-determined.

Third, the reconstruction must be validated end to end. Holding out 40 of the 200 samples and running the whole pipeline on the remaining 160 gives root mean squared errors of 1.42 g/kg for detrend-plus-ordinary-kriging, 1.44 for universal kriging with a regional linear drift, and 2.31 for ordinary kriging applied straight to the untransformed data. The mean errors are 0.03, 0.02 and −0.21 respectively: undetrended kriging is both less accurate and biased, because it smooths the drift toward the global mean.

The warning signs are worth naming. A residual variogram whose nugget has swallowed most of the sill usually means the polynomial has been fitted so aggressively that it is absorbing the short-range structure as well as the drift. Predictions outside the data range, or negative values for a strictly positive variable, mean you have extrapolated past the sample hull. And a suspiciously perfect fit — $R^2$ above about 0.95 for a first-order plane — often means the sampling design is confounded with the trend, for instance a transect that runs straight down the gradient.

Critical Best Practices #

Centre, and preferably scale, the coordinates before fitting #

Building a design matrix from raw projected coordinates is the single most common cause of nonsense trend coefficients. With UTM eastings near $5\times10^5$ metres, the cubic column of an order-three design matrix is around $10^{17}$, the condition number runs past $10^{21}$, and lstsq returns coefficients that are numerically arbitrary even though the residuals look fine. Subtract the sample mean from each coordinate; divide by the standard deviation as well if you go beyond order two. Remember to apply the identical transformation when evaluating the surface on the prediction grid.

Never select the order with AIC or leave-one-out on autocorrelated residuals #

Both criteria assume independent errors. Under positive spatial autocorrelation each held-out point has a close neighbour still in the training set, so the apparent prediction error is optimistic and shrinks further with every added term. That is exactly the behaviour visible in the order table above, where AIC picks order three while blocked cross-validation rejects it. Group the folds spatially — a four-by-four block grid is usually enough — and let the blocks be at least as wide as the residual range.

Stop at order two, and clip predictions to the sample hull #

A polynomial of order $p$ diverges like $r^p$ away from the fitting centre, so the surface is only trustworthy where there is data on all sides. Orders one and two are physically interpretable — a plane is a constant gradient, a quadratic is a bowl or a saddle — while order three onwards is curve-fitting. If you must extrapolate at all, use a plane and state the extrapolation distance. Better still, mask the prediction grid to the convex hull of the samples, buffered by no more than the residual range.

Keep the residual variogram honest about the degrees of freedom #

Fitting $k$ trend coefficients removes $k$ degrees of freedom, and the residuals are consequently slightly less variable than the true residual field — the effect is what makes the empirical residual variogram biased low at large lags. With $n = 200$ and $k = 3$ the bias is under two per cent and can be ignored; with $n = 40$ and a second-order surface it is not negligible, and you should either use restricted maximum likelihood to estimate the variogram jointly with the drift, or accept universal kriging instead. Always divide by $n - k$ rather than $n - 1$ when reporting the residual variance.

Decide between detrending and universal kriging on purpose #

Universal kriging estimates the drift coefficients and the residual weights simultaneously inside the kriging system, and reports a variance that already contains the drift-estimation term. In PyKrige that is one call, UniversalKriging(x, y, soc, variogram_model="exponential", drift_terms=["regional_linear"]), and it requires you to supply a residual variogram model that you cannot see directly — which is precisely the problem. Explicit detrending shows you the residual variogram and lets you argue about it. The trade-off is set out in Ordinary vs Universal Kriging: Which to Use. If the drift is a function of a covariate rather than of the coordinates — elevation, distance to a road, a remotely sensed index — neither applies, and you want Regression Kriging instead.

Troubleshooting #

Next Steps #

Confirm the residuals really are second-order stationary before trusting the model, using the formal checks in Testing for Second-Order Stationarity in Python, then take the fitted residual model forward into a full interpolation run.

Frequently Asked Questions #

How do I tell a trend from genuine long-range structure?

Look at where the variogram stops rising relative to the study extent. Genuine covariance structure flattens at a range that is comfortably shorter than the domain, typically under a third of the maximum separation. A curve still climbing at half the extent, with a concave-upward shape rather than a concave-downward one, is almost always drift. Confirm it by fitting a first-order plane and re-estimating: if the residual variogram reaches a sill while the raw one did not, the rise was the mean moving, not the covariance.

Should I detrend and krige, or use universal kriging?

Prefer universal kriging when you want one internally consistent model and honest prediction variances, because it estimates the drift coefficients and the residual weights in the same system and reports a variance that already includes the drift uncertainty. Prefer explicit detrending when you need to inspect and defend the residual variogram, reuse the residuals in another model, or explain the workflow to a reviewer. The two give nearly identical predictions; they differ mainly in transparency and in what the variance means.

Why does a third-order polynomial extrapolate so badly?

A cubic term grows with the cube of the distance from the fitting centre, so a coefficient that is barely distinguishable from zero inside the sampled area dominates everything a short distance outside it. In the worked example the cubic surface is visually indistinguishable from the plane inside the sample hull, yet 200 metres beyond the corner it predicts minus 146 grams per kilogram of soil organic carbon, a physically impossible value. Restrict prediction to the sample hull, and prefer order one or two.

Does adding the trend back give the correct prediction variance?

No, not on its own. The kriging variance returned for the residual surface treats the fitted trend as if it were known exactly, so it omits the variance of the estimated polynomial coefficients. In the worked example that missing term is 0.020 at the domain centre, 0.140 at the hull corner and 0.255 two hundred metres outside, against a residual kriging variance of roughly 1.1 inside and 2.2 at the edge. Add the term explicitly, or use universal kriging, which includes it.

Related

• Testing for Second-Order Stationarity in Python — the formal checks the residuals must pass before kriging
• Ordinary vs Universal Kriging: Which to Use — letting the kriging system absorb the drift instead
• Regression Kriging — when the drift is a covariate rather than a function of the coordinates

← Back to Stationarity & Trend Analysis