Python LiDAR

Measuring Swath-to-Swath Relative Accuracy

TL;DR: For each pair of overlapping flightlines, rasterize single returns of each line separately at 1 m (filters.range on PointSourceId and NumberOfReturns[1:1], writers.gdal with output_type: "mean"), keep only cells that are flat and hard in both, difference the two rasters, and report the root-mean-square difference (RMSDz) and the largest absolute difference. Values of a few centimetres are typical of well-calibrated data.

# Context and Motivation

This guide is part of Vertical Accuracy Assessment for LiDAR. Checkpoints measure absolute accuracy at a few dozen places. Relative accuracy — agreement between overlapping swaths — measures internal consistency across the whole project, needs no field survey, and is sensitive to exactly the problems checkpoints miss: a mis-calibrated boresight that tilts one line, a trajectory drift late in a flight, a line flown on a different day with a different GNSS solution. Where two swaths overlap, both measured the same ground, so any difference between them is error.

Delivery specifications such as the USGS Lidar Base Specification set limits on swath-to-swath differences per quality level. The measurement itself is straightforward with PDAL and NumPy once you restrict it to surfaces where differences mean something.

The overlap is a free accuracy test Plan view of two parallel flightline swaths overlapping in a central strip. In the strip, cells on flat hard surfaces are shaded as usable; cells on slopes and vegetation are excluded. Beside it, a cross-section shows swath A slightly above swath B over the same road surface, with the vertical difference labelled. swath A swath B usable cells ΔZ between swaths same road, two flightlines

# Prerequisites and Assumptions

  • Point data with PointSourceId populated per flightline — the standard LAS convention. If IDs were lost when tiles were merged, relative accuracy cannot be measured from the tiles.
  • Ground classified, for the slope mask.
  • PDAL 2.x, Python with NumPy and rasterio.
  • Knowledge of which lines overlap: from a flight plan, or computed from line footprints.

# Step-by-Step Implementation

# Step 1 — List the swaths in a tile

pdal info --stats --enumerate PointSourceId lists the flightline IDs present.

# Step 2 — Rasterize each swath separately

For each ID, keep single returns (NumberOfReturns == 1), which are almost always hard surfaces, and write a mean-Z raster on a common grid with fixed origin_x, origin_y, width and height so cells align exactly.

# Step 3 — Build a mask of usable cells

Keep cells where both swaths have data, both have at least a few returns, and terrain slope from a DTM is under about 10 degrees. Slopes turn small horizontal misalignment into large vertical differences that are not what you are measuring.

# Step 4 — Difference and summarize

Compute ΔZ = Z_A − Z_B over the mask, then RMSDz, mean difference and the largest absolute difference (or a high percentile, which is more robust).

# Step 5 — Map the differences

Write the difference raster. Spatial patterns — a gradient across the overlap, a step at one end — point to specific calibration problems.

# Complete Working Example

python
"""Swath-to-swath relative accuracy for one tile and one pair of flightlines."""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import pdal
import rasterio

RES = 1.0


def grid_spec(src: Path) -> dict:
    b = pdal.Pipeline(json.dumps({"pipeline": [str(src)]})).quickinfo["readers.las"]["bounds"]
    ox, oy = np.floor(b["minx"]), np.floor(b["miny"])
    return {"origin_x": ox, "origin_y": oy,
            "width": int(np.ceil((b["maxx"] - ox) / RES)),
            "height": int(np.ceil((b["maxy"] - oy) / RES))}


def swath_raster(src: Path, psid: int, out: Path, grid: dict) -> None:
    pdal.Pipeline(json.dumps({"pipeline": [
        str(src),
        {"type": "filters.expression",
         "expression": f"PointSourceId == {psid} && NumberOfReturns == 1 && Classification != 7 && Classification != 18"},
        {"type": "writers.gdal", "filename": str(out), "resolution": RES, "radius": RES * 0.71,
         "output_type": "mean,count", "data_type": "float32", "nodata": -9999, **grid},
    ]})).execute()


def slope_deg(dtm: Path) -> np.ndarray:
    with rasterio.open(dtm) as ds:
        z = ds.read(1, masked=True).filled(np.nan)
        gy, gx = np.gradient(z, ds.res[1], ds.res[0])
    return np.degrees(np.arctan(np.hypot(gx, gy)))


def relative_accuracy(src: Path, a: int, b: int, dtm: Path, max_slope: float = 10.0) -> dict:
    grid = grid_spec(src)
    ra, rb = Path(f"out/swath_{a}.tif"), Path(f"out/swath_{b}.tif")
    swath_raster(src, a, ra, grid)
    swath_raster(src, b, rb, grid)
    with rasterio.open(ra) as da, rasterio.open(rb) as db:
        za, ca = da.read(1).astype(float), da.read(2)
        zb, cb = db.read(1).astype(float), db.read(2)
        profile = da.profile
    slope = slope_deg(dtm)
    mask = (za > -9999) & (zb > -9999) & (ca >= 2) & (cb >= 2) & (slope < max_slope)
    dz = (za - zb)[mask]
    diff = np.where(mask, za - zb, -9999).astype("float32")
    profile.update(count=1)
    with rasterio.open(f"out/dz_{a}_{b}.tif", "w", **profile) as out:
        out.write(diff, 1)
    return {"pair": (a, b), "cells": int(mask.sum()),
            "mean_m": round(float(dz.mean()), 3),
            "rmsdz_m": round(float(np.sqrt(np.mean(dz ** 2))), 3),
            "p99_abs_m": round(float(np.percentile(np.abs(dz), 99)), 3),
            "max_abs_m": round(float(np.abs(dz).max()), 3)}


if __name__ == "__main__":
    Path("out").mkdir(exist_ok=True)
    print(relative_accuracy(Path("tiles/t_0431.laz"), 1102, 1103, Path("dtm/t_0431_dtm.tif")))

output_type: "mean,count" writes two bands; the count band enforces that each cell was actually measured by both swaths rather than filled by the writer’s radius.

Reading the difference map Three small difference maps of an overlap strip. A uniform offset across the strip suggests a vertical bias between lines, such as a GNSS solution difference. A gradient across the strip, from negative at one edge to positive at the other, suggests a roll or boresight error. A gradient along the strip suggests trajectory drift over time. uniform offset vertical bias between lines gradient across roll or boresight gradient along trajectory drift

# Key Parameter Table

Parameter Type Default Guidance
RES float, m 1.0 Coarser averages more returns per cell; finer shows detail
returns used expression NumberOfReturns == 1 Single returns are mostly hard surfaces
minimum count per cell int 2 Rejects cells filled by one stray return
max_slope degrees 10 Slopes convert horizontal error into vertical difference
summary statistic RMSDz, p99, max all three Max is fragile; p99 is steadier for reports

# Verification

  • Self-difference is zero. Running the function with the same ID twice must give zero differences — a check on grid alignment.
  • Overlap area plausible. The count of usable cells should correspond to the expected overlap width times tile length, reduced by vegetation and slopes.
  • Pairs sum sensibly. For three mutually overlapping lines A, B, C, mean(A−B) + mean(B−C) should approximate mean(A−C). Large inconsistencies mean the masks differ strongly between pairs.

# Gotchas and Edge Cases

Merged or reset PointSourceId. Some processing chains overwrite PointSourceId with a tile number. Check with an enumerate first; if every point in a tile has the same ID, relative accuracy must be measured from the original swath files.

Vegetation and edges. Even single returns include some canopy tops and roof edges. The slope mask removes many, but also mask by classification (ground and buildings only) if differences look noisy.

Moving surfaces. Water, vehicles, crops that were harvested between flights and construction sites differ between swaths for real reasons. Exclude them with the classification or known change areas.

Why slopes are masked Two copies of a sloping surface offset horizontally by 0.2 metres. On a 30 degree slope, that horizontal shift produces a vertical difference of about 0.12 metres at the same map location, even though neither swath has any vertical error. On flat ground the same shift produces no vertical difference. 0.12 m apparent ΔZ 30° slope, 0.2 m horizontal shift flat ground: no apparent ΔZ same 0.2 m shift

Grid alignment. Letting writers.gdal choose its own origin for each swath misaligns cells by a fraction of a metre. Always pass the same origin_x, origin_y, width and height to every swath raster.

# Frequently Asked Questions

What is swath-to-swath relative accuracy?

The agreement between overlapping flightlines of the same LiDAR collection, measured as the vertical difference between their surfaces in the overlap. It reflects calibration and trajectory quality and needs no ground survey.

Why use only single returns?

Single returns come mostly from hard, opaque surfaces such as roads, roofs and bare ground, where both swaths should record the same elevation. Multiple returns from vegetation differ between swaths because of viewing geometry, not error.

What values are acceptable?

Specifications set thresholds by quality level; well-calibrated modern airborne systems typically show differences of a few centimetres on flat hard surfaces. Consult the specification your project follows for exact limits on RMSDz and maximum difference.

Can I measure relative accuracy on tiles instead of original swaths?

Yes, if PointSourceId still identifies the flightline of each point. Some processing chains overwrite it; check the values before relying on tiles.

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... Interpolating LiDAR Elevations at Checkpoints --- title: "Interpolating LiDAR Elevations at Checkpoints" description: "Get the LiDAR elevation at a surveyed checkpoint the way accuracy... 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...