ipdw-py

Python implementation of ipdw — Inverse Path Distance Weighting.

Path distances that honour landscape barriers (e.g. land masses in coastal applications) are used as interpolation weights instead of Euclidean distances. This prevents interpolated contours from "bleeding through" land into otherwise unconnected water bodies.

Reference — Stachelek J & Madden CJ (2015). Application of Inverse Path Distance Weighting for high density spatial mapping of coastal water quality patterns. Int J Geographical Information Science. https://doi.org/10.1080/13658816.2015.1018833

---

Algorithm overview

1. Cost raster — a grid where open-water cells carry a low movement cost (default 1) and land cells carry a very high cost (default 10 000).

2. Path distances — for every observation point the accumulated least-cost path distance (via Dijkstra / MCP) to each grid cell is computed with skimage.graph.MCP_Geometric. The raw distances are converted to IDW weights:

   $$w_i(\mathbf{x}) = \left(\frac{r}{d_i(\mathbf{x})}\right)^2 \quad d_i < r, \qquad 0 \text{ otherwise}$$

3. IDW surface — the weighted average

   $$V(\mathbf{x}) = \frac{\sum_i v_i , w_i^p}{\sum_i w_i^p}$$

   is computed at every grid cell.

---

Installation

# editable install from the repo root
pip install -e ".[dev]"

Dependencies

Package	Purpose
numpy	array operations
scipy	KDTree, nearest-neighbour distance estimation
geopandas	point / polygon spatial data
rasterio	raster I/O and georeferencing
scikit-image	MCP_Geometric (Dijkstra on a grid)
shapely	geometry objects
tqdm (optional)	progress bar
matplotlib (optional)	validation plots

---

Quick start

import numpy as np
import geopandas as gpd
from shapely.geometry import Point
from rasterio.transform import from_bounds
from ipdw import costraster_gen, ipdw

# --- 1. Observation points ---
pts = gpd.GeoDataFrame(
    {"salinity": [30.0, 25.0, 28.0, 22.0]},
    geometry=[Point(2e4, 1e4), Point(6e4, 1e4),
              Point(2e4, 5e4), Point(6e4, 5e4)],
    crs="EPSG:32632",
)

# --- 2. Cost raster (all water, 1 km cells) ---
cost = np.ones((60, 80))
transform = from_bounds(0, 0, 80_000, 60_000, 80, 60)

# Add a land barrier
cost[27:33, :56] = 10_000

# --- 3. Interpolate ---
result = ipdw(
    pts, cost, transform,
    range_val=20_000,     # 20 km search radius
    param_col="salinity",
)

Two-step usage (multiple parameters, single path-distance pass)

from ipdw import pathdist_gen, ipdw_interp

weights = pathdist_gen(pts, cost, transform, range_val=20_000)

sal_surface  = ipdw_interp(pts, weights, "salinity",  transform)
temp_surface = ipdw_interp(pts, weights, "temperature", transform)

Validation

from ipdw import error_gen

stats = error_gen(sal_surface, transform, validation_pts, "salinity", plot=True)
print(stats)
# {'r2': 0.87, 'rmse': 1.23, 'pe': 99.1, 'logmse': 0.41, ...}

---

API

costraster_gen(points_gdf, polys_gdf, resolution=None, water_cost=1, land_cost=10000, extent="pnts")

Generate a cost raster from polygon barriers. Returns (cost_array, transform, crs).

pathdist_gen(points_gdf, cost_array, transform, range_val, progress=True)

Compute path-distance weight arrays (shape n_points × rows × cols).

ipdw_interp(points_gdf, weights, param_col, transform, dist_power=1, overlapped=False)

Apply the IDW formula to precomputed weights.

ipdw(points_gdf, cost_array, transform, range_val, param_col, dist_power=1, overlapped=False, progress=True)

High-level one-call interface.

error_gen(result_array, transform, validation_gdf, param_col, plot=False, title="")

Compute R², RMSE, PE and log-MSE.

---

Example

See examples/coastal_example.py for a fully worked synthetic coastal scenario demonstrating the barrier-honouring behaviour of IPDW versus standard Euclidean IDW.

python examples/coastal_example.py