IDW Interpolation with SciPy and GeoPandas

Inverse distance weighting predicts each unsampled location as $\hat{z}(s_0) = \sum_i w_i z_i / \sum_i w_i$ with $w_i = 1/d_i^{\,p}$. Build a scipy.spatial.cKDTree over your sample coordinates, query the $k$ nearest neighbours for every grid cell, raise their distances to the power $p$, and take the weighted mean. Tune $p$ 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 pillar, and shows how to tune it honestly and where it gives way to kriging.

Environment #

pip install numpy==1.26.4 scipy==1.13.0 geopandas==0.14.4 shapely==2.0.4 matplotlib==3.9.0

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.

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()

Step 2 — Build a KD-Tree for Neighbour Search #

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

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.

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 $k$ nearest samples, weight them by $1/d^{\,p}$, and take the weighted mean. Handle the exact-hit case ($d = 0$) explicitly so a grid point coincident with a sample returns that sample’s value.

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}")

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.

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 #

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.

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 $k$ 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 = 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.

Troubleshooting #

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

• Inverse Distance Weighting — parent guide on distance-decay interpolation
• Step-by-Step Ordinary Kriging with PyKrige — the geostatistical alternative with variance output
• Ordinary & Universal Kriging — when to graduate from IDW to kriging

← Back to Inverse Distance Weighting