Comparing TIN, IDW, and Kriging for DTM Interpolation

Once ground returns are classified, a bare-earth surface still has to be interpolated onto a regular grid — and the interpolation method leaves a distinct fingerprint on the resulting terrain. This guide compares the three interpolators that matter for forestry DTMs: TIN (Delaunay triangulation with linear interpolation), IDW (inverse distance weighting), and ordinary kriging. It is one focused decision inside Digital Terrain Model Generation, part of the broader Canopy Height Modeling & Terrain Extraction workflow, and it picks up exactly where ground classification in LiDAR Point Cloud Preprocessing leaves off. For the fuller pipeline treatment of resolution and validation, see the sibling guide Generating high-res DTM from ALS data.

The three methods trade off differently between fidelity, smoothness, and compute cost. TIN honours every ground point exactly and is fast, but its linear facets produce visible triangle edges. IDW is simple and robust but flattens peaks and pits and can leave “bullseye” artifacts around isolated points. Ordinary kriging is the only method that models spatial autocorrelation explicitly and yields a per-cell error estimate, at a substantial compute cost that makes it impractical on a full survey without local neighbourhoods.

When to use each method

Match the interpolator to ground-return density, terrain character, and whether you need an uncertainty surface.

The practical rule: use TIN when ground returns are dense (4+ pts/m²) and you want the surface to pass exactly through every measurement; the linear facets are invisible at that density. Use IDW when density is moderate and you need a quick, robust fill on gentle terrain, accepting mild flattening. Use kriging when returns are sparse or irregular and the extra cost buys you a statistically optimal surface plus a mapped uncertainty layer — but confine it to local search neighbourhoods so the linear system stays solvable.

The math each method solves

IDW predicts an unsampled location as a weighted mean of nearby ground elevations, with weights that decay as a power $p$ of distance $d_i$:

$$\hat{z}({\mathbf{x}}_0)=\frac{{\sum }_{i=1}^n{w}_i{z}_i}{{\sum }_{i=1}^n{w}_i},w_i=\frac{1}{d_i^p}$$

A higher power $p$ (commonly 2) concentrates influence on the nearest returns and sharpens the surface; a lower power spreads influence and smooths it. IDW has no notion of terrain trend, which is why isolated points radiate concentric bullseyes.

Ordinary kriging instead fits a semivariogram $\gamma (h)$ that describes how elevation variance grows with separation distance $h$, then solves for weights that minimize the prediction variance under an unbiasedness constraint. The empirical semivariogram is:

$$\gamma (h)=\frac{1}{2N(h)}\sum _{i=1}^{N(h)}{\left(z({\mathbf{x}}_i)-z({\mathbf{x}}_i+h)\right)}^2$$

where $N(h)$ is the number of point pairs separated by lag $h$. A model (spherical, exponential, Gaussian) is fitted to this cloud and supplies the covariances the kriging system needs. TIN needs no such model — it connects the ground points into non-overlapping Delaunay triangles and linearly interpolates within each facet, so the surface is exact at every vertex.

Minimal reproducible example

The fastest operational route is to let PDAL build the surface directly. writers.gdal with output_type="idw" grids classified ground returns with inverse distance weighting in a single stage; isolating Class 2 first keeps vegetation out of the surface.

{
  "pipeline": [
    "classified_tile.laz",
    {
      "type": "filters.range",
      "limits": "Classification[2:2]"
    },
    {
      "type": "writers.gdal",
      "filename": "dtm_idw_1m.tif",
      "resolution": 1.0,
      "output_type": "idw",
      "data_type": "float32",
      "gdalopts": "COMPRESS=DEFLATE,PREDICTOR=2,TILED=YES"
    }
  ]
}

For a TIN surface, PDAL builds a Delaunay triangulation with filters.delaunay and then rasterizes with filters.faceraster, which linearly interpolates each grid cell inside its containing triangle:

{
  "pipeline": [
    "classified_tile.laz",
    {
      "type": "filters.range",
      "limits": "Classification[2:2]"
    },
    {
      "type": "filters.delaunay"
    },
    {
      "type": "filters.faceraster",
      "resolution": 1.0
    },
    {
      "type": "writers.raster",
      "filename": "dtm_tin_1m.tif",
      "data_type": "float32"
    }
  ]
}

When you need full control — or an ordinary-kriging surface, which PDAL does not provide — extract the ground returns to NumPy and interpolate in Python. scipy.interpolate.griddata gives you both linear TIN and a nearest baseline; pykrige fits a variogram and krigs:

import numpy as np
import pdal
import json
from scipy.interpolate import griddata

def load_ground_xyz(classified_laz: str) -> np.ndarray:
    """Return an (n, 3) array of classified ground (Class 2) x, y, z."""
    pipeline = {
        "pipeline": [
            classified_laz,
            {"type": "filters.range", "limits": "Classification[2:2]"},
        ]
    }
    p = pdal.Pipeline(json.dumps(pipeline))
    p.execute()
    arr = p.arrays[0]
    return np.column_stack([arr["X"], arr["Y"], arr["Z"]])

def tin_grid(xyz: np.ndarray, res: float = 1.0) -> np.ndarray:
    """Linear (TIN) interpolation of ground z onto a regular grid."""
    xmin, ymin = xyz[:, 0].min(), xyz[:, 1].min()
    xmax, ymax = xyz[:, 0].max(), xyz[:, 1].max()
    gx, gy = np.meshgrid(
        np.arange(xmin, xmax, res),
        np.arange(ymin, ymax, res),
    )
    dtm = griddata(xyz[:, :2], xyz[:, 2], (gx, gy), method="linear")
    return dtm

def idw_grid(xyz: np.ndarray, res: float = 1.0, power: float = 2.0,
             k: int = 12) -> np.ndarray:
    """Inverse distance weighting with a k-nearest-neighbour search."""
    from scipy.spatial import cKDTree
    xmin, ymin = xyz[:, 0].min(), xyz[:, 1].min()
    xmax, ymax = xyz[:, 0].max(), xyz[:, 1].max()
    gx, gy = np.meshgrid(
        np.arange(xmin, xmax, res),
        np.arange(ymin, ymax, res),
    )
    targets = np.column_stack([gx.ravel(), gy.ravel()])
    tree = cKDTree(xyz[:, :2])
    dist, idx = tree.query(targets, k=k)
    dist = np.maximum(dist, 1e-6)          # guard against divide-by-zero
    w = 1.0 / dist**power
    z = np.sum(w * xyz[idx, 2], axis=1) / np.sum(w, axis=1)
    return z.reshape(gx.shape)

# Ordinary kriging (confine to local neighbourhoods for speed):
#   from pykrige.ok import OrdinaryKriging
#   ok = OrdinaryKriging(xyz[:, 0], xyz[:, 1], xyz[:, 2],
#                        variogram_model="spherical", nlags=12)
#   dtm, var = ok.execute("grid", grid_x, grid_y)   # var is the error surface

ground = load_ground_xyz("classified_tile.laz")
dtm_tin = tin_grid(ground)
dtm_idw = idw_grid(ground)
print("TIN NaNs:", np.isnan(dtm_tin).sum(), " IDW range:",
      round(float(np.nanmin(dtm_idw)), 2), round(float(np.nanmax(dtm_idw)), 2))

Parameter reference

The interpolation knobs that most affect DTM quality, with starting values and rationale.

Expected output and verification

All three routes produce a single-band float32 GeoTIFF of bare-earth elevation. Regardless of method, validate against independent RTK or ground-control points and compare RMSE — the surface with the lowest RMSE and a mean error nearest zero wins for that terrain. A quick internal check confirms the interpolated range stays within the observed ground-return range (extrapolation beyond it signals a bad power or variogram fit):

import numpy as np

def sanity_check(dtm: np.ndarray, ground_z: np.ndarray) -> None:
    """Interpolated elevations must not stray outside the observed ground range."""
    zmin, zmax = float(ground_z.min()), float(ground_z.max())
    finite = dtm[np.isfinite(dtm)]
    assert finite.min() >= zmin - 0.5, "DTM dips below observed ground — check IDW power / variogram"
    assert finite.max() <= zmax + 0.5, "DTM rises above observed ground — extrapolation artifact"
    print(f"ok: DTM in [{finite.min():.2f}, {finite.max():.2f}] within ground [{zmin:.2f}, {zmax:.2f}]")

sanity_check(dtm_idw, ground[:, 2])

For kriging specifically, inspect the returned variance surface: cells far from any ground return carry high prediction variance and should be treated as low-confidence or masked out entirely, exactly the wide-gap regions where all three methods are least trustworthy.

Common pitfalls

• IDW bullseyes around isolated points. With sparse ground and power=2, each isolated return radiates concentric rings. Lower the power toward 1.5, widen the neighbour count, or switch to TIN, which interpolates linearly across the gap instead.
• TIN flat triangles over voids. Delaunay spans wide sparse-return gaps with a single large flat facet that reads as false plateau terrain. Fill only narrow gaps and mask large ones rather than trusting the facet.
• Kriging on the full tile. A global kriging system is $O(n^3)$ and will not finish on millions of returns. Always confine kriging to a local search neighbourhood (a moving window of the nearest k points).
• Extrapolating past the data hull. griddata linear returns NaN outside the convex hull, but IDW and kriging will happily invent elevations there. Clip the surface to the ground-return footprint and leave the rest NoData.

Frequently Asked Questions

Related

• Digital Terrain Model Generation — the parent workflow this interpolation choice belongs to
• Generating high-res DTM from ALS data — the sibling guide on resolution and validation in depth
• LiDAR Point Cloud Preprocessing — produces the classified ground returns this step interpolates
• Canopy Height Modeling & Terrain Extraction — the full point-cloud-to-canopy pipeline

Up: Digital Terrain Model Generation · Canopy Height Modeling & Terrain Extraction