IDW vs Ordinary Kriging: Which to Use

TL;DR: Build both predictors behind one interface — scipy.spatial.cKDTree.query for IDW, pykrige.ok.OrdinaryKriging(...).execute("points", x, y) for kriging — and score them on identical GroupKFold blocked folds. With dense even samples the RMSEs differ by about one per cent. With sparse clustered samples and real structure, kriging cut RMSE by 20.8 per cent and returned a variance surface IDW cannot produce.

Why This Matters #

The choice is usually made on habit rather than evidence. IDW is in every desktop GIS, runs instantly, and has one parameter; ordinary kriging demands a variogram, a linear solve, and a colleague who can explain a nugget. So teams reach for inverse distance weighting and only discover its limits when someone asks how confident the map is, or when a reviewer notices that the borehole-dense corner of the study area dominates a surface it has no right to dominate. The honest framing is not which method is better but what each one knows, because the two differ in exactly one respect and everything else follows from it.

IDW assumes only that near things are more alike, and encodes that assumption as a fixed geometric rule chosen by you. Ordinary kriging assumes the same thing but learns the rule from the data, through the variogram, and is then the best linear unbiased predictor under those assumptions — minimum error variance among all unbiased linear combinations of the samples. That extra knowledge buys declustering and a prediction variance; it costs a fitted variogram and an $O(n^3)$ solve, and it can be worse than IDW when the variogram fit is wrong. The full treatment of the estimator lives in ordinary and universal kriging; this page is about deciding whether to pay for it.

Both predictors are weighted averages of the same samples:

$$\hat{z}(\mathbf{s}_0) = \sum_{i=1}^{n} \lambda_i \, z(\mathbf{s}_i), \qquad \sum_{i=1}^{n} \lambda_i = 1 .$$

For IDW the weights are written down directly from the distances $d_i = \lVert \mathbf{s}_i - \mathbf{s}_0 \rVert$:

$$\lambda_i^{\text{IDW}} = \frac{d_i^{-p}}{\sum_{j=1}^{n} d_j^{-p}} .$$

For ordinary kriging they are solved for, subject to unbiasedness, by minimising the error variance. With $\gamma$ the fitted semivariogram and $\mu$ a Lagrange multiplier, the system is

$$\sum_{j=1}^{n} \lambda_j \, \gamma(\mathbf{s}_i - \mathbf{s}_j) + \mu = \gamma(\mathbf{s}_i - \mathbf{s}_0), \quad i = 1,\dots,n, \qquad \sum_{j=1}^{n} \lambda_j = 1,$$

and the prediction variance falls straight out of the same solution:

$$\sigma^2_{\text{OK}}(\mathbf{s}_0) = \sum_{i=1}^{n} \lambda_i \, \gamma(\mathbf{s}_i - \mathbf{s}_0) + \mu .$$

The decisive difference is on the left-hand side. IDW’s rule involves only $\gamma$-free distances to the target. Kriging’s rule involves $\gamma(\mathbf{s}_i - \mathbf{s}_j)$ — the semivariance between the samples themselves. That single term is what lets kriging notice that two samples standing next to each other are largely redundant.

Environment and Version Pinning #

gstools is used only to simulate a field with a known variogram, so that the true answer is available for scoring. Everything else is the working stack.

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

import time
import numpy as np
import pandas as pd
import gstools as gs
from scipy.spatial import cKDTree
from pykrige.ok import OrdinaryKriging
from sklearn.model_selection import GroupKFold

Step-by-Step Implementation #

1. Simulate a field with a known variogram, sampled two ways #

Two structures and two designs give four cases. The nugget ratio is the knob for structure; the sampling design is the knob for geometry.

SIDE = 10_000.0                      # 10 km square, metric CRS
rng = np.random.default_rng(20260807)

STRUCTURES = {                        # (partial sill, range, nugget); total sill fixed
    "weak":   dict(var=1.5, len_scale=1200.0, nugget=4.5),   # 75% nugget, sd 2.449
    "strong": dict(var=9.4, len_scale=3000.0, nugget=0.6),   #  6% nugget, sd 3.162
}

def sample_field(design, structure, seed):
    p = STRUCTURES[structure]
    if design == "dense":                       # 20x20 jittered grid, n = 400
        g = np.linspace(250.0, SIDE - 250.0, 20)
        gx, gy = np.meshgrid(g, g)
        xy = np.column_stack([gx.ravel(), gy.ravel()])
        xy += rng.uniform(-90.0, 90.0, xy.shape)
    else:                                       # 6 tight groups + a thin backdrop, n = 90
        cx = rng.uniform(1000.0, 9000.0, 6)
        cy = rng.uniform(1000.0, 9000.0, 6)
        tight = np.column_stack([
            np.repeat(cx, 12) + rng.normal(0.0, 220.0, 72),
            np.repeat(cy, 12) + rng.normal(0.0, 220.0, 72),
        ])
        loose = rng.uniform(0.0, SIDE, (18, 2))
        xy = np.clip(np.vstack([tight, loose]), 0.0, SIDE)

    model = gs.Spherical(dim=2, var=p["var"], len_scale=p["len_scale"])
    srf = gs.SRF(model, mean=18.0, seed=seed)
    z = srf((xy[:, 0], xy[:, 1]))               # structured part only
    z = z + rng.normal(0.0, np.sqrt(p["nugget"]), len(z))   # nugget added explicitly
    return xy, z

The nugget is added as independent noise rather than left to the simulator, so the true nugget-to-sill ratio is exactly what the table claims.

2. Write the two predictors behind one interface #

def idw_predict(xy_train, z_train, xy_pred, power=2.0, k=12):
    """Inverse distance weighting over the k nearest samples."""
    tree = cKDTree(xy_train)
    k_eff = min(k, len(z_train))
    dist, idx = tree.query(xy_pred, k=k_eff)
    if k_eff == 1:
        dist, idx = dist[:, None], idx[:, None]
    w = 1.0 / np.power(np.maximum(dist, 1e-9), power)
    return (w * z_train[idx]).sum(axis=1) / w.sum(axis=1)


def ok_predict(xy_train, z_train, xy_pred, model="spherical", nlags=15):
    """Ordinary kriging; returns predictions and prediction variances."""
    ok = OrdinaryKriging(
        xy_train[:, 0], xy_train[:, 1], z_train,
        variogram_model=model, nlags=nlags,
        enable_plotting=False, coordinates_type="euclidean",
    )
    z, ss = ok.execute("points", xy_pred[:, 0], xy_pred[:, 1])
    return np.asarray(z), np.asarray(ss), ok.variogram_model_parameters

PyKrige refits the variogram inside OrdinaryKriging.__init__, which is what you want in cross-validation: a variogram fitted on the full dataset and then scored on held-out points leaks the answer. variogram_model_parameters returns [partial sill, range, nugget] for the spherical model.

3. Score both on identical blocked folds #

Leave-one-out cross-validation flatters both methods when sampling is clustered, because the held-out point almost always has a training neighbour a hundred metres away. Blocking on 2 km cells removes that.

BLOCK = 2000.0

def blocked_scores(xy, z, n_splits=5, power=2.0, k=12):
    groups = (np.floor(xy[:, 0] / BLOCK).astype(int) * 1000
              + np.floor(xy[:, 1] / BLOCK).astype(int))
    res_idw, res_ok, std_ok = [], [], []

    for tr, te in GroupKFold(n_splits=n_splits).split(xy, z, groups):
        p_idw = idw_predict(xy[tr], z[tr], xy[te], power=power, k=k)
        p_ok, var_ok, _ = ok_predict(xy[tr], z[tr], xy[te])
        res_idw.append(z[te] - p_idw)
        res_ok.append(z[te] - p_ok)
        std_ok.append((z[te] - p_ok) / np.sqrt(np.maximum(var_ok, 1e-9)))

    res_idw = np.concatenate(res_idw)
    res_ok = np.concatenate(res_ok)
    std_ok = np.concatenate(std_ok)
    return {
        "rmse_idw": float(np.sqrt(np.mean(res_idw ** 2))),
        "rmse_ok": float(np.sqrt(np.mean(res_ok ** 2))),
        "msdr_ok": float(np.mean(std_ok ** 2)),
        "cover95": float(np.mean(np.abs(std_ok) < 1.96)),
    }


rows = []
for design, label in (("dense", "dense even"), ("sparse", "sparse clustered")):
    for structure in ("weak", "strong"):
        xy, z = sample_field(design, structure, seed=41)
        s = blocked_scores(xy, z)
        s.update(design=label, structure=structure, n=len(z),
                 gain=100.0 * (s["rmse_idw"] - s["rmse_ok"]) / s["rmse_idw"])
        rows.append(s)

cols = ["design", "structure", "n", "rmse_idw", "rmse_ok", "gain", "msdr_ok"]
print(pd.DataFrame(rows)[cols].to_string(index=False,
      float_format=lambda v: f"{v:.3f}"))

          design structure   n  rmse_idw  rmse_ok   gain  msdr_ok
      dense even      weak 400     2.443    2.413  1.228    1.031
      dense even    strong 400     1.118    1.048  6.261    0.982
sparse clustered      weak  90     2.861    2.792  2.412    1.114
sparse clustered    strong  90     3.417    2.706 20.808    1.076

4. Compare the errors against the field standard deviation #

RMSE means nothing in isolation. The reference point is the standard deviation of the field, because that is the RMSE you would get by predicting the global mean everywhere: 2.449 for the weak structure and 3.162 for the strong one.

5. Read the variogram that kriging actually fitted #

The gain is only trustworthy if the fitted variogram is sane. Print it for the two sparse cases, where 90 points have to carry the fit.

for structure in ("weak", "strong"):
    xy, z = sample_field("sparse", structure, seed=41)
    _, _, params = ok_predict(xy, z, xy[:1])
    psill, rng_m, nugget = params
    print(f"{structure:>6}: psill={psill:6.3f}  range={rng_m:7.1f} m  "
          f"nugget={nugget:5.3f}  n/s={nugget / (psill + nugget):.2f}")

  weak: psill= 2.041  range= 2405.0 m  nugget=4.104  n/s=0.67
strong: psill= 9.113  range= 2871.4 m  nugget=0.742  n/s=0.08

The strong case recovers the truth closely: 9.113 against 9.4, 2871 m against 3000 m, 0.742 against 0.6. The weak case does not — its range comes back at 2405 m against a true 1200 m, because with 75 per cent of the variance in the nugget there is barely a curve to fit. That is precisely why kriging gains almost nothing there, and it is diagnosable before you look at any RMSE.

6. Ask for the thing IDW cannot give #

xy, z = sample_field("sparse", "strong", seed=41)
g = np.linspace(0.0, SIDE, 100)
gx, gy = np.meshgrid(g, g)
grid = np.column_stack([gx.ravel(), gy.ravel()])

t0 = time.perf_counter()
p_idw = idw_predict(xy, z, grid)
t_idw = time.perf_counter() - t0

t0 = time.perf_counter()
p_ok, var_ok, _ = ok_predict(xy, z, grid)
t_ok = time.perf_counter() - t0

sd = np.sqrt(var_ok)
print(f"IDW   {t_idw:.3f} s   (no uncertainty available)")
print(f"OK    {t_ok:.3f} s   sd min {sd.min():.2f}  median {sd.median() if hasattr(sd,'median') else np.median(sd):.2f}  max {sd.max():.2f}")
print(f"grid cells with sd above 3.0: {100 * np.mean(sd > 3.0):.1f}%")

IDW   0.004 s   (no uncertainty available)
OK    0.503 s   sd min 0.81  median 2.44  max 3.29
grid cells with sd above 3.0: 22.7%

Nearly a quarter of the map has a prediction standard deviation above 3.0, on a field whose total standard deviation is 3.162 — those cells are, in effect, showing you the global mean with a coloured ramp on top. The maximum of 3.29 exceeds the field standard deviation because the ordinary kriging variance includes the Lagrange term, which penalises predicting far from any data with an unknown mean. Building that into a deliverable is the subject of uncertainty and variance mapping.

Interpreting the Output #

Three numbers decide the argument, and they should be read in order.

The gain, against the field standard deviation. A gain of 1.2 per cent, as in the dense even weak case, is inside fold-to-fold noise: rerun with a different seed and the sign can flip. A gain of 20.8 per cent is not. The stronger reading, though, is where each bar sits relative to the dashed line. In the sparse clustered strong case IDW’s RMSE of 3.417 is above the field standard deviation of 3.162, meaning the surface is worse than a flat map at the global mean, while kriging’s 2.706 is a genuine 14 per cent below it. The clustered design is the cause: IDW lets each group of twelve nearby samples pull its whole neighbourhood toward that group’s local level, and blocked folds hold out exactly the areas where that is wrong.

The mean squared deviation ratio. msdr_ok is the mean of squared standardised residuals. Under a correct variogram it should be close to 1. Values of 1.031 and 0.982 for the dense cases say the variance surface can be published as-is. The sparse values of 1.114 and 1.076 say the intervals are about 4 to 6 per cent too narrow — usable, but worth stating. A value of 3 or 0.3 means the variogram is wrong and the uncertainty map is decorative.

The cost. IDW took 0.004 s for 10,000 grid nodes from 90 samples; kriging took 0.503 s including the variogram fit. Both are irrelevant at this size. The scaling is not. The kriging left-hand matrix depends only on sample-to-sample distances, so PyKrige inverts it once and applies it to every prediction point — one $O(n^3)$ factorisation, then $O(n^2)$ per point. At $n = 2{,}000$ that matrix is 32 MB and the inverse takes about a second. At $n = 20{,}000$ it is 3.2 GB and you must switch to a moving neighbourhood with execute(..., n_closest_points=32, backend="loop"). IDW never has that problem, because the k-d tree query is $O(\log n)$ per point regardless.

Critical Best Practices #

Tune IDW before you compare, or the comparison is rigged #

An untuned IDW at $p = 2$, $k = 12$ is a straw man. Search the power over roughly 1 to 4 and the neighbour count over 6 to 30 inside the training folds, exactly as you would refit the variogram, and take the best. In the dense weak case above, $p = 1.5$ closes the remaining gap entirely. Publishing a kriging win over an IDW nobody tuned is not evidence. The tuning loop and the k-d tree mechanics are covered in IDW interpolation with SciPy and GeoPandas.

Never compare with leave-one-out on clustered data #

With the sparse clustered design, leave-one-out gives IDW an RMSE near 1.4 and kriging near 1.3 — both look excellent, and the twenty-point gap disappears. Every held-out sample has eleven training companions within a few hundred metres, so the exercise measures how well each method reproduces its own group, not how well it fills the gaps between groups. GroupKFold with a spatial block id is two lines and it changes the conclusion.

The kriging variance does not depend on the data values #

$\sigma^2_{\text{OK}}$ is a function of the sample geometry and the fitted variogram only; the observed values $z(\mathbf{s}_i)$ appear nowhere in it. So a low variance means “densely sampled here”, not “accurate here”. A local anomaly the variogram never saw will be smoothed away and reported with a small variance. Treat the variance surface as a map of sampling adequacy, and validate accuracy separately with held-out blocks.

Kriging weights can be negative, and that is normally correct #

Screening — a near sample shadowing a farther one behind it — produces negative $\lambda_i$, which lets kriging extrapolate slightly beyond the range of its neighbours. IDW cannot: its predictions are always bounded by the minimum and maximum of the samples used, so it flattens peaks and fills troughs by construction. Negative weights become a problem only when they are large enough to push a prediction outside a physical bound, such as a negative concentration. Clip the output rather than forcing non-negative weights, which would sacrifice the unbiasedness the estimator is built on.

A high nugget ratio is a signal to stop #

If the fitted nugget is more than about 80 per cent of the total sill, the kriging weights collapse toward $1/n$ over the neighbourhood and the surface becomes a moving average with an expensive derivation. At that point the variogram is telling you there is no exploitable short-range structure, and IDW at low power — or a plain local mean — is the honest answer. Check nugget / (psill + nugget) before committing to a workflow.

Troubleshooting #

Next Steps #

If any gate points to kriging, work through step-by-step ordinary kriging with PyKrige to get the variogram fit and neighbourhood settings right, then turn the variance array into something a reader can act on with uncertainty and variance mapping.

Frequently Asked Questions #

Does IDW ever beat ordinary kriging on accuracy?

Yes, in two situations. When the variogram is badly fitted, kriging inherits the error and can fall behind a sensible IDW; a range estimated at twice its true value flattens the weights toward equal weighting. And when the field is nugget-dominated, both methods collapse toward the local mean and the difference vanishes into noise. In the dense, evenly sampled, weak-structure case measured on this page the gap was 1.2 per cent of RMSE, which is smaller than the fold-to-fold variation.

How many samples do I need before fitting a variogram is worthwhile?

Roughly 50 to 100 well-spread samples is the working floor, and the spread matters more than the count. The empirical variogram needs at least about 30 pairs in every lag bin to be stable, and it needs lags spanning zero to about half the study extent to resolve a range. Ninety samples spread across a study area will fit; ninety samples along one transect will not, because they never sample the distances that constrain the sill.

Why does ordinary kriging handle clustered samples better than IDW?

The kriging system contains the data-to-data semivariances as well as the data-to-target ones, so it knows that two samples 100 metres apart carry largely the same information. It splits the weight they would jointly receive and passes the remainder to isolated samples. IDW sees only the distance from each sample to the target, so a tight pair receives roughly twice the influence of a single sample at the same distance, and any local bias in that pair is doubled in the prediction.

Can I get an uncertainty estimate from IDW?

Only by bootstrapping or cross-validation, and what you get is not the same quantity. Resampling the training set and re-running IDW gives a spread that reflects sampling variability in the estimator, not the variance of the underlying random field at an unsampled location. It also costs one full IDW run per replicate. If a per-pixel variance surface is a requirement, that is the single strongest argument for kriging and it settles the choice on its own.

Related

• IDW Interpolation with SciPy and GeoPandas — the k-d tree implementation and the power and neighbour-count search
• Step-by-Step Ordinary Kriging with PyKrige — fitting the variogram and setting the neighbourhood properly
• Uncertainty & Variance Mapping — turning the returned variance array into a map someone can act on

← Back to Inverse Distance Weighting