Choosing a Library for Streaming Versus In-Memory Work

The pipeline ran fine for a year and then a 900 million point survey arrived and the process was killed. Nothing changed in the code. What changed is that one stage in the chain cannot stream, and until that survey nothing had been large enough for it to matter.

Streaming is the property that decides whether a point-cloud job runs on a workstation or needs a machine nobody has. It is also poorly understood, because it is not a property of a library — it is a property of each individual operation, and one non-streaming stage anywhere in a chain forces the whole cloud into memory.

This page covers which operations can stream and which fundamentally cannot, how to check, and how to write the chunked alternative correctly when streaming is unavailable. It extends the library comparison in PDAL vs laspy for point cloud automation.

What makes an operation streamable

An operation can stream if its output for a given point depends only on that point and on a bounded amount of preceding context. Assigning a classification from a range test streams: each point’s answer depends on nothing else. Reprojecting streams. Writing streams.

An operation cannot stream if it needs to see points that may arrive later, or all of them. A statistical outlier filter needs the distribution of the whole cloud before it can judge any point. Ground classification needs a neighbourhood that extends in every direction. Sorting, obviously, needs everything.

The distinction has a practical consequence that is easy to miss: a chain streams only if every stage in it does. A pipeline of four streaming stages and one blocking stage has the memory profile of the blocking stage, and the streaming ones contribute nothing beyond a false sense of safety.

Memory profile of a chain with one non-streaming stage Two pipeline chains drawn with their memory profiles beneath. The upper chain has four streaming stages and holds one chunk at a time, about two hundred megabytes, regardless of the cloud's size. The lower chain is identical except that the third stage is a statistical outlier filter, which must see the whole distribution; its memory profile rises to the full cloud size of thirty-one gigabytes at that stage and stays there. A note states that the surrounding streaming stages contribute nothing once one stage blocks. all stages stream — 200 MB peak read reproject range write 200 MB one blocking stage — 31 GB peak read reproject outlier needs everything write 31 GB The three streaming stages around it contribute nothing. A chain is as blocking as its worst stage.

Figure 1 — Why “the pipeline streams” is not a property a pipeline has unless every stage does.

Minimal reproducible solution

PDAL will tell you whether a pipeline streams, which turns a guess into a check.

import json
import subprocess


def pipeline_streams(pipeline: dict) -> dict:
    """Ask PDAL whether this pipeline can run in streaming mode.

    `--stream` fails loudly on a pipeline that cannot stream rather than
    silently falling back, which is exactly the behaviour you want in CI: a
    change that breaks streaming should break the build, not production.
    """
    proc = subprocess.run(
        ["pdal", "pipeline", "--stdin", "--stream", "--dryrun"],
        input=json.dumps(pipeline), text=True, capture_output=True)
    ok = proc.returncode == 0
    return {"streams": ok,
            "message": (proc.stderr or "").strip()[:300] if not ok else ""}


NON_STREAMING = {
    "filters.outlier",          # needs the whole distribution
    "filters.smrf",             # neighbourhood in every direction
    "filters.pmf",
    "filters.hag_nn",           # nearest ground neighbour, unbounded search
    "filters.covariancefeatures",
    "filters.sort",
    "filters.voxelcenternearestneighbor",
}


def explain_blocking(pipeline: dict) -> list[str]:
    """Name the stages that force the whole cloud into memory."""
    return [s["type"] for s in pipeline["pipeline"]
            if isinstance(s, dict) and s.get("type") in NON_STREAMING]

Keeping an explicit list alongside the runtime check is worth it because the check requires PDAL to be installed and the list can be consulted at design time — when the pipeline is being written and the decision is cheap.

When streaming is unavailable: tiling, done correctly

For a blocking operation on a cloud that does not fit, the answer is tiling with a halo. The halo is the part that goes wrong, and it goes wrong in a way that produces a plausible result with a regular grid of artefacts.

import json
import subprocess
from pathlib import Path


def tiled_blocking_op(src: str, out_dir: str, stage: dict,
                      *, tile: float = 500.0, halo: float = 0.0) -> None:
    """Run a non-streaming stage tile by tile, with a halo, then trim it.

    The halo must exceed the stage's own spatial reach — the SMRF window, the
    outlier neighbourhood radius, the covariance knn distance. Set it below
    that and every tile edge is computed against a truncated neighbourhood,
    producing a grid of artefacts at exactly the tile spacing.
    """
    reach = {"filters.smrf": stage.get("window", 18.0),
             "filters.outlier": stage.get("radius", 1.0) * 3,
             "filters.covariancefeatures": 5.0}.get(stage["type"], 20.0)
    if halo <= reach:
        raise ValueError(
            f"halo {halo} m does not exceed the stage's reach of {reach} m — "
            "expect seams at the tile spacing")

    Path(out_dir).mkdir(parents=True, exist_ok=True)
    pipeline = {"pipeline": [
        src,
        {"type": "filters.splitter", "length": tile, "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)

Deriving the required halo from the stage’s own parameters, rather than accepting one from the caller, is what makes this safe. A halo chosen by hand is a halo that will be wrong the first time somebody changes the window size.

The laspy equivalent, and what it costs

Where the operation is not a PDAL stage, the same structure has to be written by hand. It is not difficult, but it is more code than it appears, and the halo bookkeeping is the whole of it.

import numpy as np
import laspy


def chunked_statistic(path: str, chunk_points: int = 5_000_000):
    """Stream a per-point statistic that needs no neighbourhood.

    This is the easy case: each chunk is independent, memory is bounded by
    chunk size, and there is no halo to manage. Anything needing neighbours
    requires spatial tiling instead of sequential chunking, because a
    sequential chunk is not spatially compact.
    """
    total = 0
    running_sum = 0.0
    with laspy.open(path) as fh:
        for points in fh.chunk_iterator(chunk_points):
            z = np.asarray(points.z)
            running_sum += float(z.sum())
            total += z.size
    return {"mean_z": running_sum / total, "points": total}

The comment marks the distinction that matters most. chunk_iterator yields points in file order, which for most LAS files is roughly acquisition order — not spatial order. A neighbourhood operation over a sequential chunk therefore sees a scattered subset of the site, not a compact region, and produces nonsense. Spatial tiling is a different and more expensive operation, and it is the reason PDAL’s splitter earns its keep.

Streaming and in-memory work compared on what each demands Two columns. The streaming column notes that memory stays bounded regardless of file size, that operations must be expressible per chunk, that a second pass costs a second read, and that it suits filtering, reprojection and format conversion. The in-memory column notes that memory scales with the point count, that any operation is available including ones needing global context, that repeated passes are free once loaded, and that it suits neighbourhood statistics, clustering and iterative fitting. streaming memory bounded whatever the file size operations must work per chunk a second pass costs a second read filtering, reprojection, conversion in memory memory scales with the point count any operation, including global ones repeated passes are free once loaded neighbourhoods, clustering, iterative fits The question is whether the operation needs global context, not which library is faster.

Figure 3 — The operation decides the model, and the model decides the library.

Edge-case matrix

Operation Streams? If not, what to do
Range / assign filters Yes —
Reprojection Yes —
Writing LAS/LAZ Yes —
Writing COPC No Needs the whole cloud to build the octree; size the machine
Statistical outlier No Tile with a halo, or use radius method per tile
SMRF / PMF ground No Tile with a halo exceeding the window
Height above ground No Tile with a generous halo
Covariance features No Tile with a halo above the knn reach
Sorting No Rarely necessary; avoid
Sequential chunking for neighbourhoods N/A Wrong tool — chunks are not spatially compact

Verification snippet

The failure that tiling introduces is visible as a periodic pattern, and it can be detected without looking at anything.

import numpy as np


def detect_tile_seams(values: np.ndarray, transform, tile_size: float,
                      *, z_threshold: float = 4.0) -> dict:
    """Look for structure at exactly the tile spacing, which is never real.

    Compares the mean absolute gradient along rows that fall on tile
    boundaries against the overall mean. Terrain has no reason to change at
    the tile spacing, so an excess there is a processing artefact.
    """
    gy, gx = np.gradient(np.nan_to_num(values))
    grad = np.hypot(gx, gy)
    step = max(int(round(tile_size / abs(transform.a))), 1)

    boundary = grad[:, ::step]
    overall_mean = float(np.nanmean(grad))
    overall_sd = float(np.nanstd(grad))
    boundary_mean = float(np.nanmean(boundary))

    z = (boundary_mean - overall_mean) / max(overall_sd, 1e-9)
    return {"boundary_mean": boundary_mean, "overall_mean": overall_mean,
            "z_score": z,
            "seams_detected": bool(z > z_threshold),
            "note": ("gradient is elevated at the tile spacing — increase the halo"
                     if z > z_threshold else "no periodic structure detected")}
Why a halo smaller than the operation's reach produces seams Two tile diagrams. In the first, the halo is smaller than the neighbourhood the operation needs, so a point near the tile edge is evaluated against a neighbourhood that is cut off on one side, and the result differs systematically from the same point evaluated in a neighbouring tile. In the second, the halo exceeds the operation's reach, so every point inside the tile proper sees a complete neighbourhood and the trimmed tiles join seamlessly. A note gives the required halo for the common filters. halo < reach — seams halo neighbourhood is cut off at the halo edge halo > reach — clean complete neighbourhood inside the halo required halo, by filter SMRF: above the window · outlier: three times the radius covariance: above the knn reach · height above ground: generous

Figure 2 — The halo rule, and the artefact it prevents.

When to escalate

  • A required stage cannot stream and the cloud will not tile cleanly. Some operations — COPC writing, global sorting — genuinely need everything. Size the machine for the largest survey rather than engineering around it.
  • Tiling produces seams even with a generous halo. The stage’s reach is larger than documented, or the operation is not local at all. Test on a single tile against the whole-cloud result before scaling.
  • Memory is fine and the job is still slow. Streaming is about feasibility, not speed. Profile the stages, per benchmarking PDAL pipelines against laspy loops.

← PDAL vs laspy for Point Cloud Automation