Memory Management for Large Point Clouds

When a drone survey grows from a single test flight into a corridor or regional dataset, the dense point cloud stops fitting in RAM and naive Open3D or numpy pipelines begin to thrash, swap, or trigger kernel out-of-memory kills mid-reconstruction. This page shows surveying technicians and Python GIS developers how to process UAV-derived point clouds that are far larger than physical memory by streaming LAS/LAZ files in bounded chunks, validating the coordinate reference system at the I/O boundary, and applying memory backpressure so multi-day batch runs finish without crashing.

Audience prerequisites. You should be comfortable with Python 3.10+, virtual environments, and basic point-cloud concepts (XYZ, RGB, intensity, classification). The workflow here assumes a workstation with at least 16 GB RAM and an NVMe or SSD scratch disk — out-of-core streaming trades RAM for disk bandwidth, so spinning media will bottleneck throughput. This workflow is the memory-orchestration stage referenced by Automated Image Alignment & Feature Matching Workflows; it consumes the dense cloud that bundle adjustment and densification produce upstream.

# Prerequisites

Install the following into a clean virtual environment. Versions are the minimums validated for the streaming and VLR-parsing behavior used below.

laspy 2.x ships a streaming reader (chunk_iterator) and an appending writer that are the foundation of every snippet on this page. The lazrs extra is required to read and write compressed .laz; without it laspy can only handle uncompressed .las.

# Conceptual architecture

The core idea is that a point cloud never has to exist in memory in its entirety. Instead of laspy.read() (which materializes every point at once), you open the file as a reader, pull a fixed number of points per iteration, transform and filter that slice, append the survivors to an output file, and discard the slice before the next one arrives. Peak memory is then a function of chunk_size, not of total file size — a 12 GB cloud and a 1.2 TB cloud both process inside the same bounded working set.

Three guarantees must hold for this to be safe in production: the CRS is read once from the header and never silently reprojected mid-stream, every chunk is independently error-contained so one bad block cannot abort the survey, and a memory monitor applies backpressure before the OS starts swapping. The structuring conventions from structuring drone imagery for batch processing determine how flight blocks map to files here, and the datum handling mirrors coordinate transformation workflows in pyproj.

# 1. Estimate the memory budget before you load anything

Before choosing a chunk size, compute how much memory the full cloud would consume. A dense cloud from a 500-image survey easily exceeds 200 million points; at 48–64 bytes per point (XYZ, RGB, intensity, classification, optional normals) that is roughly 9.6–12.8 GB of contiguous array data before indexing or temporary buffers. This estimate tells you immediately whether streaming is mandatory.

def estimate_point_cloud_bytes(num_points: int, bytes_per_point: int = 56) -> float:
    """Return the approximate in-memory size (GiB) of a fully materialized cloud.

    bytes_per_point defaults to 56: XYZ (3x8) + RGB (3x2) + intensity (2) +
    classification (1) rounded up for struct padding. Transformations that copy
    arrays can momentarily double this, so size your headroom accordingly.
    """
    return (num_points * bytes_per_point) / (1024 ** 3)


if __name__ == "__main__":
    for n in (50_000_000, 200_000_000, 800_000_000):
        print(f"{n:>13,} points -> ~{estimate_point_cloud_bytes(n):6.1f} GiB resident")

If the resident estimate approaches half your physical RAM, switch to the streaming pipeline below rather than trying to load and decimate in one pass.

Figure 2 — Why chunk_size has a right answer rather than a large one. Past roughly two million points the time curve is flat, so every additional point in a chunk is memory spent for nothing.

# 2. Extract and validate the CRS from LAS/LAZ VLRs

Coordinate reference system mismatches are a silent cause of memory bloat: loading points without an explicit CRS lets downstream operations trigger implicit reprojections that duplicate coordinate arrays and inflate the footprint 2–3×. In laspy 2.x the CRS is not exposed as a single parse_crs() method on the header — it lives in Variable Length Records (VLRs). The portable source is the OGC WKT VLR (user_id="LASF_Projection", record_id=2112), which you parse with pyproj.CRS.from_wkt(). The broader rules for this live in managing coordinate reference systems in GDAL.

import logging
from typing import Optional

import laspy
from pyproj import CRS


def extract_crs_from_header(header: laspy.LasHeader) -> Optional[CRS]:
    """Extract a CRS from a LAS/LAZ header by scanning VLRs for the WKT record.

    Returns None when no recognizable projection VLR is present, leaving the
    caller to decide whether to proceed without reprojection or to abort.
    """
    for vlr in header.vlrs:
        if vlr.user_id == "LASF_Projection" and vlr.record_id == 2112:
            try:
                wkt = vlr.record_data.decode("utf-8").rstrip("\x00")
                return CRS.from_wkt(wkt)
            except Exception as exc:  # malformed or truncated WKT
                logging.warning("Failed to parse WKT VLR: %s", exc)
    return None

Validate the CRS once, at the I/O boundary. If the tie-point network was solved under a specific projection, forcing a late-stage datum shift introduces geometric distortion and needless memory duplication — transform at ingest or not at all.

Figure 3 — The CRS is not a field on the point cloud; it is a record beside it, in one of three places, chosen by a bit in the header. That indirection is the whole reason the validation step below exists.

# 3. Process points in bounded chunks with memory backpressure

Each chunk is filtered (and optionally reprojected) in isolation, and the function is wrapped so a single malformed block logs and raises without corrupting the output already on disk. The memory monitor reads system utilization and forces garbage collection before the OS reaches for swap.

import gc
import logging
from typing import Optional

import numpy as np
import psutil
from pyproj import Transformer


def memory_usage_percent() -> float:
    """Current system memory utilization as a percentage (0–100)."""
    return psutil.virtual_memory().percent


def apply_backpressure(threshold: float) -> None:
    """Collect garbage if memory is high; abort cleanly if it stays high."""
    if memory_usage_percent() > threshold:
        logging.warning("Memory at %.1f%%; collecting garbage.", memory_usage_percent())
        gc.collect()
        if memory_usage_percent() > threshold:
            raise MemoryError("Memory threshold breached; aborting before OOM kill.")


def process_chunk(points, transformer: Optional[Transformer] = None):
    """Filter (and optionally reproject) one chunk; return the surviving points."""
    if transformer is not None:
        # Transformer.transform returns a tuple of arrays, never a 2D array.
        points.x, points.y, points.z = transformer.transform(
            points.x, points.y, points.z
        )
    # Example filter: drop the lowest 2% of Z values (vegetation/noise floor).
    z_min = np.percentile(points.z, 2)
    return points[points.z >= z_min]

Replace the percentile filter with your own classification, statistical-outlier, or decimation step — the contract is simply “points in, fewer or transformed points out,” with no reference held to the previous chunk.

# 4. Stream the file out-of-core and serialize incrementally

The driver ties the pieces together: open the reader, read the CRS once, open an appending writer seeded with the source header, then iterate. Each processed chunk is flushed with append_points so peak memory stays bounded by chunk_size regardless of total file size. LAZ compression on write reclaims 60–80% of the disk footprint and lowers I/O pressure on the next read cycle.

import logging
from pathlib import Path

import laspy

logging.basicConfig(
    level=logging.INFO, format="%(asctime)s | %(levelname)s | %(message)s"
)


def stream_point_cloud(
    input_path: Path,
    output_path: Path,
    chunk_size: int = 750_000,
    memory_threshold: float = 82.0,
) -> int:
    """Stream a LAS/LAZ file in chunks, validate CRS, filter, and write out.

    Returns the number of points written. Peak memory is governed by chunk_size,
    so the same call handles a 12 GiB or a 1.2 TiB cloud unchanged.
    """
    if not input_path.exists():
        raise FileNotFoundError(f"Input file not found: {input_path}")

    written = 0
    with laspy.open(input_path) as reader:
        header = reader.header
        source_crs = extract_crs_from_header(header)
        if source_crs is None:
            logging.warning("No CRS VLR found; proceeding without reprojection.")
        else:
            logging.info("Validated CRS: %s", source_crs.name)

        with laspy.open(output_path, mode="w", header=header) as writer:
            for chunk in reader.chunk_iterator(chunk_size):
                apply_backpressure(memory_threshold)
                survivors = process_chunk(chunk)
                if len(survivors):
                    writer.append_points(survivors)
                    written += len(survivors)
                logging.info("Flushed chunk; running total %s points.", f"{written:,}")

    logging.info("Wrote %s points to %s", f"{written:,}", output_path)
    return written


if __name__ == "__main__":
    stream_point_cloud(
        Path("survey_block_01.laz"),
        Path("survey_block_01_processed.laz"),
        chunk_size=750_000,
        memory_threshold=82.0,
    )

For corridor-scale jobs that span many flight blocks, drive stream_point_cloud from the orchestration layer in parallel processing strategies for alignment, assigning one block per worker so RAM is bounded per process rather than across the whole survey.

# What “Bounded” Has To Mean In Practice

A streaming loop is only bounded if every allocation inside it is bounded, and the two that most often are not have nothing to do with the point data.

The first is accumulating results. It is natural to append per-chunk statistics, classified indices or filtered arrays to a list and reduce at the end — and that list grows with the file, which is precisely the property streaming was adopted to avoid. A cloud with 900 million points and a per-chunk array of surviving indices will happily reproduce the original memory profile with an extra step. Reductions must be incremental: running counts, running sums, running minima, and anything spatial written straight out to its tile rather than held.

The second is the writer’s own buffer. Appending to a compressed file is not free of state: the encoder holds a block until it is full, and some writers hold considerably more than one. On a long run with many small chunks this is invisible; on a run that writes wide records with many extra dimensions it is not. Measuring resident memory across a few thousand chunks rather than a few dozen is what exposes it, because a slow leak and a bounded loop look identical over a short window.

import gc

import psutil


class BoundedRun:
    """Assert that a streaming loop's memory is flat, not merely slow-growing.

    Samples RSS every `every` chunks and compares the late window against the
    early one. A loop that is genuinely bounded shows no trend; one that
    accumulates shows a monotone rise long before it reaches the host limit.
    """

    def __init__(self, every: int = 200, tolerance_mb: float = 64.0):
        self.every, self.tolerance_mb = every, tolerance_mb
        self.samples: list[float] = []
        self._proc = psutil.Process()
        self._n = 0

    def tick(self) -> None:
        self._n += 1
        if self._n % self.every:
            return
        gc.collect()                     # so the sample is not measuring garbage
        self.samples.append(self._proc.memory_info().rss / (1024 ** 2))

    def assert_flat(self) -> None:
        if len(self.samples) < 4:
            return                       # too short to conclude anything
        half = len(self.samples) // 2
        early = sum(self.samples[:half]) / half
        late = sum(self.samples[half:]) / (len(self.samples) - half)
        if late - early > self.tolerance_mb:
            raise MemoryError(
                f"streaming loop is accumulating: {early:.0f} MB → {late:.0f} MB")

Calling tick() inside the chunk loop and assert_flat() at the end turns “it worked on the test file” into a claim the pipeline re-checks on every run. The tolerance is deliberately generous — allocator behaviour is noisy and a 64 MB drift over a long run is not a leak — but a genuine accumulation crosses it early, while the job is still cheap to abandon.

The complementary discipline is to size the chunk from the record width, not from a point count. A LAS 1.2 point format 3 record is 34 bytes; a format 6 record with several extra dimensions can be well over a hundred. A chunk_size of two million points is a 68 MB working set in the first case and a 220 MB one in the second, from the same constant. Deriving the chunk size from header.point_format.size keeps the memory envelope stable across the heterogeneous files a real archive contains.

# Parameter deep-dive

# Verification and output inspection

Never trust a long batch run silently. After writing, assert that the point count is plausible, the output is readable, and — critically — the CRS survived the round trip unchanged. A reprojection that fired without your knowledge will show up as a CRS mismatch or an implausible bounding box.

import laspy

from pyproj import CRS


def verify_output(input_path, output_path) -> None:
    with laspy.open(input_path) as src, laspy.open(output_path) as dst:
        src_header, dst_header = src.header, dst.header

        # 1. The filter removes points, so output must be smaller but non-empty.
        assert 0 < dst_header.point_count <= src_header.point_count, (
            f"Implausible count: {dst_header.point_count} vs {src_header.point_count}"
        )

        # 2. The CRS must be byte-for-byte identical (no silent reprojection).
        src_crs = extract_crs_from_header(src_header)
        dst_crs = extract_crs_from_header(dst_header)
        assert src_crs == dst_crs, f"CRS drifted: {src_crs} -> {dst_crs}"

        # 3. The output bounding box must lie inside the input envelope.
        assert dst_header.mins[0] >= src_header.mins[0] - 1e-6
        assert dst_header.maxs[0] <= src_header.maxs[0] + 1e-6

    print("Verification passed: count, CRS, and bounds are consistent.")

Wire this assertion block into CI or a post-run hook so a regression in chunking or VLR handling fails loudly instead of producing a quietly corrupt deliverable.

# Troubleshooting

# Related

• Feature Detection Algorithms for Drone Imagery
• Optimizing Bundle Adjustment with Python
• Parallel Processing Strategies for Alignment
• Structuring Drone Imagery for Batch Processing
• Managing Coordinate Reference Systems in GDAL

← Automated Image Alignment & Feature Matching Workflows