Python LiDAR

Interpolating LiDAR Elevations at Checkpoints

TL;DR: Collect class 2 returns within a few metres of each checkpoint, build a Delaunay triangulation with SciPy, find the triangle containing the checkpoint and interpolate Z with barycentric weights. Reject the checkpoint (and say so) if it falls outside the triangulation or if the enclosing triangle is larger than a few times the nominal point spacing — the surface there is an extrapolation, not a measurement.

# Context and Motivation

This guide is part of Vertical Accuracy Assessment for LiDAR. A checkpoint almost never coincides with a LiDAR return, so an accuracy assessment has to estimate what the LiDAR says the ground elevation is at the checkpoint’s exact horizontal position. The method matters more than people expect. Nearest-neighbour lookup mixes horizontal offset into vertical error. A DEM lookup adds the DEM’s own interpolation and cell-size smoothing. Inverse-distance weighting pulls toward whichever returns happen to be clustered on one side. A triangulated irregular network of ground returns, interpolated linearly inside the enclosing triangle, is what accuracy standards describe and what gives the most defensible number.

The implementation is short; the value is in the rejection rules, which keep a handful of badly supported checkpoints from contaminating the statistic.

Barycentric interpolation inside one triangle Plan view of ground returns triangulated into a mesh. A checkpoint sits inside one triangle whose vertices have elevations 212.41, 212.58 and 212.36 metres. Lines from the checkpoint to each vertex show the barycentric weights 0.52, 0.31 and 0.17, which combine the three elevations into the interpolated value 212.45 metres. vertex A 212.41 m · weight 0.52 vertex B 212.58 m · weight 0.31 vertex C 212.36 m · weight 0.17 Z at checkpoint = Σ weight × Z = 212.45 m weights are the checkpoint's barycentric coordinates; they sum to one and are all positive inside the triangle

# Prerequisites and Assumptions

  • Ground-classified LiDAR (class 2) in the same horizontal and vertical CRS as the checkpoints.
  • Checkpoint coordinates with enough precision that their horizontal position is known to a few centimetres.
  • Python with NumPy, SciPy and PDAL bindings.
  • Nominal point spacing of the ground returns, used to judge triangle size.

# Step-by-Step Implementation

# Step 1 — Gather local ground

Crop class 2 points to a square of about 3 m around the checkpoint (more under canopy, where ground is sparse).

# Step 2 — Triangulate

scipy.spatial.Delaunay on the X, Y of those points. Local triangulations are fast and avoid the numerical trouble of triangulating millions of points.

# Step 3 — Locate and interpolate

find_simplex returns the triangle containing the checkpoint; barycentric weights from Delaunay.transform interpolate Z linearly.

# Step 4 — Apply rejection rules

Reject if the checkpoint is outside the hull, if the triangle’s longest edge exceeds about three times the nominal spacing, or if fewer than a minimum number of ground points were found.

# Step 5 — Record diagnostics

Store the triangle edge length, number of points and distance to the nearest return with each result, so reviewers can see how well supported each value is.

# Complete Working Example

python
"""TIN interpolation of LiDAR ground elevation at a checkpoint, with rejection rules."""
from __future__ import annotations

import json
from dataclasses import dataclass

import numpy as np
import pdal
from scipy.spatial import Delaunay


@dataclass
class Interp:
    z: float
    status: str
    n_points: int
    max_edge_m: float
    nearest_m: float


def local_ground(tile: str, x: float, y: float, r: float) -> np.ndarray:
    p = pdal.Pipeline(json.dumps({"pipeline": [
        tile,
        {"type": "filters.crop", "bounds": f"([{x - r}, {x + r}], [{y - r}, {y + r}])"},
        {"type": "filters.range", "limits": "Classification[2:2]"},
    ]}))
    p.execute()
    return p.arrays[0]


def interpolate(pts: np.ndarray, x: float, y: float, spacing: float,
                min_points: int = 6, edge_factor: float = 3.0) -> Interp:
    n = len(pts)
    if n < min_points:
        return Interp(np.nan, f"too few ground points ({n})", n, np.nan, np.nan)
    xy = np.column_stack([pts["X"], pts["Y"]])
    nearest = float(np.min(np.hypot(xy[:, 0] - x, xy[:, 1] - y)))
    tri = Delaunay(xy)
    s = int(tri.find_simplex(np.array([[x, y]]))[0])
    if s < 0:
        return Interp(np.nan, "outside triangulation", n, np.nan, nearest)
    v = tri.simplices[s]
    corners = xy[v]
    edges = np.linalg.norm(corners - np.roll(corners, 1, axis=0), axis=1)
    max_edge = float(edges.max())
    if max_edge > edge_factor * spacing:
        return Interp(np.nan, f"triangle too large ({max_edge:.2f} m)", n, max_edge, nearest)
    T = tri.transform[s]
    b = T[:2] @ (np.array([x, y]) - T[2])
    w = np.append(b, 1.0 - b.sum())
    return Interp(float(w @ pts["Z"][v]), "ok", n, max_edge, nearest)


if __name__ == "__main__":
    cp = {"id": "CP-017", "x": 431522.314, "y": 4471380.902, "z": 212.418}
    g = local_ground("tiles/t_0431.laz", cp["x"], cp["y"], r=3.0)
    res = interpolate(g, cp["x"], cp["y"], spacing=0.35)
    print(cp["id"], res, "dz =", None if np.isnan(res.z) else round(res.z - cp["z"], 3))
Why TIN beats nearest neighbour on slopes A profile of sloping ground with returns at irregular spacing. A checkpoint lies between two returns. The TIN estimate lies on the line between them, matching the ground. The nearest-neighbour estimate takes the elevation of the closest return, which is uphill, producing a vertical error of 0.14 metres. An inverse-distance estimate using a cluster of returns on one side is pulled toward them. TIN: on the surface nearest return: +0.14 m IDW: pulled toward the cluster the checkpoint's true ground elevation is where the dashed line meets the slope

# Key Parameter Table

Parameter Type Default Guidance
search half-width r float, m 3.0 2–3 × spacing minimum; 5–10 m under dense canopy
min_points int 6 Enough for a local TIN with a few triangles
spacing float, m nominal ground spacing Sets the triangle-size rejection threshold
edge_factor float 3.0 Longest edge allowed as a multiple of spacing
method TIN — Nearest neighbour and raster lookups add their own error

# Verification

  • Synthetic plane. Generate points on an inclined plane, add a checkpoint on the plane, and confirm interpolation returns the exact plane elevation. Nearest neighbour will not.
  • Diagnostics distribution. Plot max_edge_m for all accepted checkpoints; a long tail indicates sparse ground that the rejection threshold might be allowing through.
  • Sensitivity to radius. Rerun with twice the search radius. Accepted checkpoints should not change elevation; if they do, the local TIN was on the hull edge.
python
rng = np.random.default_rng(3)
xs, ys = rng.uniform(0, 6, 60), rng.uniform(0, 6, 60)
plane = np.zeros(60, dtype=[("X", "f8"), ("Y", "f8"), ("Z", "f8")])
plane["X"], plane["Y"], plane["Z"] = xs, ys, 100 + 0.3 * xs - 0.1 * ys
r = interpolate(plane, 3.0, 3.0, spacing=0.8)
assert abs(r.z - (100 + 0.3 * 3.0 - 0.1 * 3.0)) < 1e-9

# Gotchas and Edge Cases

Breaklines near the checkpoint. A checkpoint near a kerb, ditch or wall can land in a triangle spanning the break, averaging two surfaces. Checkpoints should be sited away from breaklines; if one is not, flag it rather than interpolating across the edge.

Misclassified ground. A few low vegetation or building-edge returns in class 2 near a checkpoint bias the TIN. Inspecting the local points for a sample of checkpoints — especially outliers — often reveals classification issues worth reporting.

Checkpoints on the hull. Where a checkpoint sits near the edge of the local point set, the enclosing triangle is long and thin. The edge-length test catches most of these; increasing the search radius resolves the rest.

When the triangle is too big Left: a checkpoint in a small, compact triangle with edges under one metre, accepted. Right: a checkpoint in a gap under dense canopy, enclosed by a long thin triangle with a 4.1 metre edge, more than three times the nominal spacing, rejected as an extrapolation. edges < 1 m: accepted 4.1 m edge > 3 × spacing: rejected

Rasters as a shortcut. Sampling a 1 m DTM at checkpoints is fine for assessing the DTM deliverable, but it measures the DTM’s interpolation and smoothing too. Report point-cloud accuracy from the TIN and DEM accuracy from the raster, labelled as such.

# Frequently Asked Questions

How should LiDAR elevation be interpolated at a checkpoint?

From a triangulated irregular network of ground-classified returns around the checkpoint, interpolating linearly within the triangle that contains it. This estimates the surface at the checkpoint’s exact position without horizontal offsets or raster smoothing.

Why not use the nearest LiDAR point?

The nearest return can be a metre or more away horizontally. On any slope, that horizontal distance creates a vertical difference that is not LiDAR error, inflating RMSEz.

When should a checkpoint be rejected?

When too few ground returns surround it, when it lies outside the local triangulation, or when the enclosing triangle is much larger than the nominal point spacing. Report rejected checkpoints and the reason rather than silently dropping them.

Can I use a DTM raster instead?

For assessing the DTM product, yes. For assessing the point cloud itself, a TIN of ground points is preferred because the raster adds its own interpolation error.

Related in Vertical Accuracy Assessment for LiDAR

↑ Vertical Accuracy Assessment for LiDAR overview --- title: "Vertical Accuracy Assessment for LiDAR" description: "Measure and report LiDAR vertical accuracy in Python: surveyed... Computing RMSEz Against Survey Checkpoints --- title: "Computing RMSEz Against Survey Checkpoints" description: "Calculate RMSEz, mean error and standard deviation of LiDAR elevations... Reporting NVA and VVA Accuracy --- title: "Reporting NVA and VVA Accuracy" description: "Compute and report non-vegetated and vegetated vertical accuracy for LiDAR... Measuring Swath-to-Swath Relative Accuracy --- title: "Measuring Swath-to-Swath Relative Accuracy" description: "Check internal consistency between overlapping LiDAR flightlines...