My Interpolated Surface Extends Far Beyond the Data
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 as ground truth, with 400 clustered samples:
distance to nearest sample cells RMSE
0 - 100 m 6,959 18.9 m
100 - 250 m 7,286 54.0 m
250 - 500 m 9,479 106.1 m
500 - 1,000 m 19,802 169.8 m
over 1,000 m 72,611 334.3 m
An eighteen-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 — barely helps. That hull covered 59.0% of the area, and the error inside it was 268 m against 286 m outside.
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.
Step-by-step solution
1. Understand why the convex hull fails
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 a clustered design: hull covers 59.0% of the area, RMSE 268 m inside and 286 m outside. Essentially no discrimination.
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 almost perfectly ordered by that distance: 18.9 m, 54.0 m, 106.1 m, 169.8 m, 334.3 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 measured data the fitted range was 4,537 m, 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 105 km² that spacing is 457 m.
4. Report what the mask removed
masked 61.2% of the study area beyond 500 m from any sample
That number is often uncomfortable and always informative. A surface covering 39% 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.
Code examples
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 2DPrinting the distance percentiles before masking is what makes the threshold a decision rather than a default. A p95 support distance of 3.6 km, as measured on the clustered design, tells you immediately that most 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} convex hull: covers 59.0%, RMSE inside 268.26, outside 285.68
distance <= 500 m: covers 38.8%, RMSE inside 72.14, outside 312.06
The discrimination is the point. The hull separates 268 from 286 — nothing. The distance mask separates 72 from 312.
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.
Explanation
Why the convex hull is the wrong 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.
The measured result is decisive. The hull covered 59.0% of the study area, and the error inside it was 268 m against 286 m outside — a discrimination of 6%.
Distance to the nearest sample answers the metric question, and it separated 72 m from 312 m on the same data.
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 — 18.9, 54.0, 106.1, 169.8, 334.3 m — trace that curve, and the last band's 334 m is approaching the DEM's own standard deviation of 205 m, which is what an uninformed prediction would achieve.
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.
Edge cases or notes
- The convex hull is not a support boundary. It discriminated error by 6% here.
- Distance to the nearest sample is the best available predictor of local error.
- 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, but it overstated error by 68% in one measurement.
- 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?
Not as a reliability mask. On a clustered sample the hull covered 59% of the area with an error of 268 m inside and 286 m outside — essentially no discrimination.
What should I mask on instead?
Distance to the nearest sample. On the same data it separated 72 m inside from 312 m outside.
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: 18.9 m within 100 m of a sample, rising to 334.3 m beyond a kilometre — an eighteen-fold range across one map.
Can I use the kriging variance as a mask?
It is a reasonable relative map of support, but it is not calibrated — it overstated the true error by 68% in one measurement.
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.