My Interpolated Surface Extends Far Beyond the Data
On this page
Problem statement
An interpolated surface fills its whole bounding box. Some of that is interpolation between samples; the rest is extrapolation into places nobody measured, rendered in the same colours with the same apparent authority.
Measured against a real elevation model of Snowdonia as ground truth, with 377 clustered samples interpolated by IDW (Example 1):
distance to nearest sample cells RMSE
0 - 100 m 9,680 21.2 m
100 - 250 m 14,020 50.3 m
250 - 500 m 15,753 101.6 m
500 - 1,000 m 25,193 142.2 m
over 1,000 m 46,243 258.7 m
A twelve-fold range of error across one map, with nothing in the output distinguishing the parts.
The obvious fix — clip to the convex hull of the samples — helps, but bluntly. The hull covered 54.0% of the area, with 123 m of error inside it against 238 m outside. A distance mask keeping the same 54% of the map kept cells with 99.7 m of error — a fifth less — and one keeping only cells within 500 m of a sample, 72 m.
Quick answer
Mask on distance to the nearest sample, not on the hull:
import numpy as np
from scipy.spatial import cKDTree
def mask_unsupported(points, targets, predicted, max_distance):
nearest, _ = cKDTree(points).query(targets, k=1)
masked = np.where(nearest <= max_distance, predicted, np.nan)
print(f" {np.isnan(masked).mean():.1%} masked beyond {max_distance} m")
return masked, nearestA defensible max_distance is the variogram range, or two to three times the mean sample spacing. Both are properties of the data rather than of the plot.
Beyond the hull: mask on distance to a sample
1. Understand why the convex hull is only a rough guide
The hull is the smallest convex polygon containing the samples. With clustered sampling it spans the gaps between clusters, so a point in the middle of a five-kilometre hole is inside the hull.
Measured on the clustered design: the hull covers 54.0% of the area, RMSE 123 m inside and 238 m outside. It separates better from worse — but choosing the same 54% by distance to the nearest sample does better still, at 99.7 m.
The hull answers "is this point surrounded by samples", which is a topological question. The useful question is "is there a sample near this point", which is metric.
2. Mask on distance to the nearest sample
The error bands at the top are ordered by that distance: 21.2 m, 50.3 m, 101.6 m, 142.2 m, 258.7 m. It is the single best predictor of local error available without ground truth.
nearest, _ = cKDTree(points).query(grid, k=1)3. Choose the threshold from the variogram range
Beyond the range, samples carry no information about the target, so the prediction is a local mean dressed as an interpolation.
On the random sample of the same surface the fitted range was 3,762 m (the flat-surface kriging fix fits it), which is generous. Two to three times the mean sample spacing — sqrt(area/n) — is the alternative when no variogram is available; with 500 samples over 99.8 km² that spacing is 447 m.
4. Report what the mask removed
masked 64.4% of the study area beyond 500 m from any sample
That is the clustered design in Example 1. The number is often uncomfortable and always informative. A surface covering 36% of its bounding box is a more honest product than one covering all of it with unmarked extrapolation.
5. Ship the distance raster alongside
If masking is not acceptable — a downstream model needs a complete surface — then ship the distance-to-nearest-sample raster as a second band and say what it means.
Hull and distance masks compared
The data used on this page
The surface is real: the Copernicus DEM GLO-30 over a 10 × 10 km window centred on the summit of Snowdon (Yr Wyddfa), warped to a 30 m grid in EPSG:32630 (x 422,880–432,870, y 5,875,440–5,885,430; 333 × 333 cells, 99.8 km²). The 500 samples are cell centres drawn at random with a fixed seed; 12,000 further cells, none of them samples, are held out as test_points/test_values. grid_points/truth are every cell centre and its elevation, and truth_fn reads the DEM bilinearly at any location, which the error tables score against. idw is the function from Spatial interpolation explained, and surface is its power-2, 12-neighbour surface on that 30 m grid. This page's examples use a clustered design instead of the random one: 8 centres placed at random, 50 samples round each with a 400 m standard deviation (seed 1), snapped to cells — 377 samples once those outside the window are dropped — as cpoints/cvalues.
import os
import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.warp import reproject
from scipy.interpolate import RegularGridInterpolator
bounds = (422_880.0, 5_875_440.0, 432_870.0, 5_885_430.0)
transform = from_origin(bounds[0], bounds[3], 30, 30)
shape = (333, 333)
if not os.path.exists("snowdon_dem_utm30.tif"):
COP = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{0}_00_DEM/"
"Copernicus_DSM_COG_10_N53_00_{0}_00_DEM.tif")
dem = np.full(shape, np.nan, "float32")
for tile in ("W005", "W004"):
part = np.full(shape, np.nan, "float32")
with rasterio.open(COP.format(tile)) as src:
reproject(rasterio.band(src, 1), part, dst_transform=transform, dst_crs="EPSG:32630",
dst_nodata=np.nan, resampling=Resampling.bilinear)
dem = np.where(np.isnan(dem), part, dem)
with rasterio.open("snowdon_dem_utm30.tif", "w", driver="GTiff", width=shape[1], height=shape[0], count=1,
dtype="float32", crs="EPSG:32630", transform=transform, nodata=np.nan) as dst:
dst.write(dem, 1)
dem = rasterio.open("snowdon_dem_utm30.tif").read(1)
xs = bounds[0] + (np.arange(shape[1]) + 0.5) * 30
ys = bounds[3] - (np.arange(shape[0]) + 0.5) * 30
gx, gy = np.meshgrid(xs, ys)
grid_points = np.column_stack([gx.ravel(), gy.ravel()])
truth = dem.ravel().astype("float64")
rng = np.random.default_rng(0)
order = rng.permutation(truth.size)
sample_idx, test_idx = order[:500], order[500:12_500]
points, values = grid_points[sample_idx], truth[sample_idx]
test_points, test_values = grid_points[test_idx], truth[test_idx]
_bilinear = RegularGridInterpolator((ys[::-1], xs), dem[::-1], bounds_error=False, fill_value=None)
truth_fn = lambda xy: _bilinear(np.column_stack([xy[:, 1], xy[:, 0]]))
from scipy.spatial import cKDTree
def idw(points, values, targets, power=2, k=12):
tree = cKDTree(points)
distance, index = tree.query(targets, k=min(k, len(points)))
distance, index = distance.reshape(len(targets), -1), index.reshape(len(targets), -1)
weights = 1.0 / np.maximum(distance, 1e-9) ** power
predicted = (weights * values[index]).sum(1) / weights.sum(1)
exact = distance[:, 0] < 1e-9
predicted[exact] = values[index[exact, 0]]
return predicted, distance[:, 0]
surface = idw(points, values, grid_points)[0].reshape(shape)
crng = np.random.default_rng(1)
centres = np.column_stack([crng.uniform(bounds[0], bounds[2], 8), crng.uniform(bounds[1], bounds[3], 8)])
raw = (centres[:, None, :] + crng.normal(0, 400, (8, 50, 2))).reshape(-1, 2)
raw = raw[(raw[:, 0] > bounds[0]) & (raw[:, 0] < bounds[2]) & (raw[:, 1] > bounds[1]) & (raw[:, 1] < bounds[3])]
cidx = np.unique(((bounds[3] - raw[:, 1]) // 30).astype(int) * shape[1] + ((raw[:, 0] - bounds[0]) // 30).astype(int))
cpoints, cvalues = grid_points[cidx], truth[cidx]
print(f"DEM {np.nanmin(dem):.0f} to {np.nanmax(dem):.0f} m; {len(points)} samples, {len(test_points):,} held out")DEM 52 to 1072 m; 500 samples, 12,000 held out
Example 1 — a masked surface with the error bands reported
import numpy as np
from scipy.spatial import cKDTree
def interpolate_with_support(points, values, grid, predict,
max_distance=None, range_m=None,
bands=((0, 100), (100, 250), (250, 500),
(500, 1000), (1000, np.inf))):
"""Predict, mask beyond the support distance, and report the bands."""
predicted = predict(points, values, grid)
nearest, _ = cKDTree(points).query(grid, k=1)
if max_distance is None:
if range_m is not None:
max_distance = range_m
else:
hull_area = _hull_area(points)
spacing = np.sqrt(hull_area / len(points))
max_distance = 3 * spacing
print(f" no range given; using 3x mean spacing "
f"({spacing:.0f} m) = {max_distance:.0f} m")
print(f" support distance: median {np.median(nearest):.0f} m, "
f"p95 {np.percentile(nearest, 95):.0f} m, "
f"max {nearest.max():.0f} m")
for lo, hi in bands:
band = (nearest >= lo) & (nearest < hi)
if band.sum():
label = f"{lo:>5.0f}-{'inf' if np.isinf(hi) else int(hi):>5}"
print(f" {label} m: {int(band.sum()):8,} cells "
f"({band.mean():6.2%})")
masked = np.where(nearest <= max_distance, predicted, np.nan)
print(f" masked {np.isnan(masked).mean():.1%} beyond {max_distance:.0f} m")
return masked, nearest
def _hull_area(points):
from scipy.spatial import ConvexHull
return float(ConvexHull(points).volume) # 'volume' is area in 2D
predict = lambda p, v, g: idw(p, v, g)[0]
cpredicted = predict(cpoints, cvalues, grid_points)
_ = interpolate_with_support(cpoints, cvalues, grid_points, predict, max_distance=500)
for lo, hi in ((0, 100), (100, 250), (250, 500), (500, 1000), (1000, np.inf)):
near = cKDTree(cpoints).query(grid_points, k=1)[0]
band = (near >= lo) & (near < hi)
print(f" RMSE {lo:>5.0f}-{hi:>6.0f} m: {np.sqrt(np.mean((cpredicted[band] - truth[band]) ** 2)):6.1f} m")
print(f" DEM standard deviation {truth.std():.0f} m") support distance: median 805 m, p95 2710 m, max 3911 m
0- 100 m: 9,680 cells ( 8.73%)
100- 250 m: 14,020 cells (12.64%)
250- 500 m: 15,753 cells (14.21%)
500- 1000 m: 25,193 cells (22.72%)
1000- inf m: 46,243 cells (41.70%)
masked 64.4% beyond 500 m
RMSE 0- 100 m: 21.2 m
RMSE 100- 250 m: 50.3 m
RMSE 250- 500 m: 101.6 m
RMSE 500- 1000 m: 142.2 m
RMSE 1000- inf m: 258.7 m
DEM standard deviation 207 m
Printing the distance percentiles before masking is what makes the threshold a decision rather than a default. A p95 support distance of 2.7 km, as measured on the clustered design, tells you immediately that much of the map is extrapolation.
Example 2 — comparing the hull and distance masks
import numpy as np
from scipy.spatial import Delaunay, cKDTree
def compare_masks(points, grid, truth, predicted, max_distance):
"""Does the hull discriminate error as well as distance does?"""
inside_hull = Delaunay(points).find_simplex(grid) >= 0
nearest, _ = cKDTree(points).query(grid, k=1)
inside_distance = nearest <= max_distance
def rmse(mask):
if not mask.any():
return np.nan
return float(np.sqrt(np.mean((predicted[mask] - truth[mask]) ** 2)))
print(f" convex hull: covers {inside_hull.mean():6.1%}, "
f"RMSE inside {rmse(inside_hull):7.2f}, "
f"outside {rmse(~inside_hull):7.2f}")
print(f" distance <= {max_distance:.0f} m: covers "
f"{inside_distance.mean():6.1%}, "
f"RMSE inside {rmse(inside_distance):7.2f}, "
f"outside {rmse(~inside_distance):7.2f}")
hull_gap = rmse(~inside_hull) - rmse(inside_hull)
distance_gap = rmse(~inside_distance) - rmse(inside_distance)
print(f" discrimination: hull {hull_gap:+.2f}, "
f"distance {distance_gap:+.2f}")
return {"hull_gap": hull_gap, "distance_gap": distance_gap}
_ = compare_masks(cpoints, grid_points, truth, cpredicted, 500)
# the same share of the map, chosen by distance instead of by the hull
inside_hull = Delaunay(cpoints).find_simplex(grid_points) >= 0
nearest, _ = cKDTree(cpoints).query(grid_points, k=1)
equal = nearest <= np.quantile(nearest, inside_hull.mean())
err = lambda m: float(np.sqrt(np.mean((cpredicted[m] - truth[m]) ** 2)))
print(f" keeping {inside_hull.mean():.1%} of the map: hull RMSE {err(inside_hull):.1f}, "
f"nearest-distance RMSE {err(equal):.1f} (threshold {np.quantile(nearest, inside_hull.mean()):.0f} m)") convex hull: covers 54.0%, RMSE inside 123.34, outside 238.15
distance <= 500 m: covers 35.6%, RMSE inside 71.66, outside 224.61
discrimination: hull +114.81, distance +152.95
keeping 54.0% of the map: hull RMSE 123.3, nearest-distance RMSE 99.7 (threshold 895 m)
The discrimination is the point. The hull separates 123 from 238; the 500 m distance mask separates 72 from 225.
Example 3 — writing the surface and its support together
import numpy as np
import rasterio
def write_surface(surface, support, profile, out_path, max_distance,
method, n_samples):
"""Two bands and the provenance, so the product explains itself."""
profile = dict(profile) | {"count": 2, "dtype": "float32",
"nodata": np.nan, "compress": "deflate",
"tiled": True}
with rasterio.open(out_path, "w", **profile) as dst:
dst.write(surface.astype("float32"), 1)
dst.write(support.astype("float32"), 2)
dst.set_band_description(1, "interpolated value")
dst.set_band_description(2, "distance to nearest sample (m)")
dst.update_tags(
method=method, n_samples=str(n_samples),
max_distance_m=str(max_distance),
note="band 1 is masked beyond max_distance_m; error grows "
"steeply with band 2",
)
dst.build_overviews([2, 4, 8])
print(f" wrote {out_path}: {np.isnan(surface).mean():.1%} masked")
return out_pathThe note in the tags is the sentence that would otherwise be lost. Band descriptions survive most processing; a README does not.
Why error rises steeply with distance
Why the convex hull is a blunt boundary
The hull answers a topological question: is this point in the convex closure of the samples. Under clustered sampling, a point in the middle of a large empty region between clusters satisfies that.
Measured, the hull covered 54.0% of the study area, with 123 m of error inside it against 238 m outside — it does separate, because on average the hull's cells are nearer the samples. But it keeps the cells between clusters too.
Distance to the nearest sample answers the metric question directly. Keeping the same share of the map by distance gave 99.7 m of error in the kept cells, a fifth less than the hull, and keeping cells within 500 m gave 72 m.
The hull remains useful for one thing: methods such as linear TIN interpolation genuinely cannot produce values outside it, so it describes their natural extent. It is not a statement about reliability.
Why error grows so steeply with distance
Interpolation weights nearby samples heavily. Within one sample spacing, the prediction is dominated by observations and the error is close to the measurement error.
Beyond the variogram range, the nearest sample carries no information, so the prediction converges on a local mean and the error converges on the standard deviation of the field.
Between those, the error rises smoothly. The measured bands — 21.2, 50.3, 101.6, 142.2, 258.7 m — trace that curve. The last band's 259 m is worse than the DEM's own standard deviation of 207 m, which is what predicting the overall mean everywhere would achieve: far from every sample, IDW predicts the elevation of the nearest cluster, and on a mountain that is a worse guess than the mean.
Why to mask rather than to smooth
The temptation is to leave the surface complete and rely on a caption. Captions do not survive: a raster gets clipped, reprojected, resampled and put in a report, and the caption stays behind.
A NaN survives all of that. It is visible in every viewer, it propagates through every calculation, and it forces a downstream user to make an explicit decision.
Where a complete surface is genuinely required, the distance band is the compromise — it travels with the data and can be applied later.
Why the threshold should come from the data
An arbitrary threshold is a different arbitrary decision from not masking at all.
The variogram range is the principled choice: it is the distance beyond which samples are uninformative, estimated from the data. Where a variogram is unavailable or unreliable, two to three times the mean sample spacing is a reasonable proxy.
Either way, record it. A masked surface without its threshold cannot be compared with another masked surface.
Kriging-variance masks and shipping the distance raster
- The convex hull is a blunt support boundary. For the same share of the map, a distance mask kept cells a fifth more accurate.
- Distance to the nearest sample orders local error better than anything else available without ground truth.
- Set the threshold from the variogram range or from the sample spacing.
- Linear TIN interpolation is naturally hull-bounded; IDW and kriging are not.
- Kriging variance is an alternative mask; check it against held-out error before relying on it — see How to fit a variogram and krige a surface.
- Report what the mask removed, as a percentage of the study area.
- Ship the distance raster if masking is not acceptable.
- Write the threshold into the file tags, not only the documentation.
Internal links
- Spatial interpolation explained — the error-versus-distance measurements
- The variogram explained — where the range comes from
- How to interpolate points to a grid with IDW in Python — where the mask is applied
- How to fit a variogram and krige a surface in Python — kriging variance as an alternative
- Sample design explained: where to measure — avoiding the gaps
- How to cross-validate an interpolated surface — why validation misses this
- Extrapolation in space explained: the area of applicability — the same idea for machine learning
- How to build a TIN and interpolate elevations — a method bounded by the hull
FAQ
Should I clip my interpolated surface to the convex hull?
Only as a rough one. On a clustered sample the hull covered 54% of the area with 123 m of error inside and 238 m outside; a distance mask keeping the same share did a fifth better.
What should I mask on instead?
Distance to the nearest sample. On the same data, cells within 500 m had 72 m of error against 225 m for the rest.
What threshold should I use?
The variogram range, or two to three times the mean sample spacing when no variogram is available.
How much error is there far from the samples?
Measured: 21.2 m within 100 m of a sample, rising to 258.7 m beyond a kilometre — a twelve-fold range across one map.
Can I use the kriging variance as a mask?
It is a reasonable relative map of support; check its calibration against held-out error before using it as a threshold — see How to fit a variogram and krige a surface.
What if I need a complete surface?
Ship the distance-to-nearest-sample raster as a second band, with the interpretation in the file tags, so a downstream user can apply the mask.
Why not just say it in the caption?
Captions do not survive clipping, reprojection and reuse. A NaN does.