How to Build a TIN and Interpolate Elevations
On this page
Problem statement
A triangulated irregular network connects sample points into triangles and interpolates linearly within each. It has one property no other common interpolator has: it refuses to extrapolate.
Measured against a real elevation model as ground truth:
samples TIN coverage TIN RMSE IDW RMSE IDW on TIN's cells
100 94.5% 106.28 m 101.77 m 102.01 m
250 94.3% 69.70 m 79.10 m 73.92 m
500 97.2% 47.56 m 54.92 m 54.91 m
1,000 99.0% 35.06 m 39.60 m 39.60 m
2,000 99.4% 24.62 m 28.71 m 28.70 m
Everything outside the convex hull of the samples is NaN, and that is the point. It also beats IDW at every sample size from 250 up — and still does when both are scored only on the cells TIN answers, so the win is not just the easier question (Example 4).
Quick answer
import numpy as np
from scipy.interpolate import griddata
targets = grid_points
surface = griddata(points, values, targets, method="linear")
coverage = np.isfinite(surface).mean()
print(f"{coverage:.1%} of targets are inside the convex hull")97.7% of targets are inside the convex hull
method="linear" is barycentric interpolation over a Delaunay triangulation. method="cubic" uses a smoother scheme and can overshoot; method="nearest" fills everything and is a different method entirely.
Triangulate, interpolate barycentrically, add breaklines
1. Build the triangulation
from scipy.spatial import Delaunay
tri = Delaunay(points)
print(f"{len(points)} points -> {len(tri.simplices)} triangles")500 points -> 983 triangles
Delaunay triangulation maximises the minimum angle across all triangles, which avoids long thin slivers where interpolation is numerically unstable.
It is unique except when four or more points are co-circular — a regular grid, for instance, where the diagonal of each square can go either way.
2. Interpolate barycentrically within each triangle
A target point inside a triangle is a weighted average of its three vertices, with weights being the barycentric coordinates. All three weights are non-negative and sum to one, which has two consequences:
- the surface is bounded by the sample values — it cannot overshoot
- the surface is exact at every sample and continuous across triangle edges
The gradient is not continuous across edges, which is why a TIN hillshade shows facets.
3. Accept the hull limitation, or fill deliberately
100 samples: 94.5% coverage
500 samples: 97.2%
2,000 samples: 99.4%
The uncovered fraction is the area outside the convex hull. With clustered samples it is far larger.
Filling it means extrapolating, which is what the TIN was declining to do. If you fill, use a different method there and record which cells came from which.
4. Add breaklines where the surface is discontinuous
A cliff, a road cutting or a river bank is a genuine discontinuity. A Delaunay triangulation ignores it and interpolates straight across, smoothing the feature away.
A constrained triangulation forces edges along the breakline so no triangle spans it. scipy does not provide constrained triangulation; triangle and CGAL bindings do.
Without it, a TIN over a levee shows a ramp instead of a wall.
5. Compare against the alternatives on your own data
The measurements above show TIN beating IDW at every size from 250 samples up, but the first comparison is not like for like — TIN declined 2.8% of the cells at 500 samples, and those are the cells furthest from data.
Score both on the cells where both produce a value. Here that barely moved IDW (54.91 m against 54.92 m at 500 samples), and TIN kept its lead.
A TIN with coverage, slivers and long triangles
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, used for the IDW column in Example 4.
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)
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 TIN surface with coverage reported
import numpy as np
from scipy.interpolate import LinearNDInterpolator, NearestNDInterpolator
from scipy.spatial import Delaunay
def tin_surface(points, values, targets, fill=None):
"""Linear TIN interpolation, with the hull coverage reported."""
triangulation = Delaunay(points)
interpolator = LinearNDInterpolator(triangulation, values)
surface = interpolator(targets)
inside = np.isfinite(surface)
print(f" {len(points):,} points -> {len(triangulation.simplices):,} "
"triangles")
print(f" {inside.mean():.1%} of targets inside the convex hull")
print(f" surface {np.nanmin(surface):.2f}..{np.nanmax(surface):.2f}, "
f"data {values.min():.2f}..{values.max():.2f}")
if fill == "nearest" and (~inside).any():
nearest = NearestNDInterpolator(points, values)
surface = np.where(inside, surface, nearest(targets))
print(f" filled {int((~inside).sum()):,} cells with nearest "
"neighbour — record these as extrapolated")
return surface, inside
_ = tin_surface(points, values, grid_points, fill="nearest") 500 points -> 983 triangles
97.7% of targets inside the convex hull
surface 53.45..1036.54, data 53.45..1036.54
filled 2,537 cells with nearest neighbour — record these as extrapolated
Printing the surface range next to the data range confirms the bounding property. A linear TIN cannot exceed its inputs; if it appears to, the interpolation is not linear.
Example 2 — finding the sliver triangles
import numpy as np
from scipy.spatial import Delaunay
def triangle_quality(points, min_angle_deg=15.0):
"""Long thin triangles are where interpolation is least reliable."""
triangulation = Delaunay(points)
vertices = points[triangulation.simplices]
a = np.linalg.norm(vertices[:, 1] - vertices[:, 0], axis=1)
b = np.linalg.norm(vertices[:, 2] - vertices[:, 1], axis=1)
c = np.linalg.norm(vertices[:, 0] - vertices[:, 2], axis=1)
with np.errstate(invalid="ignore"):
angles = np.degrees(np.arccos(np.clip(
np.stack([(b ** 2 + c ** 2 - a ** 2) / (2 * b * c),
(a ** 2 + c ** 2 - b ** 2) / (2 * a * c),
(a ** 2 + b ** 2 - c ** 2) / (2 * a * b)]),
-1, 1)))
smallest = np.nanmin(angles, axis=0)
longest = np.maximum.reduce([a, b, c])
slivers = smallest < min_angle_deg
print(f" {len(triangulation.simplices):,} triangles, "
f"minimum angle median {np.median(smallest):.1f}°")
print(f" {int(slivers.sum()):,} slivers below {min_angle_deg}° "
f"({slivers.mean():.1%})")
print(f" longest edge median {np.median(longest):.0f} m, "
f"max {longest.max():.0f} m")
if slivers.mean() > 0.1:
print(" ! many slivers — the sample geometry is nearly collinear "
"somewhere, or clustered")
return triangulation, smallest, longest
_ = triangle_quality(points) 983 triangles, minimum angle median 29.4°
182 slivers below 15.0° (18.5%)
longest edge median 636 m, max 9330 m
! many slivers — the sample geometry is nearly collinear somewhere, or clustered
Slivers appear where samples are nearly collinear — along a road transect, or around the hull edge. Even the random sample here has 18.5% of its triangles below 15°, mostly along the hull, so the function's 10% warning fires on it; compare against a random sample of the same size before reading a high share as a design problem. Interpolation inside them is dominated by two of the three vertices and is numerically sensitive.
The longest-edge statistic is also useful: a triangle spanning five kilometres is a linear guess across a gap, and its interior is extrapolation in all but name.
Example 3 — masking long triangles as unsupported
import numpy as np
from scipy.spatial import Delaunay
def tin_with_edge_limit(points, values, targets, max_edge):
"""Refuse to interpolate inside triangles that span more than max_edge."""
from scipy.interpolate import LinearNDInterpolator
triangulation = Delaunay(points)
surface = LinearNDInterpolator(triangulation, values)(targets)
vertices = points[triangulation.simplices]
edges = np.stack([
np.linalg.norm(vertices[:, 1] - vertices[:, 0], axis=1),
np.linalg.norm(vertices[:, 2] - vertices[:, 1], axis=1),
np.linalg.norm(vertices[:, 0] - vertices[:, 2], axis=1)])
longest = edges.max(axis=0)
simplex = triangulation.find_simplex(targets)
too_long = (simplex >= 0) & (longest[simplex] > max_edge)
surface = np.where(too_long, np.nan, surface)
print(f" {int(too_long.sum()):,} targets fall in triangles longer than "
f"{max_edge} m ({too_long.mean():.1%}) — masked")
print(f" final coverage {np.isfinite(surface).mean():.1%}")
return surface
_ = tin_with_edge_limit(points, values, grid_points, max_edge=1500) 5,179 targets fall in triangles longer than 1500 m (4.7%) — masked
final coverage 93.0%
Example 4 — the table at the top of the page
pool = np.setdiff1d(np.arange(truth.size), test_idx)
err = lambda pred, mask: float(np.sqrt(np.mean((pred[mask] - truth[mask]) ** 2)))
print("samples TIN coverage TIN RMSE IDW RMSE IDW on TIN's cells")
for n in (100, 250, 500, 1000, 2000):
pick = np.random.default_rng(n).choice(pool, n, replace=False)
tin = griddata(grid_points[pick], truth[pick], grid_points, method="linear")
inside = np.isfinite(tin)
idw_pred = idw(grid_points[pick], truth[pick], grid_points)[0]
print(f"{n:7,} {inside.mean():13.1%} {err(tin, inside):9.2f} m {err(idw_pred, inside | ~inside):8.2f} m "
f"{err(idw_pred, inside):9.2f} m")samples TIN coverage TIN RMSE IDW RMSE IDW on TIN's cells
100 94.5% 106.28 m 101.77 m 102.01 m
250 94.3% 69.70 m 79.10 m 73.92 m
500 97.2% 47.56 m 54.92 m 54.91 m
1,000 99.0% 35.06 m 39.60 m 39.60 m
2,000 99.4% 24.62 m 28.71 m 28.70 m
This is the TIN equivalent of a maximum search distance. A triangle whose longest edge is five kilometres does not describe the terrain between its vertices any better than a straight line does, and masking it is more honest than shading it.
Bounded values, facets and the hull
Why a TIN cannot overshoot
Barycentric coordinates within a triangle are non-negative and sum to one, so the interpolated value is a convex combination of the three vertex values.
A convex combination is bounded by its inputs. So the surface never exceeds the maximum sample or falls below the minimum, anywhere.
That is a genuine advantage over kriging, which can leave the data range through negative weights — measured, every one of 116,137 grid targets had at least one negative weight — and over splines, which overshoot at sharp changes.
Why the surface has facets
Linear interpolation is continuous across triangle edges: two triangles sharing an edge agree along it, because both interpolate the same two vertices there.
The gradient is not continuous. Each triangle has a constant slope, and adjacent triangles have different ones, so the surface is piecewise planar.
That is invisible in a contour map and obvious in a hillshade, which is a function of the gradient. A TIN hillshade shows triangular facets, and that is not a bug — it is the honest appearance of a piecewise-linear surface.
For a smooth-looking result, use cubic interpolation over the same triangulation, and accept that it can overshoot.
Why the hull limitation is a feature
Every other common interpolator produces a value everywhere in its bounding box. Some of that is extrapolation, rendered identically to the interpolation.
Measured on a clustered sample in cross-validating an interpolation, IDW's error was 21.2 m within 100 m of a sample and 258.7 m beyond a kilometre, with nothing in the output distinguishing them.
A TIN returns NaN outside the hull, which is visible in every viewer and propagates through every calculation. It is a weaker guarantee than a proper support mask — a point inside the hull can still be far from any sample — but it is free and it is honest.
Why breaklines matter
Delaunay triangulation depends only on point positions, so it happily connects a point at the top of a cliff to one at the bottom, and interpolates a ramp between them.
For terrain that removes exactly the features that matter: levees, cuttings, retaining walls, stream banks.
A constrained triangulation forces edges along specified lines, so no triangle spans them. It requires a library beyond scipy, and for engineering-grade terrain it is not optional.
Collinear points, cubic overshoot and constraints
- Everything outside the convex hull is NaN. That is the method working.
- A TIN cannot overshoot — the surface is bounded by the sample values.
- The gradient is discontinuous, so hillshades show facets.
- Duplicate points break the triangulation. Aggregate first.
- Collinear points produce slivers where interpolation is unstable.
- Mask long triangles; a five-kilometre edge is extrapolation in all but name.
griddata(method="cubic")overshoots — check the output range.- Use a constrained triangulation for breaklines;
scipydoes not provide one.
Internal links
- IDW, kriging, splines or TIN? Choosing an interpolator — the measured comparison
- Spatial interpolation explained — the error-versus-distance picture
- My interpolated surface extends far beyond the data — why the hull is not enough
- How to interpolate points to a grid with IDW in Python — the method that fills everything
- How to fit a variogram and krige a surface in Python — the method that can overshoot
- How to create a DTM from LiDAR ground points in Python — TINs in terrain modelling
- How to choose a cell size for an interpolated surface — sizing the output
- Digital elevation models explained: DEM, DSM and DTM — where TINs fit
FAQ
How do I build a TIN in Python?
scipy.spatial.Delaunay for the triangulation and scipy.interpolate.LinearNDInterpolator — or griddata(method="linear") — for the interpolation.
Why is my TIN surface full of NaN?
Everything outside the convex hull of the samples is undefined. At 100 samples that was 5.5% of the study area; at 2,000 it was 0.6%.
Is a TIN better than IDW?
On the measured data it had a lower RMSE at every sample size above 100 — but it declines to answer for the hardest cells, so the comparison is not exactly like for like.
Why does my TIN hillshade look faceted?
Because the surface is piecewise planar: each triangle has a constant gradient. The surface is continuous; its gradient is not.
Can a TIN overshoot?
Linear interpolation cannot — the value is a convex combination of three vertices. Cubic interpolation over the same triangulation can.
How do I handle cliffs and levees?
With a constrained triangulation that forces edges along breaklines. scipy does not provide one; triangle and CGAL bindings do.
Should I fill the area outside the hull?
Only knowingly, with a different method, and recording which cells were filled. Refusing to extrapolate is the method's main advantage.