Tuning the IDW Power Parameter with Cross-Validation

TL;DR: Do not accept $p = 2$. Wrap IDW in a function of p and k, sweep p over np.arange(0.5, 6.01, 0.25) with a leave-one-out loop built on scipy.spatial.cKDTree.query(xy, k=k+1), and plot RMSE against $p$. When the samples are clustered, replace the leave-one-out loop with sklearn.model_selection.GroupKFold over spatial blocks.

Why This Matters #

The power parameter in inverse distance weighting is a hyperparameter like any other, and it is almost universally left at two. That default comes from Shepard’s 1968 paper, where the squared inverse distance was a convenient choice rather than a derived one. Nothing about soil carbon, rainfall or air quality obliges the weights to decay as $d^{-2}$, and unlike ordinary kriging — which infers its weights from a fitted variogram — IDW gives you no data-driven mechanism for the decay unless you build one. Cross-validation is that mechanism, and it is cheap enough that leaving $p$ untuned is simply a decision not taken.

Two other knobs travel with it. The neighbour count $k$ decides how many samples enter the weighted mean, and the search radius decides which of them are geometrically eligible. All three control the same underlying quantity — how much of the surface at a prediction point is determined by its single nearest sample — so tuning one while the others sit at arbitrary values is how people end up concluding that “$p$ makes no difference”. The workflow below evaluates them jointly, using the same predictor described in IDW interpolation with SciPy and GeoPandas, and validates the choice with the designs set out in cross-validation strategies.

The estimator itself is one line. For a prediction location $\mathbf{s}_0$ with the $k$ nearest samples at distances $d_i$,

$$\hat{z}(\mathbf{s}_0) = \frac{\sum_{i=1}^{k} d_i^{-p} \, z_i}{\sum_{i=1}^{k} d_i^{-p}} ,$$

and the entire behaviour of $p$ follows from the ratio between any two weights, $(d_j / d_i)^{p}$. Raise $p$ and that ratio explodes: a neighbour twice as far away carries a quarter of the weight at $p = 2$ but a sixty-fourth at $p = 6$. In the limit $p \to \infty$ the estimator becomes nearest-neighbour interpolation, producing a piecewise-constant surface of Voronoi tiles — flat plateaux with a bull’s-eye of extreme values sitting on each data point. In the other limit, $p \to 0$, every weight tends to one, the estimator becomes the unweighted mean of the $k$ neighbours, and with $k$ large enough that is just the global mean drawn as a flat plane.

Environment and Version Pinning #

Nothing exotic is needed. The neighbour search comes from SciPy, the fold splitter from scikit-learn, and geopandas only carries the coordinates and the CRS.

pip install "numpy==1.26.4" "scipy==1.13.1" "pandas==2.2.2" \
            "geopandas==1.0.1" "scikit-learn==1.5.1" "shapely==2.0.5"

import numpy as np
import pandas as pd
import geopandas as gpd
from scipy.spatial import cKDTree
from sklearn.model_selection import GroupKFold

Step-by-Step Implementation #

1. Build a clustered survey and measure its spacing #

The dataset is a simulated topsoil organic carbon survey in a 5 km square: six sampling campaigns of 26 points each, plus 64 scattered points. Clustering is deliberate, because it is what real field survey looks like and it is what breaks naive validation later.

rng = np.random.default_rng(11)
SIDE = 5000.0

centres = rng.uniform(500, SIDE - 500, size=(6, 2))
parts = [c + rng.normal(0, 350, size=(26, 2)) for c in centres]
parts.append(rng.uniform(0, SIDE, size=(64, 2)))
xy = np.clip(np.vstack(parts), 0, SIDE)

# Smooth random field via random Fourier features, correlation length ~500 m,
# plus a weak east-west trend and 2 g/kg of measurement noise.
fr = np.random.default_rng(5)
M, L = 300, 500.0
Wf = fr.normal(0, 1.0 / L, size=(M, 2))
phase = fr.uniform(0, 2 * np.pi, M)
amp = fr.normal(0, 1, M)
signal = 40 + np.sqrt(2.0 / M) * (np.cos(xy @ Wf.T + phase) @ amp) * 9.0
z = signal + 0.0015 * xy[:, 0] + np.random.default_rng(3).normal(0, 2.0, len(xy))

gdf = gpd.GeoDataFrame(
    {"soc": z},
    geometry=gpd.points_from_xy(xy[:, 0], xy[:, 1]),
    crs="EPSG:32633",          # projected — distances in metres
)

tree = cKDTree(xy)
nn = tree.query(xy, k=2)[0][:, 1]          # distance to the nearest other sample
print(f"n = {len(z)}, mean = {z.mean():.2f}, sd = {z.std(ddof=1):.2f} g/kg")
print(f"nearest-neighbour distance  p10 = {np.percentile(nn, 10):.1f} m, "
      f"median = {np.median(nn):.1f} m, p90 = {np.percentile(nn, 90):.1f} m")

n = 220, mean = 44.73, sd = 10.06 g/kg
nearest-neighbour distance  p10 = 38.6 m, median = 134.0 m, p90 = 333.9 m

The spread of nearest-neighbour distances — 38.6 m at the tenth percentile against 333.9 m at the ninetieth — is the first warning that leave-one-out will flatter the model.

2. Express IDW as a function of p and k #

Keep the estimator in one small function so that every sweep below calls exactly the same code path as the final surface. Clamping the distance protects against a zero-distance sample, which would otherwise produce inf weights.

def idw_predict(tree, values, targets, p, k, eps=1e-12):
    """Predict at `targets` from `values` held at the points inside `tree`."""
    dist, idx = tree.query(targets, k=k)
    dist = np.maximum(dist, eps)
    w = dist ** (-p)
    return (w * values[idx]).sum(axis=1) / w.sum(axis=1)

3. Sweep the power parameter with leave-one-out #

Leave-one-out for IDW does not need an explicit loop. Query $k+1$ neighbours for every sample and discard the first column, which is the point itself at distance zero; the remaining $k$ columns are exactly the neighbours that would have been used had the point been absent.

def loocv(xy, z, p, k):
    dist, idx = cKDTree(xy).query(xy, k=k + 1)
    dist, idx = dist[:, 1:], idx[:, 1:]        # drop the self-match
    w = np.maximum(dist, 1e-12) ** (-p)
    pred = (w * z[idx]).sum(axis=1) / w.sum(axis=1)
    err = pred - z
    return np.sqrt((err ** 2).mean()), np.abs(err).mean(), err.mean()

rows = [dict(zip(("p", "rmse", "mae", "me"), (p,) + loocv(xy, z, p, k=12)))
        for p in np.arange(0.5, 6.01, 0.5)]
curve = pd.DataFrame(rows)
print(curve.to_string(index=False, float_format=lambda v: f"{v:.4f}"))

     p   rmse    mae      me
0.5000 4.3076 3.2358  0.4638
1.0000 3.9878 2.9551  0.4346
1.5000 3.7708 2.8005  0.4102
2.0000 3.6382 2.7341  0.3884
2.5000 3.5615 2.7073  0.3696
3.0000 3.5255 2.7111  0.3527
3.5000 3.5195 2.7319  0.3372
4.0000 3.5255 2.7627  0.3231
4.5000 3.5602 2.7948  0.3101
5.0000 3.5931 2.8307  0.2983
5.5000 3.6284 2.8638  0.2877
6.0000 3.6633 2.8921  0.2781

The coarse minimum is at $p = 3.5$. Refining on a 0.05 grid between 3 and 4 moves it slightly:

fine = np.arange(3.0, 4.001, 0.05)
scores = np.array([loocv(xy, z, p, k=12)[0] for p in fine])
best_p = float(fine[scores.argmin()])
print(f"refined p = {best_p:.2f}, RMSE = {scores.min():.4f} g/kg")
print(f"p = 2 for comparison: RMSE = {loocv(xy, z, 2.0, 12)[0]:.4f} g/kg")

refined p = 3.35, RMSE = 3.5187 g/kg
p = 2 for comparison: RMSE = 3.6382 g/kg

4. Sweep p and k together #

Tuning $p$ at a fixed $k$ answers only half the question. Because both parameters control how concentrated the weights are, the grid must be two-dimensional.

ks = [4, 8, 12, 24, 48]
ps = [0.5, 1.0, 2.0, 3.0, 4.0, 6.0]
grid = pd.DataFrame(
    [[loocv(xy, z, p, k)[0] for k in ks] for p in ps],
    index=[f"p={p}" for p in ps], columns=[f"k={k}" for k in ks],
)
print(grid.to_string(float_format=lambda v: f"{v:.3f}"))
print("\ncolumn minima:", grid.idxmin().to_dict())

        k=4   k=8  k=12  k=24  k=48
p=0.5 3.575 3.976 4.308 5.206 5.888
p=1.0 3.485 3.736 3.988 4.688 5.164
p=2.0 3.432 3.513 3.638 3.979 4.174
p=3.0 3.473 3.471 3.526 3.663 3.730
p=4.0 3.552 3.513 3.534 3.584 3.600
p=6.0 3.709 3.661 3.663 3.669 3.669

column minima: {'k=4': 'p=2.0', 'k=8': 'p=3.0', 'k=12': 'p=3.0', 'k=24': 'p=4.0', 'k=48': 'p=4.0'}

5. Re-run the sweep with spatially blocked folds #

With six tight campaigns in the sample, leaving out one point leaves its near-duplicates in the training set, and the score measures interpolation between neighbours 40 m apart rather than prediction at an unvisited location. Assign each sample to a 1 km block and hold out whole blocks — the design described in spatial block cross-validation in Python.

BLOCK = 1000.0
groups = (np.floor(xy[:, 0] / BLOCK).astype(int) * 100
          + np.floor(xy[:, 1] / BLOCK).astype(int))

def blocked_rmse(xy, z, p, k, groups, n_splits=5):
    errs = []
    for tr, te in GroupKFold(n_splits=n_splits).split(xy, z, groups=groups):
        pred = idw_predict(cKDTree(xy[tr]), z[tr], xy[te], p, min(k, len(tr)))
        errs.append(pred - z[te])
    errs = np.concatenate(errs)
    return float(np.sqrt((errs ** 2).mean()))

both = pd.DataFrame({
    "p": np.arange(0.5, 6.01, 0.5),
    "loo": [loocv(xy, z, p, 12)[0] for p in np.arange(0.5, 6.01, 0.5)],
    "blocked": [blocked_rmse(xy, z, p, 12, groups) for p in np.arange(0.5, 6.01, 0.5)],
})
print(f"{np.unique(groups).size} blocks, fold sizes "
      f"{[len(te) for _, te in GroupKFold(5).split(xy, z, groups=groups)]}")
print(both.to_string(index=False, float_format=lambda v: f"{v:.3f}"))

26 blocks, fold sizes [44, 44, 44, 44, 44]
  p   loo  blocked
0.5 4.308    7.114
1.0 3.988    6.932
1.5 3.771    6.752
2.0 3.638    6.588
2.5 3.562    6.454
3.0 3.526    6.354
3.5 3.520    6.287
4.0 3.526    6.246
4.5 3.560    6.224
5.0 3.593    6.215
5.5 3.628    6.215
6.0 3.663    6.221

6. Fit the final surface and record the choice #

Interpolate with the selected pair, then check the surface against the sample range. IDW is an exact, bounded interpolator: no predicted value may fall outside the observed minimum and maximum, and a surface that does indicates a bug.

best_p, best_k = 3.35, 12
g = np.arange(25, SIDE, 50)
gx, gy = np.meshgrid(g, g)
targets = np.column_stack([gx.ravel(), gy.ravel()])
surface = idw_predict(cKDTree(xy), z, targets, best_p, best_k).reshape(gx.shape)

print(f"surface {surface.shape}: {surface.min():.2f} to {surface.max():.2f} g/kg")
print(f"samples          : {z.min():.2f} to {z.max():.2f} g/kg")
print(f"selected p = {best_p}, k = {best_k}, blocks = {int(BLOCK)} m, "
      f"loo RMSE = {loocv(xy, z, best_p, best_k)[0]:.3f}, "
      f"blocked RMSE = {blocked_rmse(xy, z, best_p, best_k, groups):.3f}")

surface (100, 100): 16.05 to 63.51 g/kg
samples          : 16.05 to 63.51 g/kg
selected p = 3.35, k = 12, blocks = 1000 m, loo RMSE = 3.519, blocked RMSE = 6.304

Interpreting the Output #

Three numbers carry the result. The RMSE is the headline, on the units of the variable, and it is only meaningful against a baseline: predicting the global mean everywhere scores 10.08 g/kg here, and pure nearest-neighbour interpolation scores 4.10, so a tuned IDW at 3.519 is doing real work but is not far ahead of the crudest possible interpolator. The MAE of 2.725 sits well below the RMSE, which is normal and says the error distribution has a tail of a few badly-predicted points. The mean error of $+0.34$ g/kg is a small positive bias — IDW under-predicts the high values and over-predicts the low ones, as any weighted-average estimator must — and it shrinks steadily as $p$ rises because the estimate leans harder on the single nearest sample.

The interaction table is the part that changes practice. Read down each column and the minimum moves: $p = 2$ at $k = 4$, $p = 3$ at $k = 8$ and $k = 12$, $p = 4$ at $k = 24$ and $k = 48$. That diagonal is a ridge, and every cell on it scores between 3.43 and 3.60 while the off-ridge corners are far worse. The bottom-left corner — few neighbours, high exponent — collapses towards nearest-neighbour interpolation: $p = 6$ with $k = 3$ scores 3.771 against the 1-NN score of 4.10. The top-right corner — many neighbours, low exponent — collapses towards the mean: $p = 0.25$ using all 219 other samples scores 9.50 against the mean baseline of 10.08. Tune along one axis only and you walk off the ridge in one of those two directions.

The blocked comparison is the sharpest lesson on the page. The median distance from a held-out point to its nearest training point rises from 134 m under leave-one-out to 325 m under 1 km blocks, and the RMSE rises with it, from 3.519 to 6.214 g/kg. Anyone reporting the leave-one-out figure as the accuracy of the delivered map is understating the error by nearly half. Note also that the two designs disagree about the exponent — leave-one-out selects 3.35, blocked selects 5.25 — but that disagreement is nearly free: the leave-one-out choice scores 6.304 under blocked validation, 1.5% above the blocked optimum. The 77% gap between the two error estimates matters; the 1.9 difference in $p$ does not.

That is what a flat curve means in practice. Above $p \approx 4$ the blocked scores lie inside half a per cent of each other, so the exponent has stopped being the binding constraint and sample spacing has taken over. A sharp curve carries the opposite message: if RMSE doubled between $p = 2$ and $p = 4$, the weighting scheme would be doing most of the work and would deserve a finer grid and a second validation design. Report the curve so a reader can tell which case they are in.

Critical Best Practices #

Drop the self-match, and check that you dropped exactly one #

The leave-one-out shortcut relies on tree.query(xy, k=k+1) returning the point itself in column zero at distance zero. That holds only if no two samples share coordinates. With duplicated locations — repeat visits, or points snapped to a coarse grid — the self-match may land in column one and the discarded column is a genuine neighbour, leaking the answer into the prediction. Assert np.allclose(dist[:, 0], 0.0) before slicing, and de-duplicate first.

Tune the exponent and the neighbour count on the same grid #

The two parameters trade against each other along the ridge above, so a sequential search that fixes one and optimises the other converges to whichever corner it started nearest. The full grid used here is six exponents by five neighbour counts and takes under a second for 220 points; even at 100,000 points a cKDTree query dominates and can be hoisted out of the loop by querying the largest $k$ once and slicing columns.

Block the folds whenever the samples are clustered #

Leave-one-out on clustered data measures how well a sample predicts its own duplicate. Choose the block size from the sampling geometry: it must comfortably exceed the cluster diameter, or a block boundary will split a campaign and the leakage returns. Here the campaigns have a standard deviation of 350 m, so 1 km blocks contain them; blocks of 250 m would not. Confirm by printing the median distance from held-out points to their nearest training point and checking that it has actually risen.

Judge the winner against baselines, not against the other candidates #

An RMSE curve always has a minimum, even when every candidate is useless. Compute two reference scores under the same folds — the global mean, and $k = 1$ nearest neighbour — and quote them alongside. If the tuned model does not clearly beat both, the honest conclusion is that IDW has little to offer on this data, and the next step is a variogram rather than a finer grid of exponents. That comparison is set out in IDW vs ordinary kriging: which to use.

Do not report the tuning score as the map’s accuracy #

The RMSE at the selected $p$ is optimistically biased, because the same folds chose the parameter. With a curve as flat as this one the bias is small, but the clean fix is a nested design: an inner blocked split to select $(p, k)$ and an outer held-out set, never touched during tuning, to report the number that goes in the metadata.

Troubleshooting #

Next Steps #

Once $p$ and $k$ are chosen and the honest error is recorded, the useful comparison is against a model that infers its weights rather than being tuned into them — see IDW vs ordinary kriging for that decision, and cross-validation strategies for the fold designs that make either comparison fair.

Frequently Asked Questions #

Is p = 2 ever the right default?

It is a reasonable starting point and nothing more. The exponent two entered practice through Shepard’s original 1968 formulation and survives because it is cheap to type, not because any process guarantees it. On the dataset in this guide the leave-one-out minimum sits at $p = 3.35$ with RMSE 3.519, while $p = 2$ gives 3.638, a penalty of 3.4 per cent. That penalty is small here but it is measured rather than assumed, and on sparser or noisier data the same sweep can move the optimum below one.

Should I tune p and k separately or together?

Together, because they trade against each other. The best exponent rises with the neighbour count: on this dataset the leave-one-out optimum is $p = 2.0$ at $k = 4$, $p = 3.0$ at $k = 8$ and $p = 4.0$ at $k = 24$. Tuning $p$ at a fixed large $k$ and then tuning $k$ at that fixed $p$ walks along one axis of a diagonal ridge and lands off it. A two-dimensional grid over roughly a dozen exponents and half a dozen neighbour counts costs seconds and removes the problem.

Why does leave-one-out give a much lower RMSE than blocked cross-validation?

Because leaving out a single point from a clustered sample leaves its near-duplicates in the training set. In this dataset the median distance from a held-out point to its nearest training point is 134 m under leave-one-out and 325 m under 1 km blocks, and RMSE rises from 3.519 to 6.214 accordingly. The blocked number is the honest estimate for a location you have not sampled near. Use leave-one-out only when the samples are genuinely well spread relative to the prediction targets.

What does a flat RMSE curve tell me?

That the choice of exponent barely matters and the error is dominated by something else, usually sample spacing or measurement noise. On the blocked curve here every exponent from 4.0 to 6.0 lands within half a per cent of the minimum, so quoting a selected $p$ of 5.25 to two decimals would be spurious precision. A sharp curve is the opposite signal: it means the weighting is doing real work, and it is worth refining the grid and checking the choice on a second validation design.

Related

• IDW Interpolation with SciPy and GeoPandas — the predictor this page tunes, from k-d tree to raster
• IDW vs Ordinary Kriging: Which to Use — when a fitted variogram beats a tuned exponent
• Spatial Block Cross-Validation in Python — the fold design that keeps clustered samples honest

← Back to Inverse Distance Weighting