Troubleshooting Point Cloud Processing Failures

Point-cloud failures divide cleanly into two kinds, and the division decides how to look for them. Loud failures — a process killed, a pipeline raising, a file that will not open — announce themselves and are usually quick to fix once the message is read properly. Quiet failures produce a file that opens, renders and measures, and is wrong.

The quiet ones are the expensive ones, and almost all of them come from a small set of causes: a stage that silently emptied the cloud, a coordinate system that was never written, a classification that did not converge, or an attribute that was dropped in conversion. Each has a deterministic check that takes seconds.

This page aggregates the failure modes from across point cloud processing and 3D deliverables into one diagnostic sequence, with the check for each and a pointer to the page that covers it in depth.

Audience and prerequisites. Python 3.10+, PDAL and laspy available, and a failing or suspect run to look at. The checks below are ordered by cost, cheapest first.

# Prerequisites

# Conceptual architecture

A diagnostic sequence works best as a cascade: cheap checks that rule out whole families of cause, then expensive ones. For a suspect point-cloud run the order is almost always the same.

Start with the header, which is free and identifies the largest class of quiet failure — wrong scale, missing CRS, implausible bounds. Then the counts: total points, per-class counts, and how they compare with the input. Then the distribution: elevation percentiles, point density, extent. Only then the geometry, which means actually loading points and looking at neighbourhoods.

Most failures are identified in the first two steps, which together take under a second on a multi-gigabyte file.

Figure 1 — The cascade. Running it in this order is what keeps diagnosis a minute rather than an afternoon.

# Step 1: The header check

import laspy
import numpy as np


def header_diagnostics(path: str, *, expect_epsg: int | None = None,
                       max_scale: float = 0.005) -> dict:
    """Everything wrong that can be seen without decompressing a point."""
    with laspy.open(path) as fh:
        h = fh.header
    problems = []

    crs = h.parse_crs()
    if crs is None:
        problems.append("no CRS — the cloud will open at the origin")
    elif expect_epsg and crs.to_epsg() != expect_epsg:
        problems.append(f"CRS is EPSG:{crs.to_epsg()}, expected {expect_epsg}")

    if float(np.max(h.scales)) > max_scale:
        problems.append(f"scale {h.scales.tolist()} coarser than {max_scale} m")

    extent = np.asarray(h.maxs) - np.asarray(h.mins)
    if (extent <= 0).any():
        problems.append("degenerate bounds — the file may be empty or truncated")
    if extent[2] > 2000:
        problems.append(f"vertical extent {extent[2]:.0f} m — surviving outliers")
    if h.point_count == 0:
        problems.append("zero points")

    return {"points": int(h.point_count), "extent": extent.tolist(),
            "scales": h.scales.tolist(), "problems": problems}

The vertical-extent check is the highest-yield single line in this section. A photogrammetric cloud with a 2 km vertical range has surviving blunders, and those blunders will defeat every auto-scaled visualisation, every classification and every statistic downstream. The remedy is in filtering noise and outliers from dense clouds.

# Step 2: The count check

A stage that produced zero points is the commonest cause of an empty output, and pipelines do not report it — the file is written, it is simply empty, and every stage after it ran on nothing.

import json
import subprocess


def stage_counts(pipeline: dict) -> list[dict]:
    """Point count after each stage, by running truncated pipelines.

    Slower than a single run, but it turns "the output is empty" into "stage
    three emptied it", which is the whole diagnosis.
    """
    stages = pipeline["pipeline"]
    counts = []
    for i in range(1, len(stages) + 1):
        partial = {"pipeline": stages[:i] + [{"type": "filters.stats"}]}
        proc = subprocess.run(["pdal", "pipeline", "--stdin", "--metadata=/dev/stdout"],
                              input=json.dumps(partial), text=True,
                              capture_output=True)
        try:
            meta = json.loads(proc.stdout)
            n = meta["stages"]["filters.stats"]["statistic"][0]["count"]
        except Exception:
            n = None
        label = stages[i - 1].get("type", "reader") if isinstance(stages[i - 1], dict) else "reader"
        counts.append({"stage": i, "type": label, "points": n})
        if n == 0:
            break
    return counts

Common culprits: a filters.range whose limits exclude everything, a filters.crop with bounds in the wrong CRS, or a classification-based filter applied before the classification ran. All three produce an empty result and no error. Diagnosing empty output from a PDAL pipeline covers the full set.

# Step 3: The distribution check

import json

import numpy as np
import pdal


def distribution_diagnostics(path: str, sample: int = 2_000_000) -> dict:
    """Percentiles, density and class balance from a sample of the cloud."""
    pipe = pdal.Pipeline(json.dumps({"pipeline": [
        path, {"type": "filters.sample", "radius": 0.0}]}))
    pipe.execute()
    a = pipe.arrays[0]
    if a.size > sample:
        idx = np.random.default_rng(0).choice(a.size, size=sample, replace=False)
        a = a[idx]

    z = a["Z"]
    classes, counts = np.unique(a["Classification"], return_counts=True)
    area = (a["X"].max() - a["X"].min()) * (a["Y"].max() - a["Y"].min())

    problems = []
    tail = float(z.max() - np.percentile(z, 99.9))
    if tail > 50:
        problems.append(f"top 0.1 % of points span {tail:.0f} m — blunders survive")
    if len(classes) == 1:
        problems.append(f"every point is class {int(classes[0])} — "
                        "classification did not run or did not converge")
    ground = dict(zip(classes.tolist(), counts.tolist())).get(2, 0)
    if ground and ground / a.size < 0.05:
        problems.append(f"ground is {ground / a.size:.1%} of points")

    return {"z_percentiles": np.percentile(z, [1, 50, 99, 99.9]).tolist(),
            "density_per_m2": float(a.size / max(area, 1e-9)),
            "classes": dict(zip(classes.tolist(), counts.tolist())),
            "problems": problems}

# Step 4: Memory failures, which are their own category

A process killed without a Python traceback is almost always the kernel’s out-of-memory killer, and the diagnosis is structural rather than incidental: some stage in the chain cannot stream. Checking headroom before starting a blocking stage converts a crash into a useful error.

import psutil


def assert_memory_headroom(point_count: int, *, bytes_per_point: int = 64,
                           factor: float = 2.5) -> None:
    """Refuse to start a blocking stage that clearly will not fit.

    64 bytes per point is a realistic in-memory figure once coordinates,
    classification and a few computed dimensions are present; the factor
    covers the temporaries a neighbourhood filter allocates.
    """
    need = point_count * bytes_per_point * factor
    available = psutil.virtual_memory().available
    if need > available:
        raise MemoryError(
            f"this stage needs about {need / 1e9:.1f} GB and {available / 1e9:.1f} GB "
            "is available — tile the input or use a streaming alternative")

The remedies are in resolving memory errors in PDAL and laspy and the streaming analysis in choosing a library for streaming versus in-memory work.

# Step 5: Reading a PDAL error message properly

PDAL’s diagnostics are informative once their structure is understood, and opaque until then. Three shapes account for most of what a pipeline reports.

A schema or option error names a stage and a key. filters.smrf: Option 'windowsize' not recognized means exactly what it says: the option name is wrong, usually because a parameter was renamed between versions or copied from documentation for a different filter. These are cheap to fix and never indicate a data problem.

A GDAL or PROJ error surfaces through PDAL without much context. Cannot find proj.db or PROJ: internal_proj_create: no database context are environment failures, not pipeline failures — the code is correct and the container is missing its PROJ data. The same message appears in the raster stack, and the fix is documented in fixing GDAL PROJ database context errors.

A message about dimensions means a stage was asked for something not present. Dimension 'HeightAboveGround' does not exist almost always means a stage order problem: the filter that creates the dimension runs after the one that consumes it. PDAL executes stages in the order written, and the error is the first sign that the order is wrong.

import json
import subprocess


def run_pipeline(pipeline: dict) -> dict:
    """Run a pipeline and classify any failure into an actionable category."""
    proc = subprocess.run(["pdal", "pipeline", "--stdin"],
                          input=json.dumps(pipeline), text=True,
                          capture_output=True)
    if proc.returncode == 0:
        return {"ok": True}

    err = (proc.stderr or "").strip()
    if "not recognized" in err or "Schema" in err:
        category = "pipeline definition — an option name or type is wrong"
    elif "proj" in err.lower() or "GDAL" in err:
        category = "environment — PROJ or GDAL data is missing or mismatched"
    elif "does not exist" in err and "Dimension" in err:
        category = "stage order — a dimension is consumed before it is created"
    elif "std::bad_alloc" in err or "Killed" in err:
        category = "memory — a stage could not stream"
    else:
        category = "unclassified"
    return {"ok": False, "category": category, "stderr": err[:500]}

Classifying the failure before reacting to it is worth the twelve lines, because the four categories have entirely different remedies and the temptation is always to treat every failure as a data problem.

# Step 6: Building the checks into the run rather than the runbook

Every check on this page is a few lines and runs in seconds. The difference between a team that ships quiet failures and one that does not is almost never knowledge of the checks — it is whether they execute automatically.

Three integration points work well. After every write, run the header diagnostics; it is free and catches the largest class of delivery fault. Before every blocking stage, assert memory headroom; it converts an unexplained kill into a message naming the stage and the shortfall. Before every delivery, run the full cascade and compare against the previous run for the same site.

from pathlib import Path


def pipeline_with_gates(src: str, dst: str, stages: list[dict],
                        *, expect_epsg: int, expect_density: float,
                        previous: dict | None = None) -> dict:
    """A processing run with the diagnostics wired in rather than documented."""
    pipeline = {"pipeline": [src, *stages,
                             {"type": "writers.las", "filename": dst,
                              "compression": "laszip", "forward": "all"}]}

    result = run_pipeline(pipeline)
    if not result["ok"]:
        raise RuntimeError(f"{result['category']}: {result['stderr']}")

    report = gate_output(dst, expect_epsg=expect_epsg,
                         expect_density=expect_density,
                         previous_points=(previous or {}).get("points"))
    report["output_bytes"] = Path(dst).stat().st_size
    return report

The run record this produces is the same artefact the rest of this section keeps asking for. A month later, when somebody questions a deliverable, the answer to “was this checked” is a stored report rather than a recollection — and on the occasions when something did slip through, the stored reports are what make it possible to find out when it started.

A final note on discipline: resist the urge to lower a threshold when it fires. A ground-fraction gate that trips on one site is telling you something about that site, and the right response is to look rather than to relax the number. Thresholds that have been quietly widened over a year are indistinguishable from no thresholds at all.

# Parameter deep-dive

# Verification and output inspection

Bundling the cascade into one function that runs on every output is what turns these from debugging tools into a gate.

def gate_output(path: str, *, expect_epsg: int, expect_density: float,
                previous_points: int | None = None) -> dict:
    """Run the cascade and fail the job on anything in it."""
    report = {"header": header_diagnostics(path, expect_epsg=expect_epsg)}
    if report["header"]["problems"]:
        raise ValueError("; ".join(report["header"]["problems"]))

    report["distribution"] = distribution_diagnostics(path)
    if report["distribution"]["problems"]:
        raise ValueError("; ".join(report["distribution"]["problems"]))

    got = report["distribution"]["density_per_m2"]
    if got < 0.2 * expect_density:
        raise ValueError(f"density {got:.0f}/m² is far below the expected "
                         f"{expect_density:.0f}/m² — reconstruction degraded")
    if previous_points and report["header"]["points"] < 0.5 * previous_points:
        raise ValueError("point count halved against the previous run")
    return report

Comparing against the previous run is the addition worth making on a monitoring contract. Absolute thresholds catch gross failures; a comparison against last month catches the gradual degradation that absolute thresholds never trip.

# A worked diagnosis

A concrete example shows how quickly the cascade closes a question. A monthly run produces a DTM with no data over half the site, and the job reported success.

The header check passes: CRS present, scale 1 mm, bounds plausible, 214 million points. So the file is not empty and the coordinates are sound — which already rules out the format failures entirely.

The count check shows 214 million points in and 214 million out, so nothing was removed. But the class histogram from the distribution check shows 211 million points in class 1 and three million in class 2. Ground is 1.4 % of the cloud, below the five percent floor, and the gate fires.

That narrows it to the classifier in one step. Looking at the parameters recorded with the run, the SMRF window is 6 m where the previous month used 18 m — somebody changed it while investigating a building that appeared in the terrain, and the change shipped. The small window treats every gentle rise as an object, so most of the site classified as non-ground and the DTM has nothing to interpolate from.

Total elapsed time: under two minutes, none of it spent in a viewer. The three things that made it quick were that the parameters were recorded with the output, that the class histogram was computed automatically, and that the ground-fraction gate had a number attached rather than a judgement.

Figure 3 — The expensive failures are the ones that exit zero.

Figure 4 — Four inspections, and the last is the one nobody performs.

# Troubleshooting

The output file is empty and nothing raised. A stage emptied the cloud. Run the per-stage count and look for the first zero; a range filter with wrong limits or a crop in the wrong CRS are the usual causes.

The cloud opens at the origin. No CRS in the file, or one written in an encoding the reader does not support. See fixing SRS lost when writing LAS files.

The process was killed with no traceback. The kernel’s OOM killer. A stage cannot stream; tile the input or size the machine.

Every point is the same class. Either classification never ran, or a preceding filter removed everything it would have classified. Check the per-stage counts and the stage order.

The classification changes between identical runs. The classification was not reset before filtering, so each run built on the last. Begin the pipeline by assigning class zero.

Elevation statistics are dominated by a handful of points. Blunders survived filtering. Check the 99.9th percentile against the maximum; a large gap between them is the signature.

Attributes the pipeline computed are missing from the output. An intermediate stage dropped undeclared extra dimensions. Declare them explicitly on every writer.

# Related

• Diagnosing Empty Output from a PDAL Pipeline
• Fixing SRS Lost When Writing LAS Files
• Resolving Memory Errors in PDAL and laspy
• Classifying Point Clouds with PDAL and Python
• Point Cloud Formats and Interchange in Python

← Point Cloud Processing & 3D Deliverables