Chunked Reading of Large LAS Files with laspy

TL;DR: with laspy.open(path) as f: for chunk in f.chunk_iterator(2_000_000): ... yields point records of at most two million points each, with the same fields and scaled accessors as a full read. Accumulate statistics across chunks, or write filtered chunks to an open laspy.open(out, mode="w", header=...) writer. Memory stays at roughly one chunk regardless of file size.

# Context and Motivation

This guide is part of laspy and NumPy Workflows for LAS Data. Statewide deliveries, merged project files and dense corridor surveys routinely produce single LAZ files of hundreds of millions of points — far more than a laptop, and often more than a batch worker, can hold. laspy.read on such a file fails with a memory error or, worse, slows the machine to a crawl as it swaps. The chunk iterator keeps the convenience of laspy’s NumPy fields while bounding memory, which covers most “look at every point once” jobs: class histograms, extent and density checks, filtering to a smaller file, splitting by flightline or class, and fixing an attribute in place.

# Prerequisites and Assumptions

• laspy 2.x with lazrs for LAZ (pip install "laspy[lazrs]").
• Work that visits each point once and does not need neighbours across chunk boundaries.
• Enough disk for outputs, since results are written as you go.

# Step-by-Step Implementation

# Step 1 — Open and read the header

laspy.open(path) returns a reader; its header gives count, format and bounds for sizing accumulators.

# Step 2 — Choose a chunk size

One to five million points is a good range: large enough that per-chunk Python overhead is negligible, small enough to keep memory to a few hundred megabytes.

# Step 3 — Accumulate, do not collect

Update counters, histograms and grids inside the loop. Appending chunks to a list rebuilds the whole file in memory and defeats the purpose.

# Step 4 — Write as you go

For outputs, open one writer before the loop with the header you want and call write_points per chunk.

# Step 5 — Split with several writers

To split by an attribute (flightline, class), keep a dictionary of writers keyed by value and route each chunk’s subsets to the right one.

# Complete Working Example

Three jobs in one pass over a large file — a class histogram, a 10 m density grid, and a split into one file per flightline — all in bounded memory:

"""One pass over a large LAZ: class histogram, density grid, split by PointSourceId."""
from __future__ import annotations

from contextlib import ExitStack
from pathlib import Path

import laspy
import numpy as np

SRC = Path("big/county_block_7.laz")
OUT = Path("out/by_line")
CHUNK = 2_000_000
CELL = 10.0

OUT.mkdir(parents=True, exist_ok=True)

with laspy.open(SRC, laz_backend=laspy.LazBackend.LazrsParallel) as reader, ExitStack() as stack:
    h = reader.header
    x0, y0 = h.mins[0], h.mins[1]
    cols = int(np.ceil((h.maxs[0] - x0) / CELL)) + 1
    rows = int(np.ceil((h.maxs[1] - y0) / CELL)) + 1
    density = np.zeros((rows, cols), dtype=np.int32)
    classes = np.zeros(256, dtype=np.int64)
    writers: dict[int, laspy.LasWriter] = {}

    for i, chunk in enumerate(reader.chunk_iterator(CHUNK)):
        classes += np.bincount(chunk.classification, minlength=256)
        c = ((np.asarray(chunk.x) - x0) // CELL).astype(int)
        r = ((np.asarray(chunk.y) - y0) // CELL).astype(int)
        np.add.at(density, (r, c), 1)

        for psid in np.unique(chunk.point_source_id):
            if psid not in writers:
                writers[psid] = stack.enter_context(
                    laspy.open(OUT / f"line_{psid}.laz", mode="w", header=h))
            writers[psid].write_points(chunk[chunk.point_source_id == psid])

        if i % 20 == 0:
            print(f"chunk {i}: {(i + 1) * CHUNK:,} points processed")

print({int(k): int(v) for k, v in enumerate(classes) if v})
print(f"density: median {np.median(density[density > 0]) / CELL**2:.1f} pts/m²")
print(f"{len(writers)} flightline files written")

ExitStack closes every flightline writer when the block ends, which is when laspy finalizes each file’s header counts and bounds.

# Key Parameter Table

# Verification

• Counts add up. The class histogram’s total must equal header.point_count, and the per-flightline files’ counts must sum to it too.
• Headers valid. Open each output with laspy.open and check point_count and bounds are non-zero and plausible; a writer that was never closed leaves an invalid header.
• Flat memory. Watch RSS during a run; it should plateau after the first few chunks.

total = sum(laspy.open(p).header.point_count for p in OUT.glob("line_*.laz"))
assert total == laspy.open(SRC).header.point_count

# Gotchas and Edge Cases

Chunks are not spatial. Points arrive in file order, which usually follows acquisition time. A chunk is a stripe of scan lines, not a tile; anything that needs neighbours must be done with a spatial tool, not per chunk.

Many writers, many open files. Splitting by an attribute with thousands of values opens thousands of files. Group values, or write to intermediate files per group and merge later.

Append mode and LAZ. mode="a" appends to uncompressed LAS but is not available for LAZ. Keep one writer open across the loop rather than reopening per chunk.

Modifying fields in place. Chunks are copies; assigning to chunk.classification does not change the source file. To fix an attribute, write corrected chunks to a new file.

# Frequently Asked Questions

How do I read a huge LAS file with laspy?

Open it with laspy.open and loop over chunk_iterator with a chunk size of a few million points. Each chunk behaves like a small point record with the usual fields, and memory stays bounded by the chunk size.

Can I write a filtered file while iterating?

Yes. Open a writer before the loop with the header you want, and call write_points with the filtered subset of each chunk. laspy updates the header counts and bounds when the writer closes.

What chunk size should I use?

One to five million points works well. Smaller chunks increase overhead from Python and decompression calls; larger chunks increase memory without much speed gain.

How do I speed up reading LAZ?

Use the parallel lazrs backend, which decompresses on multiple cores. For repeated reads of the same data, convert once to uncompressed LAS on fast local disk, or to COPC if you need spatial queries.

# Related

• laspy and NumPy Workflows for LAS Data — the library and its data model
• Reading LAS into NumPy with laspy — the whole-file version
• Writing a LAS File from NumPy Arrays — headers for new files
• Iterating PDAL Arrays in Chunks from Python — the PDAL equivalent
• LAZ vs Uncompressed LAS for Iterative Processing — when decompression cost matters