How to Fit a Variogram and Krige a Surface in Python

Problem statement

Kriging is inverse distance weighting with the weights derived from measured spatial structure instead of from a typed-in exponent. That difference is worth about 25% accuracy on a well-behaved variable — measured on 500 elevation samples against 12,000 held-out cells:

kriging                40.9 m RMSE
TIN, linear            49.1 m
IDW, power 2           56.1 m
nearest neighbour      65.4 m

It also returns a variance for every prediction, which is the reason most people reach for it. That variance needs careful reading. On the same run:

mean kriging variance   4,450 m²  ->  predicted sd 66.7 m
actual RMSE                                        40.9 m

The kriging variance overstated the real error by 63%. It is a map of where the sample geometry is weak, not a confidence interval.

Quick answer

Ordinary kriging in about forty lines, with no dependency beyond NumPy and SciPy:

import numpy as np
from scipy.spatial import cKDTree


def spherical(h, nugget, sill, rng):
    h = np.asarray(h, dtype=float)
    inside = nugget + sill * (1.5 * h / rng - 0.5 * (h / rng) ** 3)
    return np.where(h == 0, 0.0, np.where(h <= rng, inside, nugget + sill))


def ordinary_kriging(points, values, targets, params, k=24):
    """Predictions and variances from a fitted variogram model."""
    tree = cKDTree(points)
    _, neighbours = tree.query(targets, k=min(k, len(points)))
    neighbours = neighbours.reshape(len(targets), -1)
    k = neighbours.shape[1]

    predicted = np.empty(len(targets))
    variance = np.empty(len(targets))
    cache = {}

    for i, idx in enumerate(neighbours):
        key = idx.tobytes()
        if key not in cache:
            local = points[idx]
            distances = np.sqrt(((local[:, None, :] - local[None, :, :]) ** 2).sum(-1))
            A = np.ones((k + 1, k + 1))
            A[:k, :k] = spherical(distances, *params)
            A[k, k] = 0.0                       # Lagrange multiplier row
            cache[key] = np.linalg.pinv(A)

        b = np.ones(k + 1)
        b[:k] = spherical(np.sqrt(((points[idx] - targets[i]) ** 2).sum(-1)), *params)
        w = cache[key] @ b
        predicted[i] = w[:k] @ values[idx]
        variance[i] = w @ b

    return predicted, variance

DEM 52 to 1072 m; 500 samples, 12,000 held out

A checked variogram, a local neighbourhood, the variance

1. Fit the variogram, and check the fit before using it

Kriging is only as good as the variogram. Fit a model, then look at three things: does the empirical curve reach a plateau inside your largest lag, does the fitted range sit inside it, and is the nugget a plausible fraction of the sill?

nugget      0 m²    (0.0% of the total)
sill   43,743 m²    data variance 38,896 m²
range   3,910 m     largest observed lag 5,875 m

The range is inside the data, which is what you want. The sill exceeds the data variance by 12% — expected here, because the range is a large fraction of the study extent. See The variogram explained.

2. Use a local neighbourhood, not every sample

Global kriging inverts an n × n matrix. At 500 points that is fine; at 50,000 it is impossible, and it is also unnecessary — samples beyond the range contribute weights near zero.

Take the nearest 16 to 32. The 24 used here gave 40.9 m and ran 12,000 predictions in 1.3 seconds.

3. Understand the Lagrange row

The (k+1) × (k+1) system with a row and column of ones and a zero corner is not padding. It is the constraint that the weights sum to one, which is what makes ordinary kriging unbiased when the mean is unknown.

Drop it and you get simple kriging, which needs the true mean supplied and is biased if you supply the wrong one.

4. Do not count on caching the matrix inverse

The left-hand matrix depends only on the neighbour set, not on the target, so it is tempting to cache it — adjacent cells "obviously" share neighbours.

Measured on the full 30 m grid, 110,889 targets against 500 scattered samples with k=24, that assumption is nearly worthless:

110,889 targets -> 108,650 distinct neighbour sets (2.0% cache hits), 12.0 s

With 24 neighbours drawn from a scattered sample, the set changes almost every cell. Keep the cache — it is three lines and it does help when k is small or the samples are gridded — but the thing that actually makes kriging tractable is the local neighbourhood, not the cache.

5. Read the variance as geometry, not as accuracy

mean kriging variance   4,450 m²  ->  sd 66.7 m
actual RMSE                                40.9 m

Kriging variance depends only on the sample positions and the variogram model. The measured values never enter it. Two datasets with identical sample layouts and completely different values produce identical variance maps.

It is therefore excellent for comparing across a map — this corner is better supported than that one — and unreliable as an absolute interval. Here it was 63% too wide. My kriging fit fails or returns a flat surface traces the same overstatement on these samples to the fitted curve rising too steeply at the shortest lags.

From variogram to raster, in blocks

The data used on this page

The surface is real: the Copernicus DEM GLO-30 over a 10 × 10 km window centred on the summit of Snowdon (Yr Wyddfa), warped to a 30 m grid in EPSG:32630 (x 422,880–432,870, y 5,875,440–5,885,430; 333 × 333 cells, 99.8 km²). The 500 samples are cell centres drawn at random with a fixed seed; 12,000 further cells, none of them samples, are held out as test_points/test_values. grid_points/truth are every cell centre and its elevation, and truth_fn reads the DEM bilinearly at any location, which the error tables score against. idw is the function from Spatial interpolation explained, used here for the comparison table; the kriging is this page's own code.

import os

import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.warp import reproject
from scipy.interpolate import RegularGridInterpolator

bounds = (422_880.0, 5_875_440.0, 432_870.0, 5_885_430.0)
transform = from_origin(bounds[0], bounds[3], 30, 30)
shape = (333, 333)
if not os.path.exists("snowdon_dem_utm30.tif"):
    COP = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{0}_00_DEM/"
           "Copernicus_DSM_COG_10_N53_00_{0}_00_DEM.tif")
    dem = np.full(shape, np.nan, "float32")
    for tile in ("W005", "W004"):
        part = np.full(shape, np.nan, "float32")
        with rasterio.open(COP.format(tile)) as src:
            reproject(rasterio.band(src, 1), part, dst_transform=transform, dst_crs="EPSG:32630",
                      dst_nodata=np.nan, resampling=Resampling.bilinear)
        dem = np.where(np.isnan(dem), part, dem)
    with rasterio.open("snowdon_dem_utm30.tif", "w", driver="GTiff", width=shape[1], height=shape[0], count=1,
                       dtype="float32", crs="EPSG:32630", transform=transform, nodata=np.nan) as dst:
        dst.write(dem, 1)
dem = rasterio.open("snowdon_dem_utm30.tif").read(1)

xs = bounds[0] + (np.arange(shape[1]) + 0.5) * 30
ys = bounds[3] - (np.arange(shape[0]) + 0.5) * 30
gx, gy = np.meshgrid(xs, ys)
grid_points = np.column_stack([gx.ravel(), gy.ravel()])
truth = dem.ravel().astype("float64")

rng = np.random.default_rng(0)
order = rng.permutation(truth.size)
sample_idx, test_idx = order[:500], order[500:12_500]
points, values = grid_points[sample_idx], truth[sample_idx]
test_points, test_values = grid_points[test_idx], truth[test_idx]
_bilinear = RegularGridInterpolator((ys[::-1], xs), dem[::-1], bounds_error=False, fill_value=None)
truth_fn = lambda xy: _bilinear(np.column_stack([xy[:, 1], xy[:, 0]]))
from scipy.spatial import cKDTree


def idw(points, values, targets, power=2, k=12):
    tree = cKDTree(points)
    distance, index = tree.query(targets, k=min(k, len(points)))
    distance, index = distance.reshape(len(targets), -1), index.reshape(len(targets), -1)
    weights = 1.0 / np.maximum(distance, 1e-9) ** power
    predicted = (weights * values[index]).sum(1) / weights.sum(1)
    exact = distance[:, 0] < 1e-9
    predicted[exact] = values[index[exact, 0]]
    return predicted, distance[:, 0]


surface = idw(points, values, grid_points)[0].reshape(shape)
print(f"DEM {np.nanmin(dem):.0f} to {np.nanmax(dem):.0f} m; {len(points)} samples, {len(test_points):,} held out")

DEM 52 to 1072 m; 500 samples, 12,000 held out

Example 1 — the whole pipeline, variogram to raster

import numpy as np
from scipy.optimize import curve_fit
from scipy.spatial import cKDTree


def krige_surface(points, values, targets, bin_width=250, max_lag=None, k=24):
    """Empirical variogram, spherical fit, then ordinary kriging."""
    # ---- empirical variogram -------------------------------------------
    d = np.sqrt(((points[:, None, :] - points[None, :, :]) ** 2).sum(-1))
    gamma = 0.5 * (values[:, None] - values[None, :]) ** 2
    iu = np.triu_indices(len(points), 1)
    d, gamma = d[iu], gamma[iu]

    max_lag = max_lag or float(np.percentile(d, 60))
    edges = np.arange(0, max_lag + bin_width, bin_width)
    lags, semis = [], []
    for lo, hi in zip(edges[:-1], edges[1:]):
        m = (d >= lo) & (d < hi)
        if m.sum() >= 30:
            lags.append((lo + hi) / 2)
            semis.append(gamma[m].mean())
    lags, semis = np.array(lags), np.array(semis)

    # ---- fit -----------------------------------------------------------
    params, _ = curve_fit(
        spherical, lags, semis,
        p0=[semis[0] * 0.1, semis.max(), lags.max() / 2],
        bounds=([0, 0, lags[0]], [semis.max(), semis.max() * 3, lags.max() * 3]),
        maxfev=20000,
    )
    nugget, sill, rng = params
    print(f"  nugget {nugget:,.0f}  sill {sill:,.0f}  range {rng:,.0f} m  (largest lag {lags.max():,.0f} m)")
    if rng > lags.max():
        print("  ! range beyond the largest lag — treat it as a lower bound")

    # ---- krige ---------------------------------------------------------
    predicted, variance = ordinary_kriging(points, values, targets, params, k=k)
    print(f"  predicted {predicted.min():.1f} to {predicted.max():.1f}, "
          f"mean sd {np.sqrt(variance.mean()):.1f}")
    return predicted, variance, params


import time
started = time.perf_counter()
grid_pred, grid_var, params = krige_surface(points, values, grid_points)
grid_seconds = time.perf_counter() - started
print(f"  lowest sample {values.min():.1f} m; {len(grid_points):,} targets in {grid_seconds:.1f} s")

  nugget 0  sill 43,743  range 3,910 m  (largest lag 5,875 m)
  predicted 51.3 to 1036.5, mean sd 66.5
  lowest sample 53.4 m; 110,889 targets in 12.0 s

Note that the surface reaches below the lowest sample (53.4 m). That is not a bug — see the explanation of negative weights below — but it is worth knowing before you assume a kriged surface is bounded by its inputs.

Example 2 — kriging in blocks, so it scales

import numpy as np
from scipy.spatial import cKDTree


def krige_blocks(points, values, targets, params, k=24, block=4096):
    """Krige a large target set in chunks, sharing the neighbour cache."""
    tree = cKDTree(points)
    predicted = np.empty(len(targets))
    variance = np.empty(len(targets))
    cache = {}

    for start in range(0, len(targets), block):
        chunk = targets[start:start + block]
        _, neighbours = tree.query(chunk, k=min(k, len(points)))
        neighbours = neighbours.reshape(len(chunk), -1)

        for j, idx in enumerate(neighbours):
            key = idx.tobytes()
            if key not in cache:
                local = points[idx]
                dd = np.sqrt(((local[:, None, :] - local[None, :, :]) ** 2).sum(-1))
                A = np.ones((len(idx) + 1, len(idx) + 1))
                A[:len(idx), :len(idx)] = spherical(dd, *params)
                A[len(idx), len(idx)] = 0.0
                cache[key] = np.linalg.pinv(A)

            b = np.ones(len(idx) + 1)
            b[:len(idx)] = spherical(
                np.sqrt(((points[idx] - chunk[j]) ** 2).sum(-1)), *params)
            w = cache[key] @ b
            predicted[start + j] = w[:len(idx)] @ values[idx]
            variance[start + j] = w @ b

        if start + block >= len(targets):
            print(f"  {start + len(chunk):,}/{len(targets):,} "
                  f"({len(cache):,} distinct neighbour sets cached)")
    return predicted, variance


_ = krige_blocks(points, values, grid_points, params, block=len(grid_points))

  110,889/110,889 (108,650 distinct neighbour sets cached)

Print the cache size and be honest about what it buys. On the grid measured here it was 2.0% — the neighbour set changed almost every cell. Chunking still matters for memory, and the cache pays off when k is small or the samples are themselves on a grid.

Example 3 — validating the variogram by its own cross-validation

import numpy as np


def krige_diagnostics(points, values, params, k=24):
    """Leave-one-out kriging, checking both the predictions and the variances."""
    n = len(points)
    predicted = np.empty(n)
    variance = np.empty(n)

    for i in range(n):
        keep = np.ones(n, bool); keep[i] = False
        p, v = ordinary_kriging(points[keep], values[keep],
                                points[i:i + 1], params, k=k)
        predicted[i], variance[i] = p[0], v[0]

    residual = values - predicted
    rmse = float(np.sqrt(np.mean(residual ** 2)))

    # standardised residuals should have a variance of about 1 if the
    # variogram is right — this is the real test of the model
    standardised = residual / np.sqrt(np.maximum(variance, 1e-12))
    print(f"  LOO RMSE                    {rmse:8.2f}")
    print(f"  mean standardised residual  {standardised.mean():8.3f}  (want ~0)")
    print(f"  variance of standardised    {standardised.var():8.3f}  (want ~1)")

    if standardised.var() < 0.7:
        print("  ! variances are too large — the sill is probably overestimated")
    elif standardised.var() > 1.4:
        print("  ! variances are too small — the model is overconfident")
    return {"rmse": rmse, "std_var": float(standardised.var())}


_ = krige_diagnostics(points, values, params)

  LOO RMSE                       39.45
  mean standardised residual     0.001  (want ~0)
  variance of standardised       0.291  (want ~1)
  ! variances are too large — the sill is probably overestimated

The variance of the standardised residuals is the diagnostic almost nobody runs, and it is the one that tests the variogram rather than the interpolation. Below 1 means the model's uncertainty is inflated — which is exactly the 63% overstatement measured at the top of this page.

Example 4 — the numbers quoted on this page

from scipy.interpolate import griddata

rmse = lambda p, o: float(np.sqrt(np.nanmean((p - o) ** 2)))
started = time.perf_counter()
kriged, kvar = ordinary_kriging(points, values, test_points, params)
seconds = time.perf_counter() - started
for name, pred in (("kriging", kriged),
                   ("TIN, linear", griddata(points, values, test_points, method="linear")),
                   ("IDW, power 2", idw(points, values, test_points)[0]),
                   ("nearest neighbour", griddata(points, values, test_points, method="nearest"))):
    print(f"{name:18} {rmse(pred, test_values):5.1f} m")
print(f"{len(test_points):,} predictions in {seconds:.1f} s; "
      f"mean kriging variance {kvar.mean():,.0f} m² -> sd {np.sqrt(kvar.mean()):.1f} m")
print(f"data variance {values.var():,.0f} m²; sill / variance {params[1] / values.var():.2f}")
print(f"surface sd {grid_pred.std():.1f} m against the DEM's {truth.std():.1f} m")

tree = cKDTree(points)
_, nbrs = tree.query(grid_points, k=24)
negative = 0
for idx, target in zip(nbrs, grid_points):
    local = points[idx]
    A = np.ones((25, 25)); A[:24, :24] = spherical(np.sqrt(((local[:, None] - local[None]) ** 2).sum(-1)), *params)
    A[24, 24] = 0.0
    b = np.ones(25); b[:24] = spherical(np.sqrt(((local - target) ** 2).sum(-1)), *params)
    negative += bool((np.linalg.solve(A, b)[:24] < 0).any())
print(f"targets with at least one negative weight: {negative:,} of {len(grid_points):,}")
print(f"load average {os.getloadavg()[0]:.1f}")

kriging             40.9 m
TIN, linear         49.1 m
IDW, power 2        56.1 m
nearest neighbour   65.4 m
12,000 predictions in 1.3 s; mean kriging variance 4,450 m² -> sd 66.7 m
data variance 38,896 m²; sill / variance 1.12
surface sd 199.6 m against the DEM's 206.9 m
targets with at least one negative weight: 110,889 of 110,889
load average 0.9

Why kriging beats IDW, and what its variance is not

Why kriging beats IDW

Two mechanisms, both absent from IDW.

Declustering. The left-hand matrix holds the variogram between samples. Two samples close together have low semivariance between them, and the solved weights account for that redundancy — the pair jointly gets roughly the weight one of them would get alone. IDW gives both full weight, so a cluster of samples outvotes an isolated one that carries more information.

Structure-derived decay. The rate at which influence falls off comes from the fitted variogram, not from a chosen exponent. If the field is smooth over kilometres, kriging spreads weight widely; if it decorrelates in 200 m, it concentrates.

Why the variance is not a confidence interval

Look at the kriging system: A holds variogram values between samples, b holds variogram values from samples to the target. Neither contains a single measured value. The variance w·b therefore depends only on where the samples are and what the variogram model says — not on what was measured there.

That has a useful consequence: you can compute the variance map before collecting the data, and use it to design the survey. It has an awkward one: it cannot detect that the model is wrong, and it inherits every error in the variogram.

Here the sill was 12% above the data variance, and the variance scales directly with the sill; most of the rest comes from the model's shape at short lags, which is where each prediction's nearest samples sit. Together they gave a predicted standard deviation of 66.7 m against a real RMSE of 40.9 m.

Why the surface is smoother than the data — but not bounded by it

Kriging minimises expected squared error, and the estimator that minimises squared error is a conditional mean. Means are smoother than the things they average, so a kriged surface systematically under-represents extremes. Measured over the full grid: the kriged surface had a standard deviation of 199.6 m against the DEM's 206.9 m.

What surprises people is that the surface is nevertheless not bounded by the samples. It reached 51.3 m where the lowest sample was 53.4 m.

The reason is that kriging weights can be negative. On the grid measured here, every single target — 110,889 of 110,889 — had at least one negative weight. This is the screen effect: when a sample sits directly behind a nearer one along the same bearing, the system gives the far one a small negative weight, because the near sample already carries its information and the far one now mostly contributes redundancy.

Negative weights are what let kriging extrapolate a local gradient slightly beyond the data. They are also why a kriged concentration or rainfall surface can come out negative, which is physically impossible. If that matters, either constrain the weights to be non-negative, or transform the variable (log, logit) before kriging and back-transform after.

If you need a surface with realistic variability rather than the smoothest one consistent with the data, use conditional simulation.

When kriging is not worth it

• No spatial structure. A flat variogram means kriging degenerates to the mean.
• A poor variogram fit. Garbage in; IDW with sensible parameters beats kriging with a bad model.
• Very large sample sets without a local neighbourhood, where the matrix algebra dominates.
• Categorical variables. Indicator kriging exists but is a different method.
• A strong trend. Ordinary kriging assumes a locally constant mean. With a strong regional trend, remove it first and krige the residuals.

Singular matrices, pinv and negative weights

• Duplicate coordinates make the matrix singular. Aggregate them, or use pinv as above.
• np.linalg.pinv is a safety net, not a fix. If it is doing real work, the neighbourhood has duplicates.
• Project the coordinates. Kriging weights by distance.
• A fitted range beyond your largest lag is a lower bound, not a measurement.
• Standardised residual variance well below 1 means inflated variances — usually an overestimated sill.
• Kriged surfaces are smoother than reality (sd 199.6 against the DEM's 206.9 here). Do not use them where variance matters; simulate instead.
• Kriging weights can be negative — all 110,889 targets had at least one here — so the surface can leave the range of the data and a non-negative variable can come out negative.
• Use 16–32 neighbours. More costs time and changes little.
• The nugget is a discontinuity at zero, which is why spherical() special-cases h == 0. Omit that and the surface no longer honours the data.

Internal links

• The variogram explained — fitting and checking the model kriging needs
• My kriging fit fails or returns a flat surface — when the fit goes wrong
• IDW, kriging, splines or TIN? Choosing an interpolator — whether the 27% is worth it
• Spatial interpolation explained — the error-versus-distance picture
• How to cross-validate an interpolated surface — validating properly
• How to interpolate points to a grid with IDW in Python — the simpler baseline
• Sample design explained: where to measure — designing from the variance map
• Spatial autocorrelation explained — the structure kriging exploits
• How to interpolate air quality station data — when there are too few stations for a variogram

FAQ

How do I do kriging in Python?

Fit a variogram model to the empirical variogram, then solve the ordinary kriging system for each target using its nearest 16–32 samples. It is about forty lines with NumPy and SciPy.

Is kriging more accurate than IDW?

It was 27% better here — 40.9 m against 56.1 m on the same 500 samples and 12,000 held-out cells. The gain depends on having a well-fitted variogram.

What does the kriging variance mean?

How poorly the target is surrounded by samples, given the variogram. It never sees the measured values, so it is a geometry map rather than an error estimate — here it overstated the true error by 63%.

Why is my kriged surface flatter than my data?

Because kriging is a conditional-mean estimator, and means are smooth. The kriged surface here had a standard deviation of 199.6 m against the DEM's 206.9 m. Use conditional simulation if you need realistic variability.

Why did my kriged surface produce a negative value?

Kriging weights can be negative — the screen effect — so the surface is not bounded by the data. Every target on the grid measured here had at least one negative weight. Constrain the weights or transform the variable if negatives are impossible.

How many neighbours should I use?

16 to 32. Global kriging inverts an n × n matrix for no benefit, because samples beyond the range get near-zero weight anyway.

Why is my kriging matrix singular?

Almost always duplicate sample coordinates. Aggregate them before kriging.

How do I know the variogram model is right?

Compute leave-one-out standardised residuals — the residual divided by the kriging standard deviation. Their variance should be about 1. Well below 1 means the variances are inflated.