Python LiDAR

Retiling Flightline Files into a Grid

TL;DR: Run pdal tile "swaths/*.laz" "tiles/tile_#.laz" --length 1000 --origin_x 0 --origin_y 0 to stream every flightline into 1 km tiles aligned to round coordinates. Points from overlapping swaths land in the same tile automatically, tiles are named by grid index, and memory stays low because the tile application streams. For custom logic inside a pipeline, filters.splitter does the same split in memory.

# Context and Motivation

This guide is part of Tile Indexing and Merging. Raw or lightly processed LiDAR often arrives as flightlines: long, narrow swaths, one file per pass, each overlapping its neighbours. Swaths are awkward for everything downstream. A single swath may be tens of kilometres long and too big to process in one job; any location is split across two or three files; and neither a DTM grid nor a tile index lines up with swath boundaries.

Retiling into a regular grid solves this once: each tile contains every point from every swath inside its square, tiles are of predictable size, and tile names encode their location. Every later stage — ground classification, rasterising, array jobs — then works per tile.

From swaths to tiles On the left, two long diagonal flightline swaths overlap each other, with the centreline of a third pass dashed between them. On the right, the same area is cut into a regular grid of square tiles; each tile contains points from whichever swaths cross it, and overlap areas are merged into the same tile. overlapping flightlines 0_21_22_23_2 0_11_12_13_1 0_01_02_03_0 square tiles named by grid index

# Prerequisites and Assumptions

  • PDAL 2.x with the tile application.
  • Flightlines in a projected CRS with metre or foot units; reproject first if they are in geographic coordinates.
  • Enough disk for the output, roughly the size of the input.

# Step-by-Step Implementation

# Step 1 — Choose tile size

Pick a length that keeps a tile at a comfortable point count for your processing, typically 500 m to 1.5 km for 8–20 points per square metre.

# Step 2 — Fix the grid origin

Set --origin_x and --origin_y to round coordinates (0, 0 works for a grid aligned to multiples of the tile length), so tile edges fall on round numbers and future deliveries share the same grid.

# Step 3 — Run pdal tile

Give an input glob and an output pattern containing #; PDAL replaces # with the tile’s column and row.

# Step 4 — Rename to coordinates if needed

Convert index-based names to lower-left coordinate names such as 571000_4190000.laz, which read more clearly and survive re-origining.

# Step 5 — Index the result

Build a tile index of the new tiles with pdal tindex for later queries.

# Complete Working Example

Command line, streaming all swaths:

bash
mkdir -p tiles
pdal tile "swaths/*.laz" "tiles/tile_#.laz" \
  --length 1000 --origin_x 0 --origin_y 0 \
  --buffer 0
ls tiles | head        # tile_571_4190.laz  tile_571_4191.laz ...

With origin 0 and a 1,000 m length, grid index equals the lower-left coordinate divided by 1,000, so tile_571_4190.laz covers x 571,000–572,000 and y 4,190,000–4,191,000.

The same split inside a Python pipeline, with coordinate names and a point-count report:

python
import glob
import json

import pdal

LENGTH = 1000.0
readers = [{"type": "readers.las", "filename": f, "tag": f"r{i}"}
           for i, f in enumerate(sorted(glob.glob("swaths/*.laz")))]
spec = {"pipeline": readers + [
    {"type": "filters.merge", "inputs": [r["tag"] for r in readers]},
    {"type": "filters.splitter", "length": LENGTH, "origin_x": 0, "origin_y": 0},
]}
p = pdal.Pipeline(json.dumps(spec))
p.execute()

for arr in p.arrays:                     # one array per tile
    x0 = int(arr["X"].min() // LENGTH * LENGTH)
    y0 = int(arr["Y"].min() // LENGTH * LENGTH)
    out = f"tiles/{x0}_{y0}.laz"
    w = pdal.Pipeline(json.dumps([{"type": "writers.las", "filename": out,
                                    "compression": "laszip", "forward": "all",
                                    "extra_dims": "all"}]), arrays=[arr])
    w.execute()
    print(out, len(arr))

The Python version holds all points in memory, so reserve it for modest areas; the pdal tile application streams and is the right tool for full deliveries.

Streaming, not loading Swath files are read point by point in chunks. Each point's tile index is computed from its X and Y, and the point is appended to the open writer for that tile. Memory holds only one chunk plus the open writers, not the whole delivery. swath chunk(e.g. 10k points) col = ⌊(x − x₀) / L⌋row = ⌊(y − y₀) / L⌋ writer tile_571_4190 writer tile_571_4191 writer tile_572_4190 memory: one chunk plus open writers

# Keeping Flightline Identity

Merging swaths into tiles loses the file boundary, but not necessarily the information. The PointSourceId field, set by most acquisition software to the flightline number, travels with every point, so a tile can still be split back by flightline for swath-to-swath accuracy checks or intensity normalisation, as in normalizing intensity across flightlines. Check before retiling that PointSourceId is populated and distinct per swath; if it is not, assign it from the file order with filters.assign on each reader before merging, because it cannot be recovered afterwards. The same applies to the overlap flag and GpsTime, which later tools use to identify overlapping returns.

Tile edges also create a small classic problem: ground filters and rasterisers behave poorly at tile boundaries. Retile without a buffer for storage, and add buffers at processing time as in buffered tiling to avoid edge artefacts.

# Key Parameter Table

Option Typical Notes
--length 500–1500 m Tile side in CRS units
--origin_x, --origin_y 0, 0 Aligns tiles to round coordinates
--buffer 0 for storage Add buffers at processing time instead
--out_srs optional Reproject while tiling
output pattern tiles/tile_#.laz # becomes column_row
filters.splitter length same as above In-pipeline, in-memory alternative

# Verification

  • Point conservation. The sum of points across tiles equals the sum across swaths; compare with pdal info --summary totals.
  • Bounds on the grid. Each tile’s bounds fall within its named square, with minima at or above the tile corner.
  • No slivers. Tiles along the edge of the delivery may be partly empty; a very small tile at a corner is normal, but many tiny tiles suggest a wrong origin.

# Gotchas and Edge Cases

Too many open files. Streaming retiling keeps a writer open per active tile. A delivery spanning thousands of tiles in one swath can hit the operating system’s file limit; raise it with ulimit -n, or tile in bands.

Feet and metres. In a State Plane CRS in US survey feet, --length 1000 makes 1,000-foot tiles. Choose lengths in the CRS’s own units.

Length is in CRS units Two squares compare a 1,000-unit tile in a metre-based CRS, one kilometre on a side, with a 1,000-unit tile in a US survey foot CRS, about 305 metres on a side and roughly one eleventh of the area. The same length option gives very different tiles. 1,000 m 1,000 US ft ≈ 305 m about 1/11 of the area --length is read in the CRS's own units: check the axis unit before tiling

Header CRS mismatch. All swaths must share a CRS. A single swath in a different zone lands in wildly wrong tiles; check every input’s CRS first.

# Frequently Asked Questions

How do I split LiDAR flightlines into square tiles?

Use the pdal tile application with an input glob, an output pattern containing a hash sign, and a tile length. It streams the input, so large deliveries need little memory, and overlapping swaths are merged into shared tiles.

What tile size should I use for LiDAR?

A size that keeps tiles at a manageable point count for your processing, often 500 metres to 1.5 kilometres for typical aerial densities. Denser data suits smaller tiles.

Why set a tile origin?

An origin of zero aligns tile edges to round multiples of the tile length, so names map to coordinates and future deliveries in the same area fall onto the same grid.

Does retiling lose which flightline a point came from?

No, as long as PointSourceId is populated. It travels with each point and can be used later to separate flightlines within a tile.

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... Buffered Tiling to Avoid Edge Artefacts --- title: "Buffered Tiling to Avoid Edge Artefacts" description: "Read the tile plus its neighbours, process the buffered extent, write... 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... Making Tile Outputs Idempotent --- title: "Making Tile Outputs Idempotent" description: "Design per-tile LiDAR outputs so reruns are safe: deterministic output keys… 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...