Geostatistical Modeling Spatial Statistics with Python

IDW Interpolation with SciPy and GeoPandas

Inverse distance weighting predicts each unsampled location as z^(s0)=∑iwizi/∑iwi\hat{z}(s_0) = \sum_i w_i z_i / \sum_i w_i with wi=1/di pw_i = 1/d_i^{\,p}. Build a scipy.spatial.cKDTree over your sample coordinates, query the kk nearest neighbours for every grid cell, raise their distances to the power pp, and take the weighted mean. Tune pp by leave-one-out RMSE. Coordinates must be in a projected (metric) CRS.

Why This Matters

Inverse distance weighting (IDW) is the fastest, most assumption-light way to turn scattered point measurements into a continuous surface. It needs no variogram, no covariance model, and no distributional assumptions — just distances and one tunable power. That makes it the pragmatic first pass for rainfall, elevation, sensor, and soil data before committing to a full geostatistical model. This guide builds IDW from scratch inside the Inverse Distance Weighting workflow, part of the wider Kriging, Interpolation & Surface Generation Techniques section, and shows how to tune it honestly and where it gives way to kriging.

Environment

bash
pip install numpy==1.26.4 scipy==1.13.0 geopandas==0.14.4 shapely==2.0.4 matplotlib==3.9.0
python
import numpy as np                       # 1.26.4
from scipy.spatial import cKDTree        # scipy 1.13.0
import geopandas as gpd                  # 0.14.4
from shapely.geometry import Point       # shapely 2.0.4
import matplotlib.pyplot as plt          # 3.9.0

Python 3.10 or 3.11 recommended. cKDTree ships with SciPy — no extra install.

Step-by-Step Implementation

Step 1 — Load Points in a Metric CRS

IDW weights come from Euclidean distances, so reproject to a projected CRS before anything else. Here we build a synthetic sample; in practice, read your file and drop the simulation.

python
rng = np.random.default_rng(7)
n = 150

# Synthetic field: smooth trend + local bump, in metres
x = rng.uniform(0, 10000, n)
y = rng.uniform(0, 10000, n)
z = (0.0008 * x + 0.0005 * y
     + 4 * np.exp(-((x - 6000) ** 2 + (y - 4000) ** 2) / (2 * 1800 ** 2))
     + rng.normal(0, 0.15, n))

gdf = gpd.GeoDataFrame(
    {"z": z},
    geometry=[Point(xi, yi) for xi, yi in zip(x, y)],
    crs="EPSG:32632",          # UTM zone 32N — metric
)
coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])
values = gdf["z"].to_numpy()

A cKDTree makes nearest-neighbour queries O(log⁡n)\mathcal{O}(\log n) instead of scanning all points per prediction. This is what lets IDW scale to large grids.

python
tree = cKDTree(coords)

Step 3 — Construct the Prediction Grid with GeoPandas

Span the sample extent with a regular grid. GeoPandas gives a clean way to hold the grid as geometries for later clipping or export.

python
res = 100                      # grid cell size in metres
xmin, ymin, xmax, ymax = gdf.total_bounds
gx = np.arange(xmin, xmax + res, res)
gy = np.arange(ymin, ymax + res, res)
GX, GY = np.meshgrid(gx, gy)
grid_pts = np.column_stack([GX.ravel(), GY.ravel()])

grid_gdf = gpd.GeoDataFrame(
    geometry=[Point(px, py) for px, py in grid_pts], crs=gdf.crs
)
print(f"Grid: {GX.shape} = {grid_pts.shape[0]} prediction cells")

Step 4 — Compute Inverse Distance Weights

For every grid cell, query the kk nearest samples, weight them by 1/d p1/d^{\,p}, and take the weighted mean. Handle the exact-hit case (d=0d = 0) explicitly so a grid point coincident with a sample returns that sample’s value.

python
def idw(tree, values, query_pts, k=12, power=2.0, eps=1e-12):
    """Inverse distance weighting over the k nearest samples.

    tree      : cKDTree built on sample coordinates
    values    : sample values, aligned with tree data order
    query_pts : (m, 2) array of prediction coordinates
    k         : number of nearest neighbours to use
    power     : distance-decay exponent p
    Returns an (m,) array of interpolated values.
    """
    dist, idx = tree.query(query_pts, k=k)
    if k == 1:                       # keep 2D shape for uniform handling
        dist, idx = dist[:, None], idx[:, None]

    # Exact hits: distance 0 would divide by zero — snap to that sample
    exact = dist[:, 0] < eps
    weights = 1.0 / np.power(dist, power)
    weights[exact] = 0.0
    weights[exact, 0] = 1.0

    weights /= weights.sum(axis=1, keepdims=True)
    return np.sum(weights * values[idx], axis=1)

z_grid = idw(tree, values, grid_pts, k=12, power=2.0).reshape(GX.shape)
print(f"Predicted surface: min={z_grid.min():.2f}, max={z_grid.max():.2f}")
How the IDW power redistributes weight among neighbours On the left, one prediction cell sits at the centre of a dashed search circle with five sample points drawn to scale at 120, 210, 340, 480 and 650 metres. On the right, three groups of five bars give the normalised weight each of those samples receives: at power 1 the shares are 42, 24, 15, 11 and 8 percent; at power 2 they are 65, 21, 8, 4 and 2 percent; at power 4 they are 89, 9, 1 percent and under 1 percent for the last two. The power p decides how much weight the nearest sample takes s0 (grid cell) 120 m 210 m 340 m 480 m 650 m k = 5 neighbours shown to scale; the code queries k = 12 share of total weight 120 m 210 m 340 m 480 m 650 m 42% 24% 15% 11% 8% 65% 21% 8% 4% 2% 89% 9% 1% <1% <1% p = 1 p = 2 p = 4 smoother, regional the common default nearly nearest-neighbour

Step 5 — Tune the Power Parameter by Leave-One-Out

The power controls how fast influence decays. Rather than guess, search for the power that minimises leave-one-out (LOO) RMSE. For each held-out sample, interpolate from the others.

python
def loo_rmse(coords, values, k=12, power=2.0):
    tree_full = cKDTree(coords)
    # Query k+1 so we can drop the point itself (its own nearest neighbour)
    dist, idx = tree_full.query(coords, k=k + 1)
    dist, idx = dist[:, 1:], idx[:, 1:]        # exclude self
    w = 1.0 / np.power(dist, power)
    w /= w.sum(axis=1, keepdims=True)
    pred = np.sum(w * values[idx], axis=1)
    return np.sqrt(np.mean((pred - values) ** 2))

powers = np.arange(0.5, 4.01, 0.5)
scores = [(p, loo_rmse(coords, values, k=12, power=p)) for p in powers]
best_power, best_rmse = min(scores, key=lambda t: t[1])

for p, r in scores:
    flag = "  <-- best" if p == best_power else ""
    print(f"power={p:.1f}  LOO RMSE={r:.4f}{flag}")
print(f"\nSelected power={best_power:.1f}  (RMSE={best_rmse:.4f})")

# Regenerate the surface with the tuned power
z_grid = idw(tree, values, grid_pts, k=12, power=best_power).reshape(GX.shape)

Step 6 — Visualise the Surface

python
fig, ax = plt.subplots(figsize=(7, 6))
im = ax.pcolormesh(gx, gy, z_grid, cmap="viridis", shading="auto")
ax.scatter(coords[:, 0], coords[:, 1], c=values, edgecolors="white",
           linewidth=0.4, s=25, cmap="viridis", zorder=5)
fig.colorbar(im, ax=ax, label="Interpolated value")
ax.set_title(f"IDW surface (power={best_power:.1f})")
ax.set_xlabel("Easting (m)"); ax.set_ylabel("Northing (m)")
plt.tight_layout(); plt.savefig("idw_surface.png", dpi=150); plt.show()

Interpreting the Output

  • LOO RMSE curve: the headline diagnostic. It is typically U-shaped — too low a power over-smooths (regional blur), too high a power makes the surface blocky and nearest-neighbour-like. The minimum marks the power best matched to your sample density and field roughness.
  • Best power: often lands near 2 for moderately smooth fields; expect lower values for sparse, gently varying data and higher values for dense, sharply varying data.
  • Surface range: IDW is an exact interpolator that never predicts outside the min/max of the sampled values — it cannot extrapolate peaks or troughs beyond what was observed. Flat plateaus around isolated highs are the classic “bull’s-eye” artefact.
  • No uncertainty: IDW returns only predictions. It gives no variance, standard error, or confidence interval — a fundamental difference from kriging.
The U-shaped leave-one-out RMSE curve and what each end looks like The left chart plots leave-one-out RMSE against the power p over the eight-value search grid from 0.5 to 4.0. The curve falls from a high value at p equals 0.5, reaches its minimum at p equals 2, then climbs again to p equals 4. The right column shows three transect thumbnails: at p equals 1 the prediction drifts to the local mean between samples, at p equals 2 it follows the dashed true field closely, and at p equals 4 it becomes a staircase of flat plateaus. The LOO RMSE curve is U-shaped — each end fails differently IDW prediction true field LOO RMSE (lower is better) 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 power p (search grid 0.5 → 4.0 in steps of 0.5) minimum → the tuned power p ≤ 1: weights nearly uniform → regional blur p ≥ 3: nearest sample dominates → blocky bull’s-eye plateaus p = 1 — drifts to the local mean between samples p = 2 — tracks the field; this is the LOO minimum p = 4 — flat plateaus, a nearest-neighbour map

Critical Best Practices

Always Work in a Projected CRS

Distances drive every weight. In geographic coordinates a degree of longitude shortens toward the poles, distorting distances and corrupting weights. Reproject with gdf.to_crs("EPSG:32632") (or your local UTM zone) before building the KD-tree.

Cap the Neighbourhood with k

Using all points for every prediction is slow and lets distant samples leak in. Restricting to the kk nearest (via tree.query(..., k=k)) keeps IDW local and fast. Tune k alongside the power; 8-16 neighbours is a sensible starting range.

Handle Exact Hits and Duplicates

A grid point coincident with a sample gives d=0d = 0 and a division by zero. The idw function above snaps such cells to the sample value. Deduplicate coincident samples beforehand, or they will each claim full weight.

Tune, Do Not Assume, the Power

Power 2 is a convention, not a law. The leave-one-out search in Step 5 costs almost nothing and routinely finds a better value. Report the tuned power and its RMSE so the surface is reproducible.

Contrast with Kriging

IDW’s simplicity is also its ceiling. Its weights ignore the actual spatial autocorrelation of the data and it produces no uncertainty estimate. When you need statistically optimal weights derived from a fitted variogram, and a prediction variance to map uncertainty, step up to the geostatistical approach in Ordinary & Universal Kriging. A practical rule: prototype with IDW, ship with kriging when uncertainty matters.

What kriging returns that IDW cannot Two side-by-side transects run through the same five sample points. On the left, the IDW prediction is a single line with dashed horizontal guides marking the observed maximum and minimum that it never crosses, and no uncertainty is drawn. On the right, the ordinary-kriging prediction carries a shaded plus-or-minus one standard error band that closes to zero width at every sample and swells between samples and past the ends of the transect. IDW — deterministic w = 1 / d^p — distance and a fixed p, nothing else no band drawn: there is nothing to draw observed max observed min distance along a transect Output: the prediction grid, and nothing else No variance, no standard error, no interval Never leaves the observed min/max envelope One knob: the power p, tuned by LOO RMSE Ordinary kriging — geostatistical weights solved from the fitted variogram γ(h) ±1 kriging standard error distance along a transect Output: prediction plus a kriging variance map Variance is zero at data and grows between Weights follow the modelled autocorrelation Knobs: model, nugget, sill and range

Troubleshooting

Symptom Likely cause Fix
Bull’s-eye artefacts around points Power too high; too few neighbours Lower the power; increase k
Surface over-smoothed / featureless Power too low; too many neighbours Raise the power; reduce k
RuntimeWarning: divide by zero Grid point coincides with a sample (d=0d=0) Use the exact-hit branch shown; deduplicate samples
Distances look wrong / weights odd Data in geographic CRS (degrees) Reproject to a metric CRS before querying
Predictions never reach observed peaks Expected — IDW cannot extrapolate beyond sample range Accept it, or use kriging with an appropriate variogram
Slow on large grids Querying all points per cell Cap neighbours with k; reuse a single cKDTree

Next Steps

Return to the Inverse Distance Weighting guide for weighting variants, then compare against a full geostatistical model with Step-by-Step Ordinary Kriging with PyKrige.

Frequently Asked Questions

What power parameter should I use for IDW?

A power of 2 is the common default, giving a smooth surface where influence decays with the square of distance. Higher powers (3-4) make the surface more local and blocky, approaching a nearest-neighbour map; lower powers (1) produce smoother, more regional surfaces. Do not guess — tune the power by minimising leave-one-out RMSE on your own data, since the best value depends on sample density and the roughness of the field.

How is IDW different from kriging?

IDW is a deterministic interpolator: weights depend only on distance and a fixed power, with no statistical model and no uncertainty estimate. Kriging is geostatistical: it fits a variogram to the data’s spatial autocorrelation, derives optimal weights, and returns a prediction variance. IDW is fast and assumption-light; kriging is more accurate when spatial structure is well modelled and, crucially, quantifies uncertainty.

Why must IDW coordinates be in a projected CRS?

IDW weights are computed from Euclidean distances. In geographic coordinates (degrees), a degree of longitude shrinks toward the poles, so distances are distorted and weights become physically meaningless. Reproject to a metric CRS such as UTM with GeoPandas before building the KD-tree and computing distances.

Related

← Back to Inverse Distance Weighting

Related

IDW vs Ordinary Kriging: Which to Use Tuning the IDW Power Parameter with Cross-Validation