Reading LAS into NumPy with laspy

TL;DR: las = laspy.read("tile.laz") loads every point; las.x, las.y, las.z are scaled float64 arrays, las.classification, las.intensity, las.return_number are integer arrays, and np.vstack((las.x, las.y, las.z)).T gives an (N, 3) array ready for SciPy or scikit-learn. Subset with boolean masks — las.points[las.classification == 2] — and use laspy.open first when you only need the header.

# Context and Motivation

This guide is part of laspy and NumPy Workflows for LAS Data. Most Python analysis of LiDAR starts by getting points into NumPy: a k-d tree for nearest neighbours, a histogram of heights, a scikit-learn model, a quick scatter plot in a notebook. laspy does that in one call and exposes each LAS field as a NumPy array with the right type. The two things people trip over are the difference between scaled and raw coordinates, and memory — a full read holds every point, and derived arrays multiply that.

# Prerequisites and Assumptions

• laspy 2.x with a LAZ backend: pip install "laspy[lazrs]".
• NumPy; pandas optional.
• A file that fits in memory. For larger files, see chunked reading of large LAS files with laspy.

# Step-by-Step Implementation

# Step 1 — Inspect the header first

laspy.open(path) reads only the header; check point_count and point_format before a full read.

# Step 2 — Read the file

laspy.read(path) returns a LasData with header, VLRs and points.

# Step 3 — Access fields as arrays

Standard fields use snake_case names (classification, return_number, number_of_returns, gps_time, point_source_id); extra dimensions use the names in their VLR.

# Step 4 — Subset with masks

A boolean array indexes las.points to produce a new record with only the selected points, keeping every field aligned.

# Step 5 — Build analysis arrays

Stack coordinates into an (N, 3) float array for SciPy and scikit-learn, or build a DataFrame for pandas.

# Complete Working Example

"""Read a LAZ tile into NumPy, subset it, and build analysis-ready arrays."""
from __future__ import annotations

import laspy
import numpy as np
import pandas as pd
from scipy.spatial import cKDTree

PATH = "tiles/t_0431.laz"

with laspy.open(PATH) as f:
    print(f"LAS {f.header.version}, PDRF {f.header.point_format.id}, "
          f"{f.header.point_count:,} points")
    print("dimensions:", list(f.header.point_format.dimension_names))

las = laspy.read(PATH)

# Scaled coordinates as float64 arrays.
xyz = np.vstack((las.x, las.y, las.z)).T
print("xyz shape:", xyz.shape, "dtype:", xyz.dtype)

# Ground points only, every field kept aligned.
ground = las.points[las.classification == 2]
print(f"ground: {len(ground):,} points, mean Z {np.mean(ground.z):.2f} m")

# Last returns above 2 m relative to the lowest ground, as a DataFrame.
last = (las.return_number == las.number_of_returns)
df = pd.DataFrame({
    "x": las.x[last], "y": las.y[last], "z": las.z[last],
    "intensity": las.intensity[last], "cls": las.classification[last],
})
print(df.describe().loc[["mean", "min", "max"]].round(2))

# Nearest-neighbour spacing on ground, from a k-d tree.
g_xy = np.vstack((ground.x, ground.y)).T
d, _ = cKDTree(g_xy).query(g_xy, k=2)
print(f"median ground spacing: {np.median(d[:, 1]):.2f} m")

# Extra dimensions, if present, by name.
for name in las.point_format.extra_dimension_names:
    print(name, las[name].dtype, float(np.nanmean(las[name])))

# Reading Only Part of a File

laspy.read has no spatial or attribute filter: it decompresses and loads every point, then you mask. That is fine for tiles of tens of millions of points on a workstation, and wasteful when you need a small area or a single class from a large file. There are three better routes, depending on the source.

For COPC files, use laspy.CopcReader and its query method with a bounding box and optional resolution; only the octree nodes that intersect the box are read, locally or over HTTP. For ordinary LAS or LAZ files, iterate with laspy.open(path).chunk_iterator(n) and keep only the points you need from each chunk, which bounds memory by the chunk size even though every point is still decompressed once. And when the subset is defined spatially by a polygon, or you need it repeatedly, let PDAL do the selection with filters.crop or filters.range and hand the result to Python as an array — PDAL streams, so the whole file never sits in memory.

Choosing well matters most in notebooks, where a casual laspy.read of a 3 GB tile is the most common way to exhaust a laptop’s memory and restart the kernel.

# Key Parameter Table

# Verification

• Count. len(las.points) equals las.header.point_count.
• Bounds. las.x.min() and las.x.max() match las.header.mins[0] and maxs[0] to within one scale unit; a mismatch means the header is stale.
• Types. las.classification.max() above 31 in a PDRF 0–5 file is impossible; if you see it, you are reading a different file than you think.

assert len(las.points) == las.header.point_count
assert abs(las.x.min() - las.header.mins[0]) <= las.header.scales[0]

# Gotchas and Edge Cases

Arithmetic on raw integers. las.Z - 100 subtracts 100 scale units — one metre at a 0.01 scale — not 100 metres. Always use lower-case fields for calculations.

Float32 precision. Casting UTM coordinates to float32 loses centimetres: at 4,471,000 m, float32 resolves only about 0.25 m. Keep coordinates float64, or subtract a local origin before casting.

Memory multiplies. np.vstack((las.x, las.y, las.z)) creates a new 24-byte-per-point array on top of the file’s records; building a DataFrame adds more. On large tiles, select fields and subsets before building derived arrays.

Legacy classification bits. In PDRF 0–5, classification is the 5-bit class; synthetic, key-point and withheld are separate flag fields (synthetic, key_point, withheld). In PDRF 6–10 the class is a full byte and overlap has its own flag.

# Frequently Asked Questions

How do I get an N by 3 array of coordinates from a LAS file?

Read it with laspy and stack the scaled coordinates: np.vstack of las.x, las.y and las.z, transposed. The result is a float64 array with one row per point.

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

las.X holds the raw 32-bit integers stored in the file; las.x applies the header’s scale and offset to give real-world coordinates as floats. Use the lower-case form for any calculation.

How do I read only ground points?

laspy reads the whole file, then you select with a mask such as las.points where classification equals 2. To avoid reading non-ground points at all, filter with PDAL or read in chunks and keep only matching points.

Can laspy read files larger than memory?

Not with laspy.read. Use laspy.open and chunk_iterator to process the file a fixed number of points at a time.

# Related

• laspy and NumPy Workflows for LAS Data — the wider picture
• Chunked Reading of Large LAS Files with laspy — when files do not fit
• Writing a LAS File from NumPy Arrays — the reverse direction
• Understanding LAS Point Data Record Formats — which fields each format holds
• How to Parse LAS Headers with Python — the header in detail