Universal Kriging with External Drift in PyKrige

TL;DR: Kriging with external drift in PyKrige is UniversalKriging(..., drift_terms=["specified"], specified_drift=[elev_pts]) followed by uk.execute("grid", gridx, gridy, specified_drift_arrays=[drift_grid]). The drift must exist at the observations and at every prediction cell, in the grid’s own row order, with no NaN. Fit the variogram to the regression residuals, not the raw data.

Why This Matters #

Universal kriging in its textbook form models the trend as a low-order polynomial in the coordinates. That is a confession of ignorance: it says the mean varies smoothly across the map but offers no reason why. When you have an auxiliary variable that is known everywhere — elevation for air temperature, distance to the nearest road for a traffic pollutant, modelled rainfall for a soil property — you can use it as the trend directly and let the kriging system estimate only the coefficient. The result is kriging with external drift, also written KED, and it is by some distance the most useful member of the Ordinary & Universal Kriging family for environmental work.

The practical gain is that the interpolated surface inherits the auxiliary variable’s spatial detail. Two hundred weather stations cannot resolve a mountain valley; a 30 m digital elevation model can, and if temperature really does fall with altitude then a temperature surface built through the elevation drift will show the valley too. What the method demands in return is discipline about alignment, because a single misordered array turns a defensible surface into an expensive picture. The related choice of when a trend is present at all belongs to Stationarity & Trend Analysis, and the decision between polynomial and auxiliary trends is set out in Ordinary vs Universal Kriging: Which to Use.

The Model and Why the Drift Is Needed Everywhere #

Write the observed field as a deterministic trend built from $p+1$ known basis functions plus a zero-mean, second-order stationary residual:

$$Z(\mathbf{s}) = \sum_{k=0}^{p} a_k f_k(\mathbf{s}) + \varepsilon(\mathbf{s}), \qquad f_0(\mathbf{s}) \equiv 1 .$$

For external drift with a single auxiliary variable $q$, $p = 1$ and $f_1(\mathbf{s}) = q(\mathbf{s})$. The coefficients $a_k$ are never estimated separately — they are eliminated by imposing universality constraints on the kriging weights. The predictor is the usual linear combination $\hat{Z}(\mathbf{s}_0) = \sum_i \lambda_i Z(\mathbf{s}_i)$ subject to

$$\sum_{i=1}^{n} \lambda_i f_k(\mathbf{s}_i) = f_k(\mathbf{s}_0), \qquad k = 0, \dots, p ,$$

which gives the bordered system

$$\begin{bmatrix} \boldsymbol{\Gamma} & \mathbf{F} \\ \mathbf{F}^{\mathsf{T}} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \boldsymbol{\lambda} \\ \boldsymbol{\mu} \end{bmatrix} = \begin{bmatrix} \boldsymbol{\gamma}_0 \\ \mathbf{f}_0 \end{bmatrix}, \qquad \sigma^2_{\text{KED}}(\mathbf{s}_0) = \sum_i \lambda_i \gamma(\mathbf{s}_i, \mathbf{s}_0) + \sum_k \mu_k f_k(\mathbf{s}_0) .$$

Look at the right-hand side. The constraint for $k = 1$ is literally the number $q(\mathbf{s}_0)$: the drift value at the point you are predicting. There is no way to assemble the system without it, which is the entire reason the drift has to be known at every prediction location and not merely at the samples. It is also why an auxiliary variable measured only at the sample points cannot be an external drift; that case is Regression Kriging territory instead.

Note that $\gamma$ in the system is the variogram of $\varepsilon$, not of $Z$. That distinction drives everything on this page.

The four drift terms PyKrige offers #

"external_Z" and "specified" do the same statistical job. The difference is who does the resampling: "external_Z" accepts a regularly spaced raster plus its axis coordinates and interpolates it for you, while "specified" accepts values you have already extracted and does nothing except use them. Prefer "specified" in production, because the resampling then happens in rasterio where you can inspect it, choose the resampling method, and see the nodata cells before they reach the linear algebra.

"point_log" is a genuinely different animal, intended for logarithmic drawdown around pumping wells or point emission sources, and it is not interchangeable with an external drift raster.

Environment and Version Pinning #

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

import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
from rasterio.windows import from_bounds
import skgstat as skg
from pykrige.uk import UniversalKriging

PyKrige 1.7 is the first release in which pseudo_inv=True is stable for badly conditioned drift systems, which matters here because an elevation column and a column of ones can be close to collinear when the sample elevations span a narrow band.

Step-by-Step Implementation #

The worked example is 120 automatic weather stations in a mountainous catchment, projected to EPSG:32633, with a 200 m digital elevation model covering an 80 km by 64 km frame.

1. Sample the drift at the observation points #

DEM = "dem_200m.tif"
stations = gpd.read_file("stations.gpkg")          # columns: temp_c, geometry

with rasterio.open(DEM) as src:
    if stations.crs != src.crs:
        stations = stations.to_crs(src.crs)
    coords = list(zip(stations.geometry.x, stations.geometry.y))
    elev_pts = np.array([v[0] for v in src.sample(coords)], dtype="float64")
    if src.nodata is not None:
        elev_pts[elev_pts == src.nodata] = np.nan

x = stations.geometry.x.to_numpy(dtype="float64")
y = stations.geometry.y.to_numpy(dtype="float64")
temp = stations["temp_c"].to_numpy(dtype="float64")

assert np.isfinite(elev_pts).all(), "station outside the DEM or on a nodata cell"
print(f"stations: {len(stations)}  |  CRS: {stations.crs.to_string()}")
print(f"elevation: {elev_pts.min():.1f} to {elev_pts.max():.1f} m, "
      f"NaN {np.isnan(elev_pts).sum()}")
print(f"temperature: {temp.min():.2f} to {temp.max():.2f} degC")

stations: 120  |  CRS: EPSG:32633
elevation: 41.0 to 1948.3 m, NaN 0
temperature: 6.14 to 18.86 degC

src.sample is a generator yielding one array per point with one element per band, so the v[0] extracts band 1. It does not reproject, which is why the CRS check comes first — sampling latitude and longitude against a UTM raster silently returns the nodata value for every point rather than raising.

2. Build the drift on the prediction grid #

BOUNDS = (400_000.0, 5_180_000.0, 480_000.0, 5_244_000.0)   # 80 km x 64 km

with rasterio.open(DEM) as src:
    win = from_bounds(*BOUNDS, transform=src.transform)
    dem = src.read(1, window=win, masked=True).astype("float64")
    tf = src.window_transform(win)

ny, nx = dem.shape
xs = tf.c + tf.a * (np.arange(nx) + 0.5)     # cell centres, ascending
ys = tf.f + tf.e * (np.arange(ny) + 0.5)     # tf.e < 0, so descending

gridx = xs
gridy = ys[::-1]                    # PyKrige is given ascending y
drift_grid = np.flipud(dem.filled(np.nan))   # flip the raster to match

print(f"grid: nx={nx}, ny={ny}, cell={tf.a:.0f} m")
print(f"drift_grid.shape = {drift_grid.shape}  NaN = "
      f"{int(np.isnan(drift_grid).sum())}")
print(f"gridy[0]={gridy[0]:.0f}  gridy[-1]={gridy[-1]:.0f}  (ascending)")

grid: nx=400, ny=320, cell=200 m
drift_grid.shape = (320, 400)  NaN = 0
gridy[0]=5180100  gridy[-1]=5243900  (ascending)

This is the step that goes wrong most often. A north-up GeoTIFF stores row 0 at the northern edge, so its rows run from high $y$ to low $y$. PyKrige returns a grid whose row 0 corresponds to gridy[0]. Give it a descending gridy and everything still runs, but the drift is applied upside down relative to the coordinates and the map comes out mirrored about the horizontal axis. Flipping both together, as above, keeps them consistent whichever convention you choose.

3. Fit the variogram to the residuals, not the data #

import numpy.linalg as la

A = np.column_stack([np.ones_like(elev_pts), elev_pts])
beta, *_ = la.lstsq(A, temp, rcond=None)
resid = temp - A @ beta
r2 = 1.0 - resid.var(ddof=0) / temp.var(ddof=0)

print(f"lapse rate  = {beta[1]:.6f} degC per m")
print(f"intercept   = {beta[0]:.3f} degC")
print(f"R2          = {r2:.3f}")
print(f"residual sd = {resid.std(ddof=2):.3f} degC  (n = {len(resid)})")

lapse rate  = -0.006103 degC per m
intercept   = 19.412 degC
R2          = 0.812
residual sd = 0.949 degC  (n = 120)

V = skg.Variogram(np.column_stack([x, y]), resid,
                  model="spherical", n_lags=15, maxlag=40_000,
                  normalize=False)
vrange, psill, nugget = V.parameters      # scikit-gstat order
sill = psill + nugget

print(f"range   = {vrange:,.0f} m")
print(f"psill   = {psill:.3f}  nugget = {nugget:.3f}  sill = {sill:.3f}")
# the plateau of the fitted model must equal nugget + partial sill
assert abs(V.fitted_model(3 * vrange) - sill) < 1e-9

range   = 17,940 m
psill   = 0.784  nugget = 0.118  sill = 0.902

Two things deserve attention. First, V.parameters returns [range, sill, nugget] where scikit-gstat’s “sill” is the plateau above the nugget — a partial sill. PyKrige’s list form for variogram_parameters wants exactly [psill, range, nugget], so the hand-over is a reorder rather than a conversion, but the dictionary form uses the key sill for the total plateau and derives the partial sill as sill - nugget. Mixing the two conventions silently changes the model by the size of the nugget.

Second, and more important: PyKrige computes its automatic experimental variogram from the raw z you pass in, which still contains the trend. Left to itself, UniversalKriging will fit a model to a variogram that climbs past the sample variance and never flattens, and the weights it derives from that model are wrong. Always pass the residual parameters explicitly.

4. Run the kriging #

uk = UniversalKriging(
    x, y, temp,
    variogram_model="spherical",
    variogram_parameters=[psill, vrange, nugget],   # [psill, range, nugget]
    drift_terms=["specified"],
    specified_drift=[elev_pts],
    exact_values=True,
    pseudo_inv=True,
)

z, ss = uk.execute("grid", gridx, gridy,
                   specified_drift_arrays=[drift_grid])

print(f"z  {z.shape}   ss {ss.shape}")
print(f"prediction: {z.min():.2f} to {z.max():.2f} degC, "
      f"masked {int(np.ma.getmaskarray(z).sum())}")
print(f"kriging variance: min {ss.min():.3f}, mean {ss.mean():.3f}, "
      f"max {ss.max():.3f} degC^2")

z  (320, 400)   ss (320, 400)
prediction: 6.02 to 19.31 degC, masked 0
kriging variance: min 0.124, mean 0.857, max 1.412 degC^2

execute returns masked arrays. The shape is (len(gridy), len(gridx)), which is why drift_grid had to be (ny, nx) and not its transpose — PyKrige will silently transpose a (nx, ny) array if the grid happens to be square, so a square study area hides this bug until you change the extent.

5. Write the surface back out with the right transform #

from rasterio.transform import from_origin

out_tf = from_origin(gridx[0] - tf.a / 2, gridy[-1] - tf.e / 2, tf.a, -tf.e)
profile = dict(driver="GTiff", height=ny, width=nx, count=2,
               dtype="float32", crs=stations.crs, transform=out_tf,
               nodata=np.float32(np.nan), compress="deflate")

with rasterio.open("temp_ked.tif", "w", **profile) as dst:
    dst.write(np.flipud(z.filled(np.nan)).astype("float32"), 1)
    dst.write(np.flipud(np.sqrt(ss.filled(np.nan))).astype("float32"), 2)
    dst.set_band_description(1, "temperature_c")
    dst.set_band_description(2, "kriging_sd_c")

The np.flipud on the way out is the mirror of the one on the way in: PyKrige’s row 0 is the south edge, a GeoTIFF’s row 0 is the north edge. Band 2 stores the kriging standard deviation rather than the variance, because that is the unit readers can compare against the predictions directly.

6. Cross-validate #

errs = np.empty(len(temp))
kvar = np.empty(len(temp))

for i in range(len(temp)):
    keep = np.arange(len(temp)) != i
    fold = UniversalKriging(
        x[keep], y[keep], temp[keep],
        variogram_model="spherical",
        variogram_parameters=[psill, vrange, nugget],
        drift_terms=["specified"],
        specified_drift=[elev_pts[keep]],
        exact_values=True, pseudo_inv=True,
    )
    zi, si = fold.execute("points", x[i:i + 1], y[i:i + 1],
                          specified_drift_arrays=[elev_pts[i:i + 1]])
    errs[i] = temp[i] - float(zi[0])
    kvar[i] = float(si[0])

std_resid = errs / np.sqrt(kvar)
print(f"LOO RMSE        : {np.sqrt((errs ** 2).mean()):.3f} degC")
print(f"mean kriging sd : {np.sqrt(kvar).mean():.3f} degC")
print(f"standardised    : mean {std_resid.mean():+.3f}, "
      f"sd {std_resid.std(ddof=1):.3f}")

LOO RMSE        : 0.968 degC
mean kriging sd : 0.921 degC
standardised    : mean -0.012, sd 1.043

The variogram is deliberately held fixed across folds. Refitting it inside the loop is more honest about total uncertainty but roughly doubles the runtime and makes the folds incomparable; the usual compromise is to fix the variogram and treat the result as a check on the weights rather than on the whole modelling chain.

Interpreting the Output #

z is the predicted field and ss is the kriging variance in the squared units of the target, here $^{\circ}\text{C}^{2}$. Three readings matter.

The minimum variance should sit near the nugget, because at an observation the residual is reproduced exactly and only the measurement-scale component remains. Here ss.min() is 0.124 against a nugget of 0.118; a minimum far above the nugget means no prediction cell fell close to a station, and a minimum of exactly zero means exact_values=True is doing what it says at a cell centred on a sample.

The maximum variance can legitimately exceed the sill. Here it reaches 1.412 against a sill of 0.902, and that excess is the $\sum_k \mu_k f_k(\mathbf{s}_0)$ term: uncertainty in the estimated drift coefficient, which grows as you extrapolate to elevations outside the sampled range. If you see a maximum many times the sill, find where it occurs — it will be the cells with the most extreme drift values, and it is telling you the lapse rate is being extrapolated rather than interpolated.

The standardised residuals are the honest check. Their standard deviation of 1.043 says the reported kriging variance is a fair account of the actual error. Anything above about 1.3 means the variance is optimistic, and the usual culprit is a residual variogram that was fitted to a field still containing a trend.

The right-hand panels show the failure mode that has no visual signature. When the drift and the target share only a broad regional gradient — say a pollutant that tracks urbanisation at the scale of the whole city but not street by street — the regression still reports a high $R^2$, and the map still inherits all the drift raster’s fine texture. That texture is not evidence; it is an assumption applied at a scale where nobody checked it. The residual variogram gives it away by continuing to climb, because the short-range structure the drift failed to explain is still there, and the standardised residuals give it away by coming out at 2.41 instead of 1.04, which means the reported variance understates the squared error by a factor of about six.

Critical Best Practices #

Never let PyKrige fit the variogram for you here #

UniversalKriging estimates its experimental variogram from the values you pass in, before any drift is removed. For a strongly trended field that variogram will not have a sill, the automatic fit will chase it with an implausibly long range and a huge partial sill, and the resulting weights will be far too smooth. Fit the residual model separately and pass variogram_parameters every time you use a drift term. This is the single most consequential difference between PyKrige’s universal kriging and the two-stage approach described in Regression Kriging vs Universal Kriging.

Insist on a linear relationship, and check it before you fit #

External drift kriging assumes the target is a linear function of the auxiliary variable with the same coefficient everywhere. Temperature against elevation satisfies that over a few hundred metres and stops satisfying it across an inversion layer. If the scatter plot bends, transform the drift — $\log$ of distance-to-road rather than distance-to-road is the standard fix for pollution — rather than hoping the kriging will absorb the curvature. It will not: it will push the curvature into the residuals, where it appears as a variogram that climbs without a sill.

Treat nodata as a build failure, not a warning #

A single NaN in specified_drift or specified_drift_arrays propagates through the matrix solve and returns an all-NaN surface with no exception raised, and a nodata sentinel such as -9999 read as a genuine elevation is worse, because it produces a plausible-looking number. Assert np.isfinite(...).all() on both arrays before constructing the object. Where the raster genuinely has holes inside the study area, fill them explicitly — a focal mean or a coarse-resolution fallback — and record how many cells were filled.

Do not extrapolate the drift beyond its sampled range #

The lapse rate here was estimated from stations spanning 41 m to 1948 m. Cells above 1948 m in the DEM are extrapolations of a fitted coefficient, and the kriging variance grows there but does not grow enough, because it accounts for coefficient uncertainty and not for the model being wrong. Mask predictions where the drift lies outside the sampled interval, or at minimum flag them, and say so in the metadata.

Budget memory, because there is no moving window #

OrdinaryKriging.execute accepts n_closest_points to solve a local system per target; UniversalKriging.execute does not, because the drift constraints make a local neighbourhood ill-posed in general. The vectorized backend therefore builds a dense (M, N) matrix — 128 000 cells by 120 stations here, about 123 MB in float64, which is fine, and 40 GB for a 4 000 by 4 000 grid, which is not. Tile the prediction grid into blocks of a few hundred thousand cells, slice drift_grid the same way, and concatenate.

Troubleshooting #

Next Steps #

If the residual variogram refuses to flatten no matter which auxiliary variable you try, the trend is not linear in that variable and the two-stage approach in Regression Kriging lets you fit an arbitrary regression before kriging what is left; the trade-offs are compared directly in Regression Kriging vs Universal Kriging.

Frequently Asked Questions #

Why must the drift be supplied at the prediction locations as well as the data?

The external drift enters the kriging system as an unbiasedness constraint. The weights are forced to reproduce the drift value at the target location, so that value has to be known before the system can be solved. Without it there is no equation to satisfy and PyKrige raises an error rather than guessing. This is also why the drift variable must be exhaustively mapped: an auxiliary variable known only at the sample points cannot be used as an external drift, only as a covariate in regression kriging.

Should the variogram be fitted to the data or to the residuals?

To the residuals. PyKrige fits its experimental variogram to the raw observations, which still contain the trend, so the automatic fit is biased upward and usually shows no sill at all. Regress the target on the drift, fit a model to the residuals with scikit-gstat or PyKrige itself, then pass those parameters to UniversalKriging through variogram_parameters. The residual variogram is also the diagnostic: if it never flattens, the drift has not removed the trend.

What is the difference between specified, external_Z and regional_linear drift terms?

regional_linear adds the coordinates $x$ and $y$ themselves as trend basis functions, so the trend is a plane fitted to the sample locations. external_Z takes a regularly gridded auxiliary raster plus its axis coordinates and bilinearly interpolates it wherever PyKrige needs a value. specified takes the drift values you have already extracted, once at the data points and once for every prediction location, and does no interpolation at all. Use specified when you control the sampling and want the alignment to be explicit.

How do I tell whether the kriging variance can be trusted?

Run leave-one-out cross-validation and divide each prediction error by the kriging standard deviation at that point. The standard deviation of those standardised residuals should be close to one. A value near 2.4, as happens when the drift tracks the target only at long range, means the reported variance understates the real error by roughly a factor of six. The map will still look plausible, because the drift raster supplies convincing fine detail that the sparse observations never tested.

Related

• Ordinary vs Universal Kriging: Which to Use — deciding whether a trend term is warranted at all
• Regression Kriging — fitting the trend separately when it is not linear in the drift
• Stationarity & Trend Analysis — the assumption the residual variogram is testing

← Back to Ordinary & Universal Kriging