CAD/GIS/BIM Interop

Reading LAS and LAZ Files with laspy

To read a LAS or LAZ file in Python, open it with laspy, read the header before allocating anything, and iterate the points in chunks using the lower-case x, y and z accessors so the header scale and offset are applied. The two failures that dominate first attempts are reading the raw integer arrays, which yields coordinates at the wrong origin and scale, and reading the whole file, which for an airborne delivery is an out-of-memory kill. This page is part of Point Cloud and Reality Capture Integration.

How LAS Stores a Point

A LAS file is a header followed by fixed-size point records. Each record stores X, Y and Z as 32-bit signed integers, and the header carries a scale and an offset per axis. The real coordinate is raw * scale + offset.

The two accessors and what each returns Two ways of reading the same ordinate. The upper-case accessor returns the stored 32-bit integer, which has the right shape and neither the right origin nor the right scale. The lower-case accessor applies the header scale and offset and returns the real coordinate. Nothing distinguishes them at the call site except one character. las.X — raw — the stored 32-bit integer — no scale, no offset — site appears kilometres across — looks like a unit error las.x — scaled — X * scale + offset — the real coordinate — sits where the survey says — what every consumer expects One character apart, a thousandfold and a relocation apart in the result.

This is why a file spanning kilometres can carry millimetre resolution in 32 bits, and it is why the distinction between accessors matters. las.X gives the stored integer; las.x applies the transformation. The difference does not raise, it just relocates the survey.

The header also declares the point format, which decides what attributes exist. Intensity and return number are present in every format; GPS time, colour, and near-infrared are not. Code that reads las.red on a format without colour fails, and code that assumes las.classification is meaningful succeeds on a file where nothing assigned it.

Finally, the coordinate reference system lives in the variable-length records, either as WKT or as legacy GeoTIFF keys depending on the LAS version. laspy exposes both through one accessor, which returns None when the file declares nothing.

Production-Ready Script

# laspy[lazrs]>=2.5, numpy>=1.24, Python 3.9+
from __future__ import annotations

from dataclasses import dataclass
import numpy as np
import laspy


@dataclass(frozen=True)
class CloudHeader:
    count: int
    point_format: int
    scales: tuple[float, float, float]
    offsets: tuple[float, float, float]
    mins: tuple[float, float, float]
    maxs: tuple[float, float, float]
    crs_wkt: str | None
    has_classification: bool
    has_colour: bool


def inspect(path: str) -> CloudHeader:
    """Header only — one seek, and it decides every later choice."""
    with laspy.open(path) as reader:
        h = reader.header
        dims = {d.name for d in h.point_format.dimensions}
        crs = h.parse_crs()
        return CloudHeader(
            count=h.point_count,
            point_format=h.point_format.id,
            scales=tuple(h.scales), offsets=tuple(h.offsets),
            mins=tuple(h.mins), maxs=tuple(h.maxs),
            crs_wkt=crs.to_wkt() if crs is not None else None,
            has_classification="classification" in dims,
            has_colour={"red", "green", "blue"} <= dims,
        )


def read_filtered(
    path: str,
    *,
    classes: set[int] | None = None,
    bbox: tuple[float, float, float, float] | None = None,
    chunk: int = 2_000_000,
) -> np.ndarray:
    """Chunked read returning an (n, 3) float64 array of REAL coordinates."""
    meta = inspect(path)
    if meta.crs_wkt is None:
        raise ValueError(f"{path}: no CRS in the variable-length records")
    if classes and not meta.has_classification:
        raise ValueError(f"{path}: point format {meta.point_format} has no classification")

    kept: list[np.ndarray] = []
    with laspy.open(path) as reader:
        for points in reader.chunk_iterator(chunk):
            mask = np.ones(len(points), dtype=bool)
            if classes:
                mask &= np.isin(points.classification, list(classes))
            # .x/.y/.z apply the header scale and offset; .X/.Y/.Z do not.
            x, y, z = np.asarray(points.x), np.asarray(points.y), np.asarray(points.z)
            if bbox:
                minx, miny, maxx, maxy = bbox
                mask &= (x >= minx) & (x <= maxx) & (y >= miny) & (y <= maxy)
            if mask.any():
                kept.append(np.column_stack((x[mask], y[mask], z[mask])))

    if not kept:
        return np.empty((0, 3), dtype=float)
    return np.vstack(kept)


if __name__ == "__main__":
    meta = inspect("survey.laz")
    print(f"{meta.count:,} points, format {meta.point_format}, "
          f"classification={meta.has_classification}")
    ground = read_filtered("survey.laz", classes={2})
    print(f"ground points: {len(ground):,}")
Reading a file larger than memory Four stages. The header is read alone, costing one seek and settling whether the file can be handled at all. The coordinate reference system is resolved or the read is refused. Points are then iterated in fixed-size chunks, and filtering happens inside the chunk so only retained points accumulate. Peak memory is the chunk plus the result, never the file. Read the header count, format, scale, bounds 1 one seek, no points Resolve the CRS or refuse the file 2 never assume the project system Iterate in chunks fixed size 3 peak memory is the chunk Filter inside the loop class and bounding box 4 only survivors accumulate

Key implementation notes:

  • inspect never reads a point. On a 400-million-point file that is the difference between a decision and a memory kill.
  • The CRS check raises rather than defaulting. An unlabelled cloud assumed into the project system is the point-cloud equivalent of a DXF assumed into millimetres.
  • Filtering happens inside the chunk loop, so peak memory is one chunk plus the retained points rather than the whole file.
  • np.isin handles a set of classes in one vectorised pass; a Python-level membership test per point dominates the runtime.
  • The classification availability check fails early with a useful message instead of raising an attribute error deep in the loop.

Compatibility Matrix

Component Supported range Notes
laspy >=2.5 chunk_iterator, parse_crs
LAZ support lazrs or laszip backend LAS works without one; LAZ does not
LAS versions 1.2 – 1.4 1.4 adds WKT CRS and extended point formats
Point formats 0 – 10 attribute availability varies; check before reading
numpy >=1.24 column stacking and boolean masking

Fallback Strategies

1. LaspyException on a LAZ file. No compression backend. Install laspy[lazrs], and assert the backend at start-up so the failure surfaces at deploy rather than on the first compressed delivery.

What each point format actually carries Four LAS point data record formats and the attributes each provides. Position, intensity and return number are universal; classification, GPS time and colour are not. Code that reads an attribute a format does not carry raises, and code that trusts a classification nothing assigned succeeds and returns nonsense. Format Classification GPS time Colour 0 yes — — 2 yes — yes 3 yes yes yes 7 extended yes yes Check availability from the header before reading; the failure is an attribute error mid-loop.

2. Coordinates are absurdly large. Raw accessors. Change las.X to las.x; the magnitude of the error is the reciprocal of the header scale, typically a thousand.

3. No CRS in the file. Common in processing intermediates and in PLY or PCD conversions. Obtain it from the survey report and record it explicitly rather than assuming — and prefer keeping the cloud in LAS so the metadata has somewhere to live.

4. Classification is all zeros. The producer never ran a classifier. Filtering on class 2 then returns nothing, which looks like a bug in the filter. Check the distinct classification values before relying on them.

5. The bounding box filter returns nothing. The bbox is in a different coordinate system from the cloud. Compare the filter bounds against meta.mins and meta.maxs before assuming the file is empty of the area of interest.

FAQ

What is the difference between las.x and las.X?

The lower-case accessors return real coordinates: the stored integer multiplied by the header scale and added to the header offset. The upper-case accessors return the raw stored integers. Reading the raw values gives a cloud with the correct shape at the wrong origin and the wrong scale — the classic symptom is a site that appears to be thousands of kilometres across.

Why does laspy refuse to open my LAZ file?

LAZ is compressed and laspy needs a compression backend to decode it. Install laspy[lazrs] or laspy[laszip]. Without one, LAS opens and LAZ raises — which in a container is a dependency that was present in development and absent in the image.

How do I read only the ground points?

Filter on the classification field, which uses the ASPRS class numbering — 2 is ground. Do it inside the chunk loop so only ground points accumulate. Treat the classification as the producer’s claim rather than as ground truth, and validate the resulting surface before building terrain from it.