⚙️ How-To LiDAR & point clouds Files & formats LiDAR & Point Clouds · step 3 of 14 Updated

How to Read a LAS or LAZ File in Python with laspy

Problem statement

laspy.read() is one line and gets you an object with more than a dozen attributes, two of which are called X and x and differ by a factor of a hundred.

Three things go wrong reliably:

  • Memory. A 13-million-point cloud is 311 MB of coordinates alone as float64, before any attribute.
  • Coordinates. las.X is the raw stored integer; las.x applies the header's scale and offset.
  • Backends. Reading .laz requires a compression backend that is not installed by default, and the error message names neither.

Quick answer 

import laspy
import numpy as np

las = laspy.read("tile.laz")

x = np.asarray(las.x)                       # scaled coordinates, float64
y = np.asarray(las.y)
z = np.asarray(las.z)
classification = np.asarray(las.classification)
return_number = np.asarray(las.return_number)
number_of_returns = np.asarray(las.number_of_returns)

print(f"{len(las.points):,} points, LAS {las.header.version}, "
      f"point format {las.header.point_format.id}")
print(f"x {x.min():.2f}..{x.max():.2f}   z {z.min():.2f}..{z.max():.2f}")
180,756 points, LAS 1.2, point format 1
x -17342522.12..-17340909.54   z 22.10..197.34

An x coordinate of −17 million is a Web Mercator easting. If your coordinates look like that and you intend to measure anything, reproject first — see laspy returns the wrong coordinates.

A laspy object exposing raw integer X, Y, Z alongside scaled x, y, z, plus classification, return number and intensity arrays.
Uppercase is stored, lowercase is scaled. Everything else is a plain array of the same length.

Step-by-step solution

1. Install a LAZ backend

laspy alone reads uncompressed .las. For .laz it needs a backend:

pip install "laspy[lazrs]"      # pure Rust, no system dependency
pip install "laspy[laszip]"     # bindings to the reference implementation

Without one, opening a .laz raises LaspyException: No LazBackend selected, cannot decompress data. See A LAZ file will not open in Python.

2. Read the header first when you only need metadata 

with laspy.open("tile.laz") as reader:
    header = reader.header
    print(f"{header.point_count:,} points")
    print(f"bounds {header.mins} .. {header.maxs}")
    print(f"CRS {header.parse_crs()}")

This reads the header only. On an 80 MB file it is instant, against seconds to decompress the points — which matters when surveying a directory of hundreds of tiles.

3. Use x, not X 

print(las.X.min(), las.X.max())    # -802212 -640954     raw int32
print(las.x.min(), las.x.max())    # -17342522.12 -17340909.54

The relationship is x = X * scale + offset, with scale and offset in the header — 0.01 and −17,334,500 here. Using X by mistake gives coordinates that are wrong by a factor of 100 and look like a projection problem.

4. Stream large files rather than loading them 

with laspy.open("big.laz") as reader:
    for points in reader.chunk_iterator(1_000_000):
        process(np.asarray(points.z))

chunk_iterator yields point records without materialising the whole file. For anything that reduces — a histogram, a grid, a count — memory stays flat.

5. Convert to what you actually need, once

laspy attributes are lazily evaluated views. Wrapping them in np.asarray once, at the top, avoids repeated conversion inside loops and makes the memory cost explicit.

Reading a whole file into memory against streaming chunks and reducing, with memory flat in the streaming case.
Anything that reduces to a grid or a statistic should stream. Only interactive work needs the whole cloud.

Code examples

Example 1 — a reader that reports and validates 

import laspy
import numpy as np

CLASS_NAMES = {0: "never classified", 1: "unclassified", 2: "ground",
               3: "low vegetation", 4: "medium vegetation",
               5: "high vegetation", 6: "building", 7: "low noise",
               9: "water", 11: "snow", 18: "high noise"}


def read_cloud(path, require_projected=True):
    """Read a LAS/LAZ file into plain arrays, with the traps checked."""
    las = laspy.read(path)
    header = las.header

    crs = header.parse_crs()
    print(f"  LAS {header.version}, point format {header.point_format.id}, "
          f"{header.point_count:,} points")
    print(f"  scales {header.scales}, offsets {header.offsets}")
    print(f"  CRS {crs.to_string() if crs else 'NOT SET'}")

    if crs is None:
        print("  ! no CRS in the header — you must supply it externally")
    elif require_projected and not crs.is_projected:
        raise ValueError(f"{path} is in a geographic CRS ({crs.to_string()}); "
                         "reproject before any metric computation")
    elif crs and crs.to_epsg() == 3857:
        print("  ! Web Mercator: distances are inflated by 1/cos(latitude)")

    data = {
        "x": np.asarray(las.x), "y": np.asarray(las.y), "z": np.asarray(las.z),
        "classification": np.asarray(las.classification),
        "return_number": np.asarray(las.return_number),
        "number_of_returns": np.asarray(las.number_of_returns),
        "intensity": np.asarray(las.intensity),
    }
    if "gps_time" in las.point_format.dimension_names:
        data["gps_time"] = np.asarray(las.gps_time)
    if "point_source_id" in las.point_format.dimension_names:
        data["point_source_id"] = np.asarray(las.point_source_id)

    codes, counts = np.unique(data["classification"], return_counts=True)
    for code, count in sorted(zip(codes.tolist(), counts.tolist()),
                              key=lambda t: -t[1]):
        print(f"    {code:2d} {CLASS_NAMES.get(code, 'reserved'):18} "
              f"{count:11,}  {count / len(data['x']):6.2%}")

    print(f"  {sum(v.nbytes for v in data.values()) / 1e6:.0f} MB in memory")
    return data, header
  LAS 1.2, point format 1, 180,756 points
  scales [0.01 0.01 0.01], offsets [-1.73345e+07  8.07826e+06  4.86e+02]
  CRS EPSG:3857
  ! Web Mercator: distances are inflated by 1/cos(latitude)
     1 unclassified        139,352  77.10%
     2 ground               41,404  22.90%
  10 MB in memory

The Web Mercator warning is worth the three lines. It is the most common CRS on point clouds served for the web and the one that silently breaks every metric result.

Example 2 — reading a directory of tiles with a bounding-box filter 

import glob
import laspy
import numpy as np


def read_bbox(pattern, bounds, max_points=None):
    """Read only the tiles that intersect a box, and only points inside it."""
    left, bottom, right, top = bounds
    kept_x, kept_y, kept_z, kept_c = [], [], [], []
    tiles_read = tiles_skipped = 0

    for path in sorted(glob.glob(pattern)):
        with laspy.open(path) as reader:
            h = reader.header
            if (h.maxs[0] < left or h.mins[0] > right or
                    h.maxs[1] < bottom or h.mins[1] > top):
                tiles_skipped += 1
                continue

        tiles_read += 1
        las = laspy.read(path)
        x, y = np.asarray(las.x), np.asarray(las.y)
        inside = (x >= left) & (x <= right) & (y >= bottom) & (y <= top)
        if inside.any():
            kept_x.append(x[inside]); kept_y.append(y[inside])
            kept_z.append(np.asarray(las.z)[inside])
            kept_c.append(np.asarray(las.classification)[inside])
        if max_points and sum(len(a) for a in kept_x) >= max_points:
            break

    print(f"  {tiles_read} tiles read, {tiles_skipped} skipped by header bounds")
    if not kept_x:
        raise ValueError("no points in the requested box")
    print(f"  {sum(len(a) for a in kept_x):,} points inside")
    return (np.concatenate(kept_x), np.concatenate(kept_y),
            np.concatenate(kept_z), np.concatenate(kept_c))

Testing the header bounds before decompressing is the whole optimisation. LAZ has no spatial index, so this header-level filter is the cheapest approximation to one — and for a tiled delivery it is usually enough.

Example 3 — writing a filtered cloud back out 

import laspy
import numpy as np


def write_subset(source_path, out_path, mask, extra_dims=None):
    """Write a subset, preserving the header, scales and CRS."""
    las = laspy.read(source_path)
    mask = np.asarray(mask, dtype=bool)
    if len(mask) != len(las.points):
        raise ValueError(f"mask has {len(mask):,} entries, "
                         f"file has {len(las.points):,} points")

    out = laspy.LasData(las.header)
    out.points = las.points[mask]

    for name, values in (extra_dims or {}).items():
        out.add_extra_dim(laspy.ExtraBytesParams(name=name, type=values.dtype))
        setattr(out, name, values[mask])

    out.write(out_path)
    print(f"  {mask.sum():,} of {len(mask):,} points "
          f"({mask.mean():.1%}) -> {out_path}")

Constructing LasData from the source header keeps the scales, offsets and CRS. Building a fresh header instead means re-deriving all of them, and a wrong scale silently changes every coordinate.

Explanation

Why X and x both exist

LAS stores coordinates as signed 32-bit integers with a per-axis scale and offset in the header. That gives uniform absolute precision — every coordinate here is exact to 1 cm — where floating point would give relative precision and lose resolution at large coordinate values.

laspy exposes both: X is the stored integer, x is X * scale + offset as a float. The lowercase form is what you almost always want.

The one place X matters is writing: assigning to las.x re-derives the integers using the existing scale, and a value outside the representable range silently wraps.

Why LAZ needs a backend

LAZ compression is a specific chunked encoding, implemented in the reference laszip C++ library and reimplemented in Rust as lazrs. Neither ships with laspy, which is a pure-Python package.

laspy[lazrs] is usually the better choice: no system libraries and no compiler. laszip is the reference implementation, marginally faster on some files, and needs the shared library present.

Why memory disappears so fast

Coordinates alone, as float64, are 24 bytes per point. The 12,968,770-point cloud used here is 311 MB before a single attribute — classification, returns and intensity add more.

Two mitigations. Cast to float32 where 7 significant digits suffice, which halves the coordinate cost but is not safe for large projected coordinates: at an easting of 17 million, float32 resolves about a metre. And stream rather than load for anything that reduces.

Why header-level filtering matters so much

A LAZ file has no spatial index, so selecting points in a box means decompressing everything and testing.

For a tiled delivery — 287 files here — the header of each tile records its bounds, so most tiles can be rejected without decompression. That turns an all-tiles read into a few-tiles read for the cost of opening each header, which is microseconds.

For a single large file, the equivalent does not exist, which is what COPC and EPT were designed to fix.

Opening a LAS file for its header only against reading every point, with the header answering most questions.
For a tiled delivery, the header is the poor relation of a spatial index — and usually enough.

Edge cases or notes

  • las.x is scaled; las.X is raw. Off by 100 or 1,000.
  • Install a LAZ backend: pip install "laspy[lazrs]".
  • parse_crs() can return None. Older files often carry no CRS at all.
  • Check for Web Mercator. It is common on web-served clouds and breaks every distance.
  • chunk_iterator for large files. 13 M points is 311 MB of coordinates.
  • float32 is unsafe for large projected coordinates — about 1 m resolution at 17 million.
  • Preserve the source header when writing, or the scale changes silently.
  • Filter on header bounds first when reading many tiles.

FAQ

How do I read a LAZ file in Python?

laspy.read("file.laz"), with a compression backend installed — pip install "laspy[lazrs]". Then use las.x, las.y, las.z for scaled coordinates.

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

X is the raw stored 32-bit integer; x applies the header's scale and offset. They differ by a factor of 100 or 1,000.

How do I read a huge point cloud without running out of memory?

Use laspy.open(path).chunk_iterator(n) and reduce as you go. A 13-million-point cloud is 311 MB of coordinates alone.

How do I get the CRS?

las.header.parse_crs(). It can return None — many older files carry no CRS, and you must supply it from the survey metadata.

Can I read only part of a file?

Not by geometry, in plain LAS or LAZ — there is no spatial index. For a tiled delivery, filter on each tile's header bounds first. For a single large file, use COPC or EPT.

How do I write a filtered point cloud?

Build laspy.LasData from the source header and assign the masked points, so the scales, offsets and CRS are preserved.

Should I cast coordinates to float32?

Only for small local coordinates. At a Web Mercator easting of 17 million, float32 resolves about a metre.