Tiling and Buffering Large LiDAR Surveys with PDAL

Every neighbourhood operation in a LiDAR pipeline — ground classification, outlier removal, interpolation — needs points beyond the pixel it is computing. Process a tile in isolation and those operations run out of data at the edge, producing a visible ridge of wrong values along every tile boundary that survives all the way into the finished terrain and canopy products. The fix is a halo: process a tile with a margin of borrowed neighbouring points, then crop the margin away. This guide covers exactly that, as one operational step inside LiDAR Point Cloud Preprocessing, part of Canopy Height Modeling & Terrain Extraction.

When to use this approach

Tiling is not optional above a certain survey size, but how you tile depends on what you are running.

Use filters.splitter for anything that becomes a raster: an aligned square grid mosaics without resampling. Reserve filters.chipper for point-domain work — segmentation, statistics — where equal work per chip matters more than a tidy grid.

Sizing the halo

The halo has to be at least as wide as the widest neighbourhood any stage will reach for. Take the maximum over the stages in the pipeline:

$$b=\max \left(w_{\text{smrf}},k_{\text{outlier}}\cdot s,2r_{\text{crown}},d_{\text{interp}}\right)$$

where $w_{\text{smrf}}$ is the ground filter’s largest morphological window, $s$ the mean point spacing, $r_{\text{crown}}$ the largest crown radius if you will segment trees, and $d_{\text{interp}}$ the void-fill search distance. For a typical forestry pipeline — SMRF at an 18 m window, crowns up to 8 m radius, a 6 m fill distance — a 30 m halo covers everything, and 50 m is a safe default that costs little.

Minimal reproducible example

Two pipelines do the work. The first splits a survey into a buffered grid; the second processes one buffered tile and crops on write.

Splitting with an explicit halo — filters.splitter takes the buffer directly, so neighbouring points are duplicated into each tile rather than dropped:

{
  "pipeline": [
    "survey.laz",
    {
      "type": "filters.splitter",
      "length": 1000,
      "buffer": 50,
      "origin_x": 500000,
      "origin_y": 5400000
    },
    {
      "type": "writers.las",
      "filename": "tiles/tile_#.laz",
      "forward": "all",
      "compression": "laszip"
    }
  ]
}

Pinning origin_x and origin_y to a round coordinate is what makes the outputs mosaic without resampling later. Leave the origin unset and PDAL anchors the grid to the data’s own bounding box, so two surveys of the same forest tile differently and every multi-temporal comparison inherits a half-pixel shift.

Processing one buffered tile and cropping back on write:

import json
import pdal


def process_tile(in_laz: str, out_tif: str, bounds: tuple, buffer_m: float = 50.0):
    """Classify + rasterize a buffered tile, writing only the nominal extent.

    bounds is the NOMINAL (unbuffered) extent as (minx, maxx, miny, maxy).
    """
    minx, maxx, miny, maxy = bounds
    pipeline = {
        "pipeline": [
            in_laz,
            {"type": "filters.outlier", "method": "statistical",
             "mean_k": 12, "multiplier": 2.2},
            {"type": "filters.smrf", "slope": 0.2, "window": 18.0,
             "scalar": 1.2, "threshold": 0.45},
            {"type": "filters.range", "limits": "Classification[2:2]"},
            # Crop AFTER every neighbourhood operation, never before.
            {"type": "filters.crop",
             "bounds": f"([{minx},{maxx}],[{miny},{maxy}])"},
            {"type": "writers.gdal", "filename": out_tif, "resolution": 1.0,
             "output_type": "min", "data_type": "float32",
             "bounds": f"([{minx},{maxx}],[{miny},{maxy}])",
             "gdalopts": "COMPRESS=DEFLATE,PREDICTOR=2,TILED=YES"},
        ]
    }
    pdal.Pipeline(json.dumps(pipeline)).execute()

The ordering is the whole point: filters.crop sits after the classifier and before the writer. Crop first and you have thrown away the context the classifier needed; do not crop at all and neighbouring rasters overlap, so the mosaic double-counts along every seam.

Deriving the nominal bounds for each tile from its own header, with the halo subtracted:

import json
import subprocess


def nominal_bounds(tile_laz: str, buffer_m: float = 50.0) -> tuple:
    """Read a buffered tile's extent and shrink it back to the nominal grid cell."""
    info = json.loads(subprocess.run(
        ["pdal", "info", "--metadata", tile_laz],
        capture_output=True, text=True, check=True).stdout)
    md = info["metadata"]
    return (md["minx"] + buffer_m, md["maxx"] - buffer_m,
            md["miny"] + buffer_m, md["maxy"] - buffer_m)

Parameter reference

A 50 m halo on a 1 km tile adds about 20% more points to read and process. On a 250 m tile the same halo adds 96% — the overhead is why very small tiles are a false economy.

Expected output and verification

Correct tiling produces outputs that abut exactly once. The two failures worth testing for are gaps (halo cropped too aggressively) and overlaps (not cropped at all), and both are visible in the extents alone:

import glob

import rasterio
from itertools import combinations


def check_mosaic(tifs: list[str], tol: float = 1e-6) -> dict:
    """Verify tile rasters share a grid, do not overlap, and leave no gaps."""
    boxes, res, crs = [], None, None
    for path in tifs:
        with rasterio.open(path) as src:
            boxes.append(src.bounds)
            if res is None:
                res, crs = abs(src.transform.a), src.crs
            assert abs(abs(src.transform.a) - res) < tol, f"{path}: resolution differs"
            assert src.crs == crs, f"{path}: CRS differs"
            # Pixel edges must land on the same global grid.
            assert abs(src.bounds.left / res - round(src.bounds.left / res)) < 1e-3, \
                f"{path}: origin off-grid — tiles will not mosaic cleanly"

    overlaps = 0
    for a, b in combinations(boxes, 2):
        if a.left < b.right - tol and b.left < a.right - tol and \
           a.bottom < b.top - tol and b.bottom < a.top - tol:
            overlaps += 1
    return {"n_tiles": len(boxes), "overlapping_pairs": overlaps, "res_m": res}


report = check_mosaic(sorted(glob.glob("dtm/*.tif")))
assert report["overlapping_pairs"] == 0, "tiles overlap — the halo was never cropped"

The stronger check is on values rather than extents: mosaic the tiles, then sample the finished raster along a tile boundary and along a line well inside a tile, and compare the distribution of local slope. A seam shows up as a slope spike confined to the boundary column.

Common pitfalls

• Cropping before classification. The most common ordering mistake, and it silently reproduces exactly the artifact the halo was added to prevent.
• Letting PDAL choose the grid origin. Two acquisitions over the same forest then tile differently, and every change-detection raster inherits a sub-pixel misalignment that looks like real growth.
• A halo narrower than the ground filter window. An 18 m SMRF window with a 10 m halo still starves at the edge; the artifact is smaller but not gone.
• Forgetting bounds on writers.gdal. Without it the raster extent follows the surviving points, so a tile with no returns in its top-left corner produces a raster smaller than its tile and the mosaic gets a hole.
• Assuming vendor tiles are already buffered. They almost never are. Delivery tiles are storage units, not processing units.

Frequently Asked Questions

Related

• LiDAR Point Cloud Preprocessing — the stage this tiling discipline serves
• Filtering Noise & Outliers from LiDAR with PDAL — one of the neighbourhood operations that needs the halo
• Digital Terrain Model Generation — where an un-buffered seam first becomes visible
• Mosaicking CHM Tiles Without Edge Seams — the raster-side counterpart to this page

Up: LiDAR Point Cloud Preprocessing · Canopy Height Modeling & Terrain Extraction