Orchestrating PDAL Workflows with Airflow

TL;DR: Define a discover_tiles @task that returns a list of S3 keys, define a process_tile @task that runs a PDAL SMRF-plus-DTM pipeline for one key, then wire them with process_tile.expand(tile_key=discover_tiles()) — Airflow fans out one retry-protected task instance per tile and passes each result path through XCom.

# Context and Motivation

This guide is part of Airflow DAG Orchestration, which explains the design patterns; here we build one concrete DAG end to end and run it. If you have ever kicked off a for tile in tiles: subprocess.run(...) loop and then watched it die on tile 340 with no way to resume, this is the fix. Airflow turns that loop into a graph where every tile is an independently scheduled, independently retryable unit of work, and where a single failed tile does not discard the 339 that already succeeded.

The workload is deliberately narrow: given a prefix full of LAZ tiles in object storage, produce one 1-metre DTM GeoTIFF per tile. Each tile runs the same two-part PDAL pipeline — classify ground with SMRF, then interpolate a raster with writers.gdal — the same operations covered in depth by SMRF Ground Classification and DTM Raster Generation. What Airflow adds is the scheduling, retry, and result-tracking scaffolding around that pipeline.

# Prerequisites and Assumptions

This walkthrough assumes tiles are in EPSG:26918 (NAD83 / UTM 18N) — a common US survey CRS distinct from the WGS84 UTM codes used elsewhere. Adjust the spatialreference value to match your data.

# Step-by-Step Implementation

# Step 1 — Discover the tiles

The discovery task lists LAZ keys and returns them as a list. That return value is what mapping iterates over, so it must be a genuine list.

@task
def discover_tiles(bucket: str, prefix: str) -> list[str]:
    import boto3
    s3 = boto3.client("s3")
    pages = s3.get_paginator("list_objects_v2").paginate(Bucket=bucket, Prefix=prefix)
    keys = [o["Key"] for p in pages for o in p.get("Contents", [])
            if o["Key"].endswith(".laz")]
    if not keys:
        raise ValueError(f"no tiles under s3://{bucket}/{prefix}")
    return keys

Note the import boto3 inside the function. Keeping heavy imports out of module scope is what stops the scheduler from paying that cost on every DAG parse.

# Step 2 — Process one tile with PDAL

This task is written for a single tile. Dynamic mapping will call it once per key; the function itself never loops.

@task(pool="pdal_pool", retries=3)
def process_tile(tile_key: str, bucket: str, out_prefix: str) -> dict:
    import json, posixpath, pdal
    name = posixpath.basename(tile_key).replace(".laz", ".tif")
    out_key = f"{out_prefix}{name}"
    pipe = {"pipeline": [
        {"type": "readers.las",
         "filename": f"/vsis3/{bucket}/{tile_key}",
         "spatialreference": "EPSG:26918"},
        {"type": "filters.smrf", "slope": 0.2, "window": 16.0, "threshold": 0.45},
        {"type": "filters.range", "limits": "Classification[2:2]"},
        {"type": "writers.gdal",
         "filename": f"/vsis3/{bucket}/{out_key}",
         "resolution": 1.0, "output_type": "idw", "nodata": -9999},
    ]}
    count = pdal.Pipeline(json.dumps(pipe)).execute()
    if count == 0:
        raise RuntimeError(f"{name}: no ground points classified")
    return {"tile": tile_key, "raster": out_key, "ground_points": count}

# Step 3 — Expand over the tile list

Mapping is one line. .expand() binds the mapped argument; constant arguments use .partial().

tiles = discover_tiles(BUCKET, RAW_PREFIX)
results = process_tile.partial(bucket=BUCKET, out_prefix=DTM_PREFIX).expand(tile_key=tiles)

Airflow reads the tiles list at runtime and creates process_tile[0], process_tile[1], … one per key. Each is scheduled, retried, and logged on its own.

# Step 4 — Reduce the results

A plain (non-mapped) task receives the full list of mapped return values. It only runs once every mapped instance succeeds, making it the natural place for batch-level assertions.

@task
def summarise(results: list[dict]) -> dict:
    total = sum(r["ground_points"] for r in results)
    return {"tiles": len(results), "ground_points": total}

# Step 5 — Trigger and watch it run

Unpause the DAG in the UI (or airflow dags trigger pdal_tile_walkthrough), then open the grid view. Each mapped square is one tile; hover to see its state, click to read its log and XCom. A green summarise means the whole batch is done.

# Complete Working Example

Save as dags/pdal_tile_walkthrough.py. It is self-contained and runnable against any S3 prefix of LAZ tiles.

"""pdal_tile_walkthrough.py — a per-tile PDAL DTM workflow in Airflow 2.x."""
from __future__ import annotations

import json
import posixpath
from datetime import datetime, timedelta

from airflow.decorators import dag, task

BUCKET = "lidar-demo"
RAW_PREFIX = "survey-a/raw/"
DTM_PREFIX = "survey-a/dtm/"

default_args = {"retries": 3, "retry_delay": timedelta(minutes=3),
                "execution_timeout": timedelta(minutes=30)}


@dag(dag_id="pdal_tile_walkthrough", start_date=datetime(2024, 7, 1),
     schedule=None, catchup=False, default_args=default_args,
     max_active_tasks=8, tags=["pdal", "dtm"])
def pdal_tile_walkthrough():

    @task
    def discover_tiles(bucket: str, prefix: str) -> list[str]:
        import boto3
        s3 = boto3.client("s3")
        pages = s3.get_paginator("list_objects_v2").paginate(
            Bucket=bucket, Prefix=prefix)
        keys = [o["Key"] for p in pages for o in p.get("Contents", [])
                if o["Key"].endswith(".laz")]
        if not keys:
            raise ValueError(f"no tiles under s3://{bucket}/{prefix}")
        return keys

    @task(pool="pdal_pool")
    def process_tile(tile_key: str, bucket: str, out_prefix: str) -> dict:
        import pdal
        name = posixpath.basename(tile_key).replace(".laz", ".tif")
        out_key = f"{out_prefix}{name}"
        pipe = {"pipeline": [
            {"type": "readers.las",
             "filename": f"/vsis3/{bucket}/{tile_key}",
             "spatialreference": "EPSG:26918"},
            {"type": "filters.smrf", "slope": 0.2,
             "window": 16.0, "threshold": 0.45},
            {"type": "filters.range", "limits": "Classification[2:2]"},
            {"type": "writers.gdal",
             "filename": f"/vsis3/{bucket}/{out_key}",
             "resolution": 1.0, "output_type": "idw", "nodata": -9999},
        ]}
        count = pdal.Pipeline(json.dumps(pipe)).execute()
        if count == 0:
            raise RuntimeError(f"{name}: no ground points classified")
        return {"tile": tile_key, "raster": out_key, "ground_points": count}

    @task
    def summarise(results: list[dict]) -> dict:
        total = sum(r["ground_points"] for r in results)
        empty = [r["tile"] for r in results if r["ground_points"] == 0]
        if empty:
            raise ValueError(f"empty tiles: {empty}")
        print(f"{len(results)} DTMs, {total:,} ground points")
        return {"tiles": len(results), "ground_points": total}

    tiles = discover_tiles(BUCKET, RAW_PREFIX)
    results = process_tile.partial(
        bucket=BUCKET, out_prefix=DTM_PREFIX).expand(tile_key=tiles)
    summarise(results)


dag = pdal_tile_walkthrough()

# Key Parameter Table

# Verification

After the run, confirm three things. First, the grid view: every process_tile square should be green, and the count of squares should equal the number of LAZ tiles in the prefix. Second, the XCom of summarise should report a tiles count matching that number and a plausible ground_points total. Third, list the output prefix to confirm one GeoTIFF per tile:

aws s3 ls s3://lidar-demo/survey-a/dtm/ | grep -c '\.tif$'

Spot-check one raster’s validity with pdal info on its source, or open it in QGIS and confirm the elevation range is sane for the terrain. A DTM whose values are all -9999 means every cell was nodata — usually a sign the SMRF stage classified nothing, which the count == 0 guard should already have caught as a task failure.

# Gotchas and Edge Cases

1. Mapping over a non-list collapses the fan-out. .expand() needs a real list. If discover_tiles returns a generator, a dict, or a single string, Airflow either errors or produces one oddly-shaped instance. Always materialise the keys with a list comprehension, as the example does, and raise on empty results so an empty prefix fails loudly instead of expanding to zero instances.

2. The worker must actually have PDAL. A DAG that parses fine can still fail every task with ModuleNotFoundError: No module named 'pdal' because the scheduler parsed it but the worker lacks the library. PDAL, GDAL, and PROJ must live in the worker’s environment. If your workers are heterogeneous, run the pipeline through PDAL Docker Containers via DockerOperator so the runtime travels with the task instead of depending on each host.

3. Top-level code throttles the scheduler. Anything outside a task body runs on every DAG parse. A boto3.client("s3").list_objects_v2(...) at module scope will hammer S3 and slow the scheduler for every DAG in the deployment. Keep the module body to declarations; push all I/O and heavy imports inside @task functions. This is the single most common performance mistake when porting a script to Airflow.

4. Mapped instance count has a ceiling. Dynamic mapping is capped at max_map_length (default 1024). A national dataset with 40,000 tiles will hit that wall. Batch the tile list into groups of, say, 200 and map process_batch over the groups, looping the tiles inside each task — trading some per-tile granularity for a mapping count Airflow can track comfortably. The same fan-out logic then scales through AWS Batch Processing when the compute exceeds one cluster.

# Frequently Asked Questions

How many tiles can Airflow dynamic task mapping expand to?

The default ceiling is 1024 mapped instances per task, set by the max_map_length configuration value. You can raise it, but very large expansions inflate the scheduler’s bookkeeping; for tens of thousands of tiles, batch the tile list into groups and map over the groups instead of individual files.

Where do I install PDAL so the Airflow worker can import it?

PDAL and its Python bindings must exist in the same environment the worker executes tasks in. Install it into the worker image with conda (conda install -c conda-forge pdal python-pdal), or run the task through DockerOperator against an image that already contains PDAL so the worker itself stays thin.

Why is my Airflow scheduler slow after adding a PDAL DAG?

The scheduler re-parses every DAG file on a short interval. Any expensive code at module level — a boto3 client call, an S3 listing, or a heavy import outside a task — runs on every parse and throttles the scheduler. Move all such work inside @task functions so it only runs when a task executes.

Can I add a reprojection or validation step to each tile?

Yes. Insert extra PDAL stages into the pipeline list inside process_tile — for example a filters.reprojection before writers.gdal, following the PDAL stage chaining order. Everything stays within one task, so the mapping and retry behaviour is unchanged.

How do I re-run only the tiles that failed?

In the grid view, select the failed process_tile instances and clear them; Airflow re-queues just those instances while leaving the successful ones untouched. That selective replay is the core reason to orchestrate with Airflow rather than a shell loop.

# Related

• Airflow DAG Orchestration — parent guide to modelling PDAL workflows as DAGs
• Batch Automation and Cloud Integration for PDAL — running PDAL beyond a single workstation
• SMRF Ground Classification — the classification stage each tile runs
• DTM Raster Generation — the writers.gdal step that produces each DTM
• PDAL Docker Containers — package PDAL so any worker can run the task