Python LiDAR

Building a Tile Index with pdal tindex

TL;DR: pdal tindex create index.gpkg -f GPKG --lyr_name tiles --fast_boundary tiles/*.laz — then immediately check that the feature count equals the file count, that one CRS is present, and that no polygon is much larger than a grid cell.

# Context and Motivation

This guide is part of Tile Indexing, Buffering and Merging. The command is one line; the value is in what you do with the result in the next five minutes, because a tile index is the cheapest opportunity you will get to find out what is actually in a delivery.

pdal tindex reads headers rather than points, so indexing a thousand tiles takes seconds. Each file becomes one feature whose geometry is either the header bounding box — fast, and enough for buffered reads — or a computed boundary that follows irregular coverage, which is slower and far more informative when the flight lines do not fill their tiles.

Fast boundary against computed boundary One tile whose points cover only a diagonal swath. The header bounding box is the full rectangle, which overstates coverage but is instant to obtain. The computed boundary follows the swath edge, showing the true extent, and costs a full pass over the points to produce. --fast_boundary computed boundary header bbox — instant, overstates coverage true extent — one pass over the points fast form for a working index; slow form once, for a coverage map

# Prerequisites and Assumptions

Requirement Detail
PDAL 2.4+ with the tindex application
GDAL/OGR for the output driver and for querying afterwards
One CRS tiles in mixed coordinate systems produce a meaningless layer
Write access tindex create appends, so the target must not already exist

# Step-by-Step Implementation

# Step 1 — Build the index

bash
rm -f index/tiles.gpkg
pdal tindex create index/tiles.gpkg -f GPKG --lyr_name tiles \
  --fast_boundary tiles/*.laz

The rm is not defensive tidiness — create appends, so without it a rebuild doubles every feature.

# Step 2 — Count features against files

bash
ls tiles/*.laz | wc -l
ogrinfo -so index/tiles.gpkg tiles | grep "Feature Count"

# Step 3 — Check the coordinate systems agree

bash
ogrinfo -so index/tiles.gpkg tiles | grep -A3 "Layer SRS"

# Step 4 — Look for oversized polygons

A feature whose area is much larger than the nominal grid cell holds points that belong somewhere else.

# Step 5 — Add the attributes you will filter on later

tindex writes the path and geometry; point counts, acquisition dates and processing status are worth joining in from a header pass.

# Complete Working Example

python
"""Build a tile index and immediately audit it."""
from __future__ import annotations

import json
import logging
import subprocess
from pathlib import Path

LOG = logging.getLogger("tindex_build")


def build(tiles: list[Path], out: Path, layer: str = "tiles") -> None:
    out.parent.mkdir(parents=True, exist_ok=True)
    if out.exists():
        out.unlink()
    subprocess.run(["pdal", "tindex", "create", str(out), "-f", "GPKG",
                    "--lyr_name", layer, "--fast_boundary",
                    *[str(t) for t in tiles]], check=True)


def audit(index: Path, expected: int, cell_m: float, layer: str = "tiles") -> dict:
    info = subprocess.run(["ogrinfo", "-so", str(index), layer],
                          check=True, capture_output=True, text=True).stdout
    count = int([l for l in info.splitlines() if "Feature Count" in l][0].split(":")[1])

    sql = (f"SELECT COUNT(*) AS n FROM {layer} "
           f"WHERE (ST_MaxX(geom)-ST_MinX(geom)) > {cell_m * 1.5}")
    oversized = subprocess.run(
        ["ogrinfo", "-q", "-dialect", "SQLITE", "-sql", sql, str(index)],
        check=True, capture_output=True, text=True).stdout

    result = {
        "features": count,
        "files": expected,
        "duplicated": count > expected,
        "oversized_report": oversized.strip().splitlines()[-1:] or ["n = 0"],
    }
    if count != expected:
        raise AssertionError(
            f"index holds {count} features for {expected} files — "
            "the index was appended to, or the glob missed a subdirectory"
        )
    return result


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    tiles = sorted(Path("tiles").glob("*.laz"))
    build(tiles, Path("index/tiles.gpkg"))
    print(json.dumps(audit(Path("index/tiles.gpkg"), len(tiles), cell_m=1000.0), indent=2))

# Key Parameter Table

Option Effect
-f GPKG Output driver; GeoPackage avoids shapefile’s ten-character field names
--lyr_name Layer name, needed by every later SQL query
--fast_boundary Header bounding box instead of a computed hull
-t_srs Reproject the index geometry only; the tiles are untouched
--stdin Read the file list from standard input, for very large campaigns
--write_absolute_path Store absolute rather than relative paths; usually the wrong choice
What to put in the index besides the geometry Eight fields worth carrying on every index feature: the file path, the point count, the coordinate system, the acquisition date, the nominal extent identifier, the processing status, the run identifier that last touched it, and the timestamp of that change. The first three come from the header; the rest turn the index into a work ledger. location path, relative to the index num_points from the header, not a count of rows srs_wkt so a mixed campaign is visible at a glance acquired delivery date, for provenance status pending · running · done · failed pdal tindex writes the first field and the geometry; everything else is a header pass you run once

# Verification

Feature count equals file count. Asserted in the example; catches both the double-append and a glob that missed a directory.

One CRS. Two in an index means the campaign cannot be merged as it stands.

No wildly oversized polygons. A tile whose bounding box spans several grid cells will break every neighbour query that follows.

Paths resolve. Open two at random and confirm they exist from the directory the index will be read from.

What the boundary choice costs Indexing 1,240 tiles two ways. With fast boundaries the run reads headers only and completes in eleven seconds. With computed boundaries it reads every point in every file and takes just over four hours. The output differs only in the shape of the polygons. 1,240 tiles, 186 GB total --fast_boundary 11 s — headers only computed boundary 4 h 12 m — every point in every file both produce a usable index; only one can be re-run when a delivery changes, which makes the fast form the default and the slow form occasional.

# Gotchas and Edge Cases

create appends. The most common tile-index bug, and it produces an index that works for lookups and reports twice the coverage.

Relative paths are relative to the working directory, not the index. Build the index from the directory it will be read from, or resolve paths explicitly at read time.

Shapefile truncates field names. Ten characters, silently — the DBF limit that also bites tile-attribute exports.

A file that fails to open is skipped. tindex reports it and continues, so the feature count check is what tells you a tile is unreadable.

The index describes the moment it was built. It records what each header said at build time, which is exactly what makes it fast and exactly why it goes stale. Rebuild it whenever the delivery changes rather than patching individual features, and keep the build script in version control beside the pipelines that read it — a tile index nobody can regenerate is a liability the first time somebody asks whether it is still correct.

# Frequently Asked Questions

Why does my index have twice as many features as files?

Because pdal tindex create appends to an existing layer rather than replacing it. The resulting index works for lookups and reports double the coverage, which is why the build script should delete the output first.

Should I use fast boundaries or computed ones?

Fast for anything you rebuild. On 1,240 tiles the header-only form takes eleven seconds and the computed form takes over four hours, because the latter reads every point. Compute boundaries once, for a coverage map you will show someone.

Why GeoPackage rather than shapefile?

Field names. Shapefile truncates them to ten characters silently, so acquisition_date and acquisition_sensor become the same field. GeoPackage also carries a spatial index by default, which neighbour queries on a large campaign need.

What happens to a file that cannot be opened?

It is reported and skipped, and the index is built from the rest. That is reasonable behaviour and it is why comparing the feature count with the file count matters — the skipped file is otherwise invisible.

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... 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...