Python LiDAR

Buffered Tiling to Avoid Edge Artefacts

TL;DR: Query the index for every tile intersecting the target buffered by the neighbourhood your processing needs, read them all, crop to the buffered extent, process, then crop back to the nominal extent before writing. The buffer exists to be discarded.

# Context and Motivation

This guide is part of Tile Indexing, Buffering and Merging. Any operation that looks at a point’s neighbours — ground classification, outlier removal, height above ground, rasterization — behaves differently at a tile edge, because half the neighbourhood is in another file. Buffered tiling is the general fix, and the same pattern serves every one of those stages.

The width of the buffer is not a matter of taste. It is the reach of the widest neighbourhood operation in the pipeline: the search radius for a rasterizer, the window size for a morphological ground filter, the neighbour distance for an outlier filter. Take the largest, add a margin, and use that. Too narrow leaves a seam; too wide costs read time and nothing else, which makes erring generous the cheap mistake.

Three extents, and only one of them is written Concentric rectangles around one tile. The outer rectangle is what gets read — the tile plus its neighbours out to the buffer distance. The middle is what gets processed, so that edge cells have full neighbourhoods. The inner rectangle is the nominal tile extent, and it is the only part written to the output. written processed read read: tile + neighbours process: nominal + buffer write: nominal extent only every buffer point is read, informs the edge, and is then discarded

# Prerequisites and Assumptions

Requirement Detail
A tile index from building a tile index
PDAL 2.4+ with filters.merge and filters.crop
Read access to neighbours each worker must be able to open the tiles around its own
A known neighbourhood reach the largest radius or window in the pipeline

# Step-by-Step Implementation

# Step 1 — Determine the buffer from the pipeline

Rasterizer radius 1.4 m, SMRF window 33 m, outlier neighbour distance about 2 m — the largest is 33, so a buffer of 40 m is safe and cheap.

# Step 2 — Ask the index for neighbours

sql
SELECT b.location FROM tiles a, tiles b
WHERE a.location LIKE '%tile_0431%'
  AND ST_Intersects(ST_Buffer(a.geom, 40), b.geom)
  AND b.location <> a.location

# Step 3 — Read the tile and its neighbours together

json
{"pipeline": ["tile_0431.laz", "tile_0430.laz", "tile_0432.laz",
              {"type": "filters.merge"},
              {"type": "filters.crop", "bounds": "([511960, 513040], [4782960, 4784040])"}]}

# Step 4 — Process as usual

Classification, filtering and rasterization all now see full neighbourhoods at the tile edge.

# Step 5 — Crop back before writing

json
{"type": "filters.crop", "bounds": "([512000, 513000], [4783000, 4784000])"}

Skip this and every point in the buffer is written twice — once by this tile and once by its neighbour.

# Complete Working Example

python
"""Process one tile with a neighbourhood buffer, writing only its own extent."""
from __future__ import annotations

import json
import logging
import subprocess
from pathlib import Path

import pdal

LOG = logging.getLogger("buffered")


def neighbours(index: Path, tile: str, buffer_m: float, layer: str = "tiles") -> list[str]:
    sql = (f"SELECT b.location FROM {layer} a, {layer} b "
           f"WHERE a.location LIKE '%{tile}%' "
           f"AND ST_Intersects(ST_Buffer(a.geom, {buffer_m}), b.geom) "
           f"AND b.location <> a.location")
    out = subprocess.run(["ogr2ogr", "-f", "CSV", "/vsistdout/", str(index),
                          "-dialect", "SQLITE", "-sql", sql],
                         check=True, capture_output=True, text=True).stdout
    return [ln.strip() for ln in out.splitlines()[1:] if ln.strip()]


def process(tile: Path, index: Path, out: Path,
            bounds: tuple[float, float, float, float], buffer_m: float = 40.0) -> int:
    xmin, ymin, xmax, ymax = bounds
    bx0, by0, bx1, by1 = xmin - buffer_m, ymin - buffer_m, xmax + buffer_m, ymax + buffer_m

    others = neighbours(index, tile.stem, buffer_m)
    LOG.info("%s: reading %d neighbour(s)", tile.name, len(others))

    stages: list = [{"type": "readers.las", "filename": str(tile)}]
    stages += [{"type": "readers.las", "filename": p} for p in others]
    stages += [
        {"type": "filters.merge"},
        {"type": "filters.crop", "bounds": f"([{bx0}, {bx1}], [{by0}, {by1}])"},
        {"type": "filters.outlier", "method": "statistical", "mean_k": 12, "multiplier": 2.5},
        {"type": "filters.range", "limits": "Classification![7:7]"},
        {"type": "filters.smrf", "window": 33, "slope": 0.2, "threshold": 0.6, "cell": 1.0},
        # Crop back BEFORE writing, or the buffer is duplicated into the output.
        {"type": "filters.crop", "bounds": f"([{xmin}, {xmax}], [{ymin}, {ymax}])"},
        {"type": "writers.las", "filename": str(out), "compression": "laszip",
         "forward": "all"},
    ]
    return pdal.Pipeline(json.dumps({"pipeline": stages})).execute()


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    n = process(Path("tiles/tile_0431.laz"), Path("index/tiles.gpkg"),
                Path("out/tile_0431.laz"),
                bounds=(512000.0, 4783000.0, 513000.0, 4784000.0))
    LOG.info("wrote %d points", n)
What the neighbour query returns A three by three block with the target tile in the centre. A forty metre buffer around it intersects all eight surrounding tiles, so the read list is nine files rather than one. At a corner of the campaign the same query returns four, and the pipeline needs no special case for it. target the 40 m buffer read list: 9 files write extent: 1 tile at a campaign corner the query returns four files, and the pipeline needs no special case

# Key Parameter Table

Choice Value Why
buffer width ≥ largest neighbourhood reach Ground filter window usually dominates
first crop buffered extent Bounds the work; the neighbours may be much larger than needed
second crop nominal extent Prevents duplicate points across tile outputs
neighbour query buffer the target only One spatial predicate instead of thousands
read cost buffer area ÷ tile area A 40 m buffer on a 1 km tile is about 16% more points

# Verification

Output extent equals the nominal extent. pdal info --summary on the result; a bounding box larger than the tile means the second crop is missing.

Total output points equal total input points. Summed across the campaign, with overlap removed. More means duplication; far fewer means a crop is too tight.

Edges are seamless. The seam test from building a seamless DTM mosaic applies to the raster products.

Choosing the buffer width Four buffer widths on a one kilometre tile. With no buffer the seam error is eleven centimetres and no extra points are read. At ten metres the error falls to two centimetres for four percent more points. At forty metres it is below a millimetre for sixteen percent more. At a hundred metres nothing further is gained and the read cost is forty-four percent. 1 km tile, SMRF window 33 m no buffer seam error 11 cm 10 m 2 cm · 4% more points read 40 m below 1 mm · 16% more points read 100 m no better · 44% more past the largest neighbourhood the curve flattens and you pay for nothing

# Gotchas and Edge Cases

Forgetting the second crop. Every buffered point is written twice, and the merged campaign has duplicate returns in every overlap strip.

Buffering the crop without reading the neighbours. Adds empty space, not context. The read list must include the neighbours.

Edge tiles have no neighbours on one side. Expected at the campaign boundary. The seam test should exclude the outer edge.

Classification differing across the seam. If each tile classifies independently, the ground surface can still disagree slightly in the overlap. Buffering the classification, as here, is what prevents it.

# Frequently Asked Questions

How wide should the buffer be?

At least the largest neighbourhood reach in the pipeline — the rasterizer search radius, the ground filter window, the outlier neighbour distance, whichever is greatest. A 33 metre SMRF window makes 40 metres a safe choice, which costs about sixteen percent more points read on a one kilometre tile.

What happens if I forget the final crop?

Every point in the buffer is written twice, once by the tile that owns it and once by its neighbour. The merged campaign then has duplicate returns in every overlap strip, which inflates density metrics and produces doubled returns in anything computed downstream.

Is a wider buffer ever harmful?

Only to your read budget. Beyond the largest neighbourhood in the pipeline the seam error stops improving while the extra points keep costing time, so erring generous is the cheap mistake and erring narrow is the expensive one.

What about tiles at the edge of the campaign?

They genuinely have no neighbours on one side, and that is not a defect. The seam test should exclude the outer boundary, and any coverage report should distinguish “no data collected” from “processing lost it”.

Related in Tile Indexing, Buffering and Merging

↑ Tile Indexing, Buffering and Merging overview --- title: "Tile Indexing, Buffering and Merging" description: "Building a queryable tile index from LiDAR headers, using it to find... Merging Processed Tiles into One LAZ --- title: "Merging Processed Tiles into One LAZ" description: "Crop each tile to its nominal extent before merging or the overlap strips... Building a Tile Index with pdal tindex --- title: "Building a Tile Index with pdal tindex" description: "The one-line command that turns a directory of LAZ files into a queryable...