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.

# 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")

# 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.

# 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.

# Related

• Memory Management for Large Point Clouds
• Reducing RAM Usage During Dense Matching
• Generating DSM and DTM from Point Clouds with PDAL

← Memory Management for Large Point Clouds