Resolving Memory Errors in PDAL and laspy

The job has run every week for a year. This week the process disappears after eleven minutes with no traceback, no error in the log, and an exit status of 137. Nothing in the code changed; the survey was larger.

Exit 137 is the kernel’s out-of-memory killer, and it is the most common failure in point-cloud processing at scale. It is also entirely predictable — the memory a stage needs is a function of the point count, and checking it before starting converts an unexplained kill into a message naming the problem.

This page covers sizing a job, identifying which stage is responsible, and the two remedies that actually work. It is the operational companion to choosing a library for streaming versus in-memory work.

What actually consumes the memory

A point in memory is not a point on disk. On disk, a LAZ point with colour is around 8 bytes after compression. In memory, PDAL holds it as a PointView entry with every dimension materialised as a fixed-width field — X, Y, Z as doubles, classification, intensity, RGB, plus any computed dimensions — which is 40 to 80 bytes depending on the point format and what has been added.

That is a factor of five to ten between file size and memory footprint, and it is where the surprise comes from. A 3 GB LAZ is not a 3 GB memory requirement; it is closer to 20 GB.

Neighbourhood filters then add their own structures on top. A KD-tree over n points costs roughly another 24 to 32 bytes per point, and several filters build one. A statistical outlier filter over 250 million points therefore wants around 20 GB for the points and another 7 GB for the tree, before any temporaries.

From compressed file size to peak memory for a neighbourhood filter A four-step escalation for a survey of two hundred and fifty million points. The compressed LAZ file is three gigabytes. Materialised in memory with all dimensions it becomes twenty gigabytes. Adding a KD-tree for a neighbourhood filter brings it to twenty-seven gigabytes. Temporaries during the filter peak at about thirty-four gigabytes. A marker shows a typical thirty-two gigabyte workstation, which the final step exceeds. A note states that the factor between file size and peak memory is commonly ten or more. 3 GB LAZ on disk 20 GB materialised 27 GB + KD-tree 34 GB + temporaries 32 GB workstation A factor of ten between file size and peak memory is normal, not pathological.

Figure 1 — Where a three-gigabyte file becomes a thirty-four-gigabyte job.

Minimal reproducible solution

Estimate before running, and refuse with a useful message rather than being killed.

import laspy
import psutil

BYTES_PER_POINT = {6: 44, 7: 56, 8: 60}      # by LAS point format, materialised
TREE_BYTES_PER_POINT = 28
TEMPORARY_FACTOR = 1.25


def estimate_peak_bytes(path: str, *, builds_tree: bool = True) -> dict:
    """Peak memory a blocking stage will need for this file."""
    with laspy.open(path) as fh:
        h = fh.header
        n = int(h.point_count)
        per_point = BYTES_PER_POINT.get(h.point_format.id, 48)

    points = n * per_point
    tree = n * TREE_BYTES_PER_POINT if builds_tree else 0
    peak = int((points + tree) * TEMPORARY_FACTOR)
    return {"points": n, "points_bytes": points, "tree_bytes": tree,
            "peak_bytes": peak, "peak_gb": peak / 1e9}


def require_headroom(path: str, *, builds_tree: bool = True,
                     margin: float = 1.2) -> dict:
    """Fail before starting, with the numbers, rather than being killed."""
    est = estimate_peak_bytes(path, builds_tree=builds_tree)
    available = psutil.virtual_memory().available
    if est["peak_bytes"] * margin > available:
        raise MemoryError(
            f"{est['points']:,} points need about {est['peak_gb']:.1f} GB and "
            f"{available / 1e9:.1f} GB is available — tile the input "
            f"(suggested {max(2, int(est['peak_bytes'] * margin / available) + 1)} tiles)")
    return est

Suggesting the tile count in the error message is a small touch that saves the next person a calculation, and it makes the failure self-documenting in a log.

The two remedies

Stream, where the stage allows it. A pipeline of streamable stages runs in bounded memory regardless of input size. This is always the first thing to check, because it costs nothing when it applies.

Tile with a halo, where it does not. Split the cloud spatially, process each tile with a buffer larger than the stage’s neighbourhood reach, and trim the buffer before merging.

import json
import math
import subprocess
from pathlib import Path


def tiles_needed(path: str, *, target_gb: float = 8.0) -> int:
    """How many tiles keep each one inside a memory budget."""
    est = estimate_peak_bytes(path)
    return max(1, math.ceil(est["peak_gb"] / target_gb))


def tiled_run(src: str, out_dir: str, stage: dict, *, reach_m: float,
              target_gb: float = 8.0, extent_m: float = 2000.0) -> None:
    """Tile, run a blocking stage per tile with a halo, write trimmed tiles."""
    n = tiles_needed(src, target_gb=target_gb)
    side = extent_m / math.sqrt(n)
    halo = max(reach_m * 2.0, 10.0)
    if side <= halo * 3:
        raise ValueError(f"tiles of {side:.0f} m are too small for a {halo:.0f} m "
                         "halo — the buffer would dominate; use a larger machine")

    Path(out_dir).mkdir(parents=True, exist_ok=True)
    pipeline = {"pipeline": [
        src,
        {"type": "filters.splitter", "length": side, "buffer": halo},
        stage,
        {"type": "writers.las", "filename": f"{out_dir}/tile_#.laz",
         "compression": "laszip", "forward": "all"},
    ]}
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)

The guard against tiles that are too small matters. As tiles shrink, the halo becomes a larger fraction of each one, and at some point most of the work is redundant buffer processing — at which point renting a larger machine for an hour is cheaper than the compute being wasted.

Reducing memory by changing the shape of the work, not the machine Two columns. The works column lists streaming the pipeline so memory stays bounded, cropping to a tile before any in-memory stage, selecting only the attributes actually needed rather than every field, and processing by tile with a halo where a neighbourhood operation requires context. The does not work column lists adding more memory, which moves the threshold without removing it, raising the swap size, which converts an error into an unbounded slowdown, and reducing the point count by random thinning, which changes the answer rather than the cost. changes the shape of the work stream the pipeline: memory stays bounded crop to a tile before any in-memory stage select only the attributes needed tile with a halo for neighbourhood work only moves the threshold more memory — the same failure, later more swap — an error becomes an unbounded slowdown random thinning — changes the answer, not the cost A memory error is a signal that the work is shaped wrong, not that the machine is small.

Figure 3 — Four fixes and three non-fixes.

Edge-case matrix

Symptom Likely cause Remedy
Exit 137, no traceback Kernel OOM killer Estimate first; tile or stream
std::bad_alloc Allocation failed inside PDAL Same causes, different reporting
MemoryError in Python laspy read of a whole file Use the chunk iterator
Slow with heavy disk activity Swapping, not yet killed Reduce concurrency; swapping is worse than tiling
Fails only on some inputs Those are larger Size from the file, not from experience
Works alone, fails under a scheduler Concurrent jobs sharing memory Admission control on total, not per-job
Memory grows across a loop Arrays retained between iterations Delete and collect explicitly between tiles
COPC write fails at the end The octree build needs everything Not tileable; size the machine

The scheduler row is worth designing for rather than discovering. Two jobs that each fit comfortably will kill each other when run together, and the right control is a shared budget rather than a per-job limit — the same admission-control pattern used in orchestrating photogrammetry jobs with Python schedulers.

Verification snippet

import gc
import resource


def peak_rss_gb() -> float:
    """Peak resident memory of this process so far, in gigabytes."""
    return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024 / 1024


def run_and_record(fn, estimate_gb: float, *, tolerance: float = 1.5) -> dict:
    """Run a stage and compare observed peak memory against the estimate."""
    gc.collect()
    before = peak_rss_gb()
    fn()
    after = peak_rss_gb()
    observed = after - before

    note = "estimate is sound"
    if observed > estimate_gb * tolerance:
        note = (f"used {observed:.1f} GB against an estimate of {estimate_gb:.1f} GB — "
                "raise the per-point figure for this stage")
    elif observed < estimate_gb / tolerance:
        note = "estimate is conservative; the budget could be tightened"
    return {"observed_gb": observed, "estimate_gb": estimate_gb, "note": note}

Feeding the observed figure back into the per-point constants is what keeps the estimator useful. The numbers in this page are reasonable defaults; the numbers for a specific pipeline, with its specific extra dimensions, are better, and they cost one run to measure.

Redundant halo work as tiles shrink A curve of the share of processed points that fall in the halo rather than in the tile proper, against tile side length, for a forty metre halo. At a one thousand metre tile the halo is about sixteen percent of the work. At five hundred metres it is thirty-four percent. At two hundred metres it is ninety-six percent, and at one hundred metres the halo exceeds the tile itself. A shaded region beyond sixty percent marks where renting a larger machine costs less than the wasted compute. a bigger machine is cheaper 100 m 200 500 1000 2000 m tile side length work in the halo 96 % halo at 200 m 16 % at 1 km — comfortable Tiling is not free, and past a point it is more expensive than the memory it saves.

Figure 2 — The limit of tiling as a strategy, for a 40 m halo.

Reducing the footprint before adding hardware

Two changes often remove the problem entirely and cost nothing.

Drop dimensions you are not using. A cloud carried through a pipeline with intensity, GPS time, scan angle, four computed features and RGB costs far more per point than one carrying coordinates and classification. A filters.ferry that keeps only what the stage needs can halve the footprint of a neighbourhood filter, and the full set can be restored on the final write from the source file.

Decimate for the stages that do not need full density. A ground classifier does not benefit from 400 points per square metre; it works just as well on a 10 cm-spaced subsample and then transfers the classification back to the full cloud by nearest neighbour. That is a tenfold reduction in the expensive stage with no measurable change in the result, and it is standard practice in lidar workflows that has not fully crossed into photogrammetric ones.

Both changes are reversible and both are cheaper than tiling, so they belong ahead of it in the order of remedies.

When to escalate

  • The stage cannot be tiled and does not fit. COPC writing and global sorting genuinely need everything. Size the machine for the largest expected survey rather than engineering around it.
  • Memory grows steadily across a loop that should be bounded. That is a retention bug in the loop, not a sizing problem. Delete arrays explicitly and collect between iterations before adding hardware.
  • The estimate is right and the budget is still exceeded under load. Concurrency is the variable. Move to a shared admission budget rather than reducing every job’s footprint.

← Troubleshooting Point Cloud Processing Failures