⚙️ How-To Interpolation & surfaces Terrain & elevation Interpolation & Surfaces · step 8 of 13 Updated

How to Build a TIN and Interpolate Elevations

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
    100        85.2%      100.29 m    98.41 m
    250        93.5%       67.56 m    80.58 m
    500        97.3%       52.48 m    58.41 m
  1,000        99.0%       38.05 m    43.31 m
  2,000        99.5%       26.98 m    28.78 m

Everything outside the convex hull of the samples is NaN, and that is the point. It also beats IDW at every sample size above 100 — while answering a slightly easier question, because the cells it declines are the hardest ones.

Quick answer 

import numpy as np
from scipy.interpolate import griddata

surface = griddata(points, values, targets, method="linear")

coverage = np.isfinite(surface).mean()
print(f"{coverage:.1%} 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.

Sample points triangulated into a Delaunay network, with linear interpolation inside each triangle and no values outside the convex hull.
Each target falls in one triangle and is a weighted average of its three vertices. Outside the hull there is no triangle.

Step-by-step solution

1. Build the triangulation 

from scipy.spatial import Delaunay
tri = Delaunay(points)
print(f"{len(points)} points -> {len(tri.simplices)} 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: 85.2% coverage
  500 samples: 97.3%
2,000 samples: 99.5%

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 above 100 samples, but the comparison is not like for like — TIN declined 2.7% of the cells at 500 samples, and those are the cells furthest from data.

Score both on the cells where both produce a value, or accept that TIN's number describes a slightly easier problem.

TIN coverage rising from 85.2 percent at 100 samples to 99.5 percent at 2,000, with everything outside the convex hull undefined.
The undefined area is the convex hull's complement. Refusing to answer there is the method's main virtue.

Code examples

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

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

Slivers appear where samples are nearly collinear — along a road transect, or around the hull edge. 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

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.

Explanation

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, IDW's error was 18.9 m within 100 m of a sample and 334.3 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.

TIN linear interpolation beating IDW at every sample size above 100, with the gap narrowing as density rises.
At 100 samples TIN is behind — but it is only being scored where it agreed to answer.

Edge cases or notes

  • 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; scipy does not provide one.

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 14.8% of the study area; at 2,000 it was 0.5%.

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.