Estimating PDAL Memory from Point Layout

TL;DR: Memory in standard mode ≈ point count × bytes per point × safety factor. Get the count from the header (pdal info --summary), get bytes per point by reading one point through the pipeline’s own stages and checking arrays[0].dtype.itemsize, add 8 bytes for every double a later filter creates, and multiply by 1.5–2 for neighbourhood stages that build a k-d tree.

# Context and Motivation

This guide is part of Memory Management in PDAL. Out-of-memory failures are the most common way PDAL batch jobs die, and they are expensive: a job killed by the kernel after forty minutes of SMRF has wasted forty minutes and gives no Python traceback. The fix is to know, before launching, how much memory each tile needs — then choose instance sizes, tile sizes and concurrency to fit.

PDAL makes this predictable. In standard mode it holds every point in a table whose row width is the sum of the dimension sizes, so memory scales linearly with point count and with the number of dimensions. Everything else — spatial indexes, raster buffers, Python copies — is a multiple on top of that base.

# Prerequisites and Assumptions

• PDAL 2.x and the Python bindings.
• The pipeline you intend to run, and one representative input tile.
• Standard (non-streaming) execution. Streaming mode bounds memory by chunk size instead, which is covered in running a PDAL pipeline in streaming mode.

# Step-by-Step Implementation

# Step 1 — Read the point count without loading points

pdal info --summary tile.laz reads only the header and returns summary.num_points. It is instant even for multi-gigabyte files.

# Step 2 — Measure bytes per point empirically

Run the pipeline’s own stages on the first point only, using filters.head with count: 1, and read arrays[0].dtype.itemsize. This counts every dimension the reader produces and every dimension the filters add — no guessing about extra bytes or PDRF layouts.

# Step 3 — Add stage overhead

Neighbourhood stages (filters.smrf, filters.outlier, filters.hag_nn, filters.covariancefeatures) build a spatial index: budget roughly 0.5 to 1× the point table again. writers.gdal holds its raster buffers: rows × columns × 8 bytes × number of output types.

# Step 4 — Add the Python copy

pipeline.arrays returns NumPy arrays that share or copy PDAL’s data depending on version and usage. If your code touches arrays, budget one extra copy of the table.

# Step 5 — Apply a safety factor and decide

Multiply by 1.2 for allocator slack and compare against the worker’s memory, leaving room for the OS and anything else on the machine.

# Complete Working Example

"""Estimate peak memory for a PDAL pipeline on a given tile."""
from __future__ import annotations

import copy
import json
import subprocess
from pathlib import Path

import pdal

NEIGHBOURHOOD = {"filters.smrf", "filters.pmf", "filters.outlier", "filters.hag_nn",
                 "filters.hag_delaunay", "filters.covariancefeatures", "filters.normal",
                 "filters.cluster", "filters.dbscan", "filters.csf"}


def point_count(path: Path) -> int:
    out = subprocess.run(["pdal", "info", "--summary", str(path)],
                         capture_output=True, text=True, check=True).stdout
    return int(json.loads(out)["summary"]["num_points"])


def bytes_per_point(spec: dict) -> int:
    probe = copy.deepcopy(spec)
    stages = [s for s in probe["pipeline"] if not str(s.get("type", "")).startswith("writers.")]
    stages.insert(1, {"type": "filters.head", "count": 1})
    p = pdal.Pipeline(json.dumps({"pipeline": stages}))
    p.execute()
    return p.arrays[0].dtype.itemsize


def raster_bytes(spec: dict, extent_m: float = 1000.0) -> int:
    total = 0
    for s in spec["pipeline"]:
        if s.get("type") == "writers.gdal":
            cells = (extent_m / float(s.get("resolution", 1.0))) ** 2
            kinds = len(str(s.get("output_type", "all")).split(","))
            total += int(cells * 8 * (6 if s.get("output_type", "all") == "all" else kinds))
    return total


def estimate(spec: dict, tile: Path, python_copy: bool = True) -> dict:
    n = point_count(tile)
    bpp = bytes_per_point(spec)
    table = n * bpp
    types = {s.get("type") for s in spec["pipeline"] if isinstance(s, dict)}
    index = table * (1.0 if types & NEIGHBOURHOOD else 0.0)
    extra = table if python_copy else 0
    peak = 1.2 * (table + index + extra + raster_bytes(spec))
    return {"points": n, "bytes_per_point": bpp,
            "table_gb": round(table / 1e9, 2), "peak_gb": round(peak / 1e9, 2)}


if __name__ == "__main__":
    spec = json.loads(Path("dtm.json").read_text())
    print(estimate(spec, Path("tiles/t_0431.laz")))

For a 1 km² tile with 42 million points, a PDRF 6 input and a pipeline that adds HeightAboveGround, the output looks like:

{'points': 42000000, 'bytes_per_point': 56, 'table_gb': 2.35, 'peak_gb': 8.54}

# Key Parameter Table

# Verification

The estimate is only useful if it matches reality. Measure peak memory for a few tiles with the method in measuring peak memory of a PDAL pipeline and compare: a good estimate is within 20 percent and errs high. If measured peaks are consistently lower, reduce the index factor; if higher, look for a stage the estimate does not know about.

# Gotchas and Edge Cases

Point density varies across tiles. Tiles over forest or with flightline overlap can hold twice the points of tiles over fields. Estimate from the largest tile in the batch, not the median, or size workers per tile.

Extra bytes dimensions. A vendor file with ten extra-byte dimensions carries them into PDAL’s table. The empirical itemsize catches this; an estimate from the PDRF alone would not.

Merged inputs. Pipelines that merge neighbours hold all of them at once; sum the counts of every reader.

Streaming changes everything. If every stage streams, memory is bounded by chunk_size × bytes per point, independent of tile size. Check pipeline.streamable before sizing for standard mode.

# Frequently Asked Questions

How many bytes per point does PDAL use in memory?

It depends on the dimensions present. A typical point format 6 tile uses around 48 bytes per point in PDAL’s table before any filters add dimensions, because coordinates are stored as 8-byte doubles. Measure it for your data with a one-point probe.

Why is PDAL’s memory use much larger than the LAZ file size?

LAZ compresses points to a few bytes each, and even uncompressed LAS stores coordinates as scaled 4-byte integers. In memory, PDAL expands coordinates to doubles and adds any computed dimensions, so the table is often five to ten times the LAZ size.

How much extra memory does SMRF or outlier filtering need?

Neighbourhood filters build a spatial index over all points, which typically costs half to one times the point table again for the duration of the stage. Budget for it whenever such a stage is present.

Can I avoid the Python copy of the points?

Yes, by not accessing pipeline.arrays when you only need the pipeline to write files. Execute and let the writer produce the output; read the result later only if you need it.

# Related

• Memory Management in PDAL — the execution model behind these numbers
• Measuring Peak Memory of a PDAL Pipeline — checking the estimate
• Diagnosing PDAL Out-of-Memory Failures — when the estimate was wrong
• LAZ vs Uncompressed LAS for Iterative Processing — on-disk versus in-memory size
• Understanding LAS Point Data Record Formats — the formats being expanded