Streaming LAS Tiles to Avoid Out-of-Memory Kills

You finish a dense reconstruction, kick off a post-processing script to filter or reproject the point cloud, and the process vanishes — no traceback, just Killed in the shell and a line in dmesg about the OOM killer reaping your Python interpreter. The culprit is almost always a single line like las = laspy.read("cloud.laz"), which pulls a multi-gigabyte dense cloud fully into RAM before you have touched a point. A 4 cm corridor survey easily runs to 500 million points; at roughly 30 bytes per point once numpy expands the record, that is 15 GB of resident memory for the raw XYZ alone, before any working copy. This page shows how to process clouds far larger than RAM by streaming them in bounded chunks and tiling them spatially, so peak memory stays flat regardless of file size.

Why reading the whole file into RAM kills the process

A LAS/LAZ file stores points as fixed-width records, and the convenient laspy.read() call decompresses and materialises every record into numpy arrays in one shot. The OOM kill happens because that materialisation is not the only copy: reprojecting the coordinates allocates a second array, a boolean filter mask allocates a third, and the filtered output allocates a fourth. Peak memory is a multiple of the already-large raw cloud, and the Linux OOM killer terminates the process the instant that peak exceeds available RAM — without a Python exception you can catch, because the kill is a SIGKILL from the kernel, not a MemoryError.

The fix is to never hold the whole cloud at once. Two complementary strategies achieve that. Chunked streaming reads a fixed number of points at a time, processes them, writes them out, and releases them before reading the next chunk, so peak memory is set by the chunk size, not the file size. Spatial tiling splits the cloud into bounded ground tiles that each fit comfortably in RAM, which is essential when the operation needs spatial locality — gridding to a DSM, ground classification, or neighbourhood filtering — because a purely sequential chunk stream has no spatial coherence. Both are backpressure mechanisms: they trade RAM for disk bandwidth, which is exactly the trade the broader memory management for large point clouds workflow is built around.

Resident memory over time: whole-file read versus tiled streaming A time series of resident memory for the same forty gigabyte cloud processed two ways. The whole-file read climbs steadily from the start, crosses the host memory limit part-way through, and terminates at a marked kill point where the kernel out-of-memory killer sends signal nine. The tiled streaming run instead shows a low sawtooth: memory rises for each tile and drops back as the tile is released, never approaching the limit, and continues past the point where the other run died through to completion. host memory limit elapsed time resident memory SIGKILL (9) — OOM killer laspy.read() chunk_iterator() — one tile at a time run completes Peak memory of the streaming run is set by tile size, not by file size — the two lines would not converge on a larger file.

Figure 1 — The failure is not gradual. A whole-file read climbs monotonically until the kernel intervenes, and everything computed up to that moment is lost; the streaming run’s peak is a property of the tile, so the same code survives a file ten times larger.

Minimal reproducible solution

laspy 2.x exposes chunk_iterator, which yields fixed-size blocks of points and pairs with an appending writer so you never hold more than one chunk plus one output buffer. The routine below reprojects and filters a multi-GB LAZ to a new file with flat memory, whatever the input size.

import laspy
import numpy as np


def stream_filter_laz(src: str, dst: str, z_min: float, z_max: float,
                      points_per_chunk: int = 2_000_000) -> int:
    """Copy src -> dst keeping only points in [z_min, z_max], streaming.

    Peak memory is bounded by points_per_chunk, not the file size, because
    each chunk is read, filtered, written, and released before the next.
    """
    kept = 0
    with laspy.open(src) as reader:
        # Reuse the source header so scale, offset, and CRS VLRs are preserved.
        with laspy.open(dst, mode="w", header=reader.header) as writer:
            for chunk in reader.chunk_iterator(points_per_chunk):
                z = chunk.z                       # decoded float z for this chunk only
                mask = (z >= z_min) & (z <= z_max)
                if mask.any():
                    writer.append_points(chunk[mask])
                    kept += int(mask.sum())
                # chunk and mask go out of scope here and are freed before
                # the next iteration reads the following block.
    return kept


n = stream_filter_laz("dense.laz", "filtered.laz", z_min=180.0, z_max=340.0)
print(f"wrote {n} points with flat memory")

The two lines that matter are laspy.open(src) instead of laspy.read(src), and chunk_iterator, which returns a generator of ScaleAwarePointRecord blocks. Because the writer appends, the output is built incrementally on disk and never fully resident either. When the operation needs spatial neighbourhoods rather than a sequential pass, delegate the tiling to PDAL, whose filters.splitter cuts the cloud into a grid of bounded tiles you then process one at a time.

import json
import subprocess
from pathlib import Path


def tile_with_pdal(src: str, out_dir: str, tile_m: float = 250.0) -> list[Path]:
    """Split a large cloud into square ground tiles with PDAL streaming mode.

    filters.splitter runs in PDAL's streaming pipeline, so the full cloud is
    never resident; each tile is written as it fills.
    """
    Path(out_dir).mkdir(parents=True, exist_ok=True)
    pipeline = {
        "pipeline": [
            src,
            {"type": "filters.splitter", "length": tile_m},
            {"type": "writers.las",
             "filename": f"{out_dir}/tile_#.laz",     # # expands to a tile index
             "compression": "laszip",
             "forward": "all"},                        # carry header/CRS forward
        ]
    }
    subprocess.run(["pdal", "pipeline", "--stdin", "--stream"],
                   input=json.dumps(pipeline), text=True, check=True)
    return sorted(Path(out_dir).glob("tile_*.laz"))


tiles = tile_with_pdal("dense.laz", "tiles", tile_m=250.0)
for t in tiles:
    n = stream_filter_laz(str(t), str(t.with_suffix(".filt.laz")), 180.0, 340.0)
    print(f"{t.name}: kept {n} points")
Why a streamed tile needs a halo before it is filtered Two tiles of a point cloud shown side by side. On the left, a tile processed with hard edges: a neighbourhood search near the boundary finds only points inside the tile, so the ground filter sees a truncated neighbourhood and misclassifies points along the seam, drawn as a broken line of wrong classifications down the edge. On the right, the same tile read with a surrounding halo of extra points: the neighbourhood search is complete, classification is correct across the seam, and the halo points are discarded after filtering so no point is written twice. Hard tile edges k-NN search neighbourhood truncated at the seam ground filter misclassifies a stripe of points Tile read with a halo written tile halo neighbourhood complete across the seam halo points are filtered, then dropped before writing Halo width must exceed the largest neighbourhood radius any downstream filter uses — otherwise the seam artefact simply moves inward.

Figure 2 — Streaming is not only a memory technique: a tile read without a halo silently changes the answer along every seam, because neighbourhood operators cannot see the points they need.

Edge-case matrix

These are the input conditions that turn a “just read it” script into an OOM kill, the symptom each produces, and how the streaming approach handles it.

Input variant Downstream symptom if unchecked Expected handling
300M+ point dense LAZ laspy.read OOM-kills before first point processed laspy.open + chunk_iterator; never materialise the whole cloud
Reproject + filter in one pass 3–4× RAM multiple from intermediate arrays Transform inside the chunk loop so only one chunk’s copies exist
Neighbourhood filter / gridding Sequential chunks lack spatial locality; wrong results Tile spatially with PDAL filters.splitter, process per tile
.laz without lazrs backend laspy raises on compressed read Install laspy[lazrs]; confirm before streaming
Chunk size set too large Peak memory still exceeds RAM, OOM returns Lower points_per_chunk until peak RSS fits with headroom
Output header not forwarded Tiles lose CRS/scale, misalign downstream Reuse source header; PDAL forward: all
Many tiny tiles Per-tile overhead dominates, throughput collapses Size tiles so each holds a few million points, not thousands
Where an out-of-memory kill leaves no traceback A comparison of two termination paths. On the left, a Python MemoryError path: the allocation fails inside the interpreter, an exception propagates, the except block runs, a traceback is written to the log and the process exits with a code the orchestrator can read. On the right, the kernel out-of-memory killer path: the kernel sends signal nine directly to the process, so no exception is raised, no except block runs, nothing is flushed to the log, and the orchestrator sees only exit status 137. A note explains that this is why the guard must be a pre-flight check on available memory rather than a try-and-catch. MemoryError (catchable) allocation fails in the interpreter except block runs, log is flushed orchestrator retries the tile SIGKILL (not catchable) kernel picks the largest RSS signal 9 — no exception, no flush partial output left on disk orchestrator sees only exit 137 A try/except cannot defend the right-hand path. The guard has to run before the allocation: check available memory, size the tile to fit, and refuse to start otherwise.

Figure 3 — The two ways a memory-bound stage ends, and why only one of them can be handled in Python. Exit code 137 is the shell’s way of reporting signal 9, and it is the signature to look for when a job vanishes without a traceback.

Verification snippet

Confirm two things after streaming: that peak memory actually stayed bounded, and that no points were silently lost — the sum of tile point counts must equal the filtered total. Measuring peak RSS is what turns “I think it streamed” into proof.

import laspy
import psutil
import os


def verify_stream(src: str, tiles: list, expected_kept: int) -> None:
    peak_mb = psutil.Process(os.getpid()).memory_info().rss / (1024 ** 2)
    # Point coverage: tiles must account for every kept point.
    tile_total = 0
    for t in tiles:
        with laspy.open(str(t)) as r:
            tile_total += r.header.point_count
    assert tile_total == expected_kept, (
        f"point loss: tiles={tile_total} expected={expected_kept}")
    # Memory: peak must stay well under a full-file read would have needed.
    src_points = laspy.open(src).header.point_count
    assert peak_mb < 0.1 * src_points * 30 / 1024**2 + 2048, (
        f"peak {peak_mb:.0f} MB too high — lower points_per_chunk")
    print(f"OK: {tile_total} points, peak {peak_mb:.0f} MB")

The coverage assertion is the one that catches a silent bug: an off-by-one in a chunk mask or a tile boundary that drops points shows up as a mismatch here long before it corrupts a downstream DSM. Log peak RSS per run so a regression that reintroduces a full-file read is caught immediately.

When to escalate

This page keeps post-processing of an existing cloud inside memory bounds. Escalate when the constraint is elsewhere:

  • RAM blows up during dense matching, before a LAS file even exists. That is an alignment-stage memory problem, not a post-processing one — cap workers and stream descriptors per memory management for large point clouds and the dense-matching guidance in reducing RAM usage during dense matching.
  • Streaming is correct but you then need a rasterised surface. Tiling is the right input to gridding — feed the tiles straight into generating DSM and DTM from point clouds with PDAL rather than re-reading the whole cloud.
  • The process is killed but this streaming code is not in the traceback path. Confirm the OOM origin by sampling peak RSS around each stage; the kill may be in a concurrent worker whose memory ceiling was never capped.

← Memory Management for Large Point Clouds