Tile Indexing, Buffering and Merging

A LiDAR campaign arrives as thousands of files and no map. Somewhere in that directory are the four tiles covering the road corridor a client asked about, and finding them by opening headers one at a time is the kind of task that quietly consumes an afternoon and then has to be repeated next week. A tile index solves it once: a vector layer with one polygon per file, carrying the file path and whatever else is worth knowing, queryable by any GIS or by ogr2ogr in a shell script. Every workflow in this section — buffered processing, spatial fan-out, merging — starts by asking the index a question. This topic belongs to Batch Automation and Cloud Integration for PDAL.

The index is also where two things that are easy to conflate get separated. A tile’s nominal extent is the rectangle it owns in the delivery grid. Its actual extent is the bounding box of the points inside it, which is usually a little larger because of flight-line overlap and sometimes a lot larger because someone merged two tiles and forgot. Buffered processing needs the first; correctness checks need the second; and a delivery in which they differ systematically is telling you something before you have processed a single point.

# Prerequisites

# Core Workflow Architecture

1. Enumerate the files. A glob, a manifest, or an object listing. Freeze it — an index built while files are still arriving describes a moment that has passed.
2. Read each header. pdal tindex opens the header only: bounds, point count, CRS, version. On a thousand tiles this is seconds, not minutes.
3. Build a polygon per file. By default the header bounding box, which is the actual extent. The nominal grid extent, if you have it, is worth carrying as separate fields.
4. Attach attributes. Path, point count, CRS, acquisition date, processing status — whatever later queries will filter on.
5. Write a vector layer. GeoPackage for anything with more than a handful of fields.
6. Query it. Spatial queries for coverage, attribute queries for status, and self-joins for neighbour lookup.

# Full Implementation

"""Build a tile index with attributes, and query it for neighbours."""
from __future__ import annotations

import json
import logging
import subprocess
from pathlib import Path

import pdal

LOG = logging.getLogger("tindex")


def header_info(path: Path) -> dict:
    """Header-only read: bounds, count and CRS without touching the points."""
    p = pdal.Pipeline(json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(path), "count": 1}]}))
    p.execute()
    qi = p.quickinfo["readers.las"]
    b = qi["bounds"]
    return {
        "path": str(path),
        "points": int(qi["num_points"]),
        "minx": b["minx"], "maxx": b["maxx"],
        "miny": b["miny"], "maxy": b["maxy"],
        "srs": (qi.get("srs") or {}).get("horizontal", "")[:64],
    }


def build_index(tiles: list[Path], out: Path, layer: str = "tiles") -> int:
    """Write a GeoPackage tile index using pdal tindex."""
    out.parent.mkdir(parents=True, exist_ok=True)
    if out.exists():
        out.unlink()  # tindex appends; a stale index is worse than none

    cmd = ["pdal", "tindex", "create", str(out),
           "--lyr_name", layer, "-f", "GPKG", "--fast_boundary"]
    cmd += [str(t) for t in tiles]
    subprocess.run(cmd, check=True)
    LOG.info("indexed %d tiles into %s", len(tiles), out)
    return len(tiles)


def neighbours(index: Path, tile_name: str, buffer_m: float = 10.0,
               layer: str = "tiles") -> list[str]:
    """Every tile whose polygon intersects a buffered target tile."""
    sql = (
        f"SELECT b.location FROM {layer} a, {layer} b "
        f"WHERE a.location LIKE '%{tile_name}%' "
        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 [line.strip() for line in out.splitlines()[1:] if line.strip()]


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    tiles = sorted(Path("tiles").glob("*.laz"))
    build_index(tiles, Path("index/tiles.gpkg"))
    LOG.info("neighbours of tile_0431: %s", neighbours(Path("index/tiles.gpkg"), "tile_0431"))

# Code Breakdown

The index is deleted before rebuilding. pdal tindex create appends to an existing layer, so re-running it against a stale file silently doubles every polygon. Deleting first makes the operation idempotent, which matters more than it sounds when the command sits inside a nightly job.

--fast_boundary uses the header bounding box. The alternative computes a concave hull from the points, which is far more informative for irregular coverage and far slower. Use the fast form for a working index and the slow form once, for a coverage map.

The neighbour query buffers the target, not every tile. Buffering one polygon and testing intersection is a single spatial predicate. Buffering all of them first would be correct and much slower on a large index.

quickinfo never reads points. The count: 1 reader plus quickinfo is a header read; it is what makes indexing a thousand tiles take seconds.

Paths go in as given. Absolute paths break when storage moves; relative paths break when the working directory changes. Pick one, write it down, and make the index-building script the only thing that decides.

# Parameter Reference Table

# Validation and Integrity Checks

Every file appears exactly once. A count of index features against a count of input files catches both the double-append and a glob that missed a subdirectory.

Every tile has a CRS, and they agree. A campaign whose index contains two coordinate systems has a real problem that will otherwise surface at merge time — see fixing CRS mismatches.

Coverage has no holes. Dissolve the polygons and compare the result against the project boundary. A hole is a missing delivery.

Actual extents do not wildly exceed nominal ones. A tile whose bounding box is twice its grid cell usually contains points that belong to a different tile, and every downstream buffered read will be wrong about what it needs.

# The Index as a Job Ledger

The obvious use of a tile index is spatial lookup. The more valuable one, on a campaign that takes days to process, is as the record of what has been done.

Add three attributes and the index becomes a work queue: status, run_id and updated_at. A worker claims a tile by writing running with its own run identifier, writes done when the output object is durable, and leaves failed with the error when it is not. Restarting a job then becomes a query — every tile whose status is not done — rather than a decision about which files to reprocess, and the answer is the same whether the interruption was one crashed worker or a whole cluster going away.

This is the same idempotency argument that appears in scaling PDAL tile processing with AWS Batch, with the state in a queryable layer instead of in marker objects. The two are complements rather than alternatives: the marker object is the authoritative per-tile fact, because it sits next to the output and cannot disagree with it, while the index is the aggregate view that answers “how far through are we” without listing a bucket. Where they disagree, the marker wins and the index is stale.

Two cautions. A GeoPackage is a SQLite file and does not tolerate many concurrent writers, so on a large fan-out the workers should write markers and a single collector should update the index, rather than every worker opening it. And the index must never become the only record of the run: it is a convenience built from things that are true elsewhere, and rebuilding it from the outputs should always be possible.

# Nominal Grids and Delivery Reality

A delivery grid is a promise about where tiles begin and end, and the point cloud inside a tile is only approximately bound by it. Three discrepancies show up often enough to plan for.

Overlap. Adjacent tiles frequently share a strip of points, either because the delivery was cut on flight lines rather than on the grid, or because the supplier buffered them deliberately. Merging tiles that overlap double-counts every point in the strip, which inflates density metrics and produces duplicate returns in any classification that follows. Cropping each tile to its nominal extent before merging costs one stage and removes the problem entirely.

Undershoot. A coastal or boundary tile may contain points across only part of its nominal cell, with the rest genuinely outside the project. Treating the nominal extent as the processing extent then produces a large NoData region that looks like a failure and is not. Carrying both extents in the index lets a report distinguish “no data was collected here” from “processing lost it”.

Gross mismatch. A tile whose bounding box is several times its grid cell almost always holds points that belong to other tiles, usually because two deliveries were concatenated. This is worth catching at index time, because every buffered read afterwards will fetch the wrong neighbours — the tile’s own bounding box already covers them, so the neighbour query returns nothing and the edge effects reappear silently.

The general rule is to store both extents, use the nominal one for partitioning work and the actual one for validation, and treat a systematic difference between them as information about the delivery rather than as noise to be normalised away.

# Performance Tuning

Index headers, not points. --fast_boundary is the difference between seconds and hours on a large campaign, and for the buffered-read use case the bounding box is all that is needed.

Build once per delivery, not per job. The index is an artefact with a lifetime. Rebuild when files change; otherwise read it.

Use a spatial index on the layer. GeoPackage creates one by default; a shapefile needs ogrinfo -sql "CREATE SPATIAL INDEX ON tiles". Without it a neighbour query on 20,000 tiles is a full scan per lookup.

Keep the index next to the data. An index in one bucket and tiles in another produces path fields that are correct for exactly one consumer.

# Common Errors and Troubleshooting

Every tile appears twice. tindex create appended to an existing file. Delete before rebuilding.

Neighbour queries return nothing. Either the layer has no spatial index and the SQL dialect fell back, or the buffer distance is in different units from the geometry — a lat/long index buffered by 10 “metres” buffers by ten degrees.

Field names are truncated. Shapefile’s ten-character limit. Write GeoPackage.

Paths in the index do not resolve. Absolute paths from a machine that no longer exists. Store paths relative to the index and resolve them at read time.

A tile is in the index but unreadable. The index records what the header said at build time. Add a readable attribute set by a periodic verification pass rather than assuming.

# Frequently Asked Questions

Why build a tile index instead of globbing the directory?

Because the questions are spatial. Which tiles cover this corridor, which tiles neighbour this one, which tiles are missing from the coverage — none of them can be answered by a filename pattern, and answering them by opening headers each time turns a one-second query into a several-minute one.

Why does my index have every tile twice?

Because pdal tindex create appends to an existing layer. Re-running it against a file that already exists silently doubles every polygon, which is why the index-building script should delete the output first and be idempotent by construction.

What is the difference between nominal and actual tile extent?

The nominal extent is the rectangle the tile owns in the delivery grid; the actual extent is the bounding box of the points inside it. They differ because of flight-line overlap, and a tile whose actual extent is far larger than its nominal one usually holds points belonging to a neighbour.

Should the index store absolute or relative paths?

Relative to the index, resolved at read time. Absolute paths break when storage moves, which it always does; relative paths break only when someone separates the index from the data, which is a mistake worth making visible.

# Related

• Batch Automation and Cloud Integration for PDAL — the section this workflow belongs to
• Building a Tile Index with pdal tindex — the command, its options and the checks to run on the result
• Buffered Tiling to Avoid Edge Artefacts — using the index to give every tile its neighbours
• Merging Processed Tiles into One LAZ — putting the results back together without duplicating the overlaps
• Building a Seamless DTM Mosaic from Tiles — the raster-side equivalent of the same problem
• Scaling PDAL Tile Processing with AWS Batch — the fan-out that consumes this index