How to Read a LAS or LAZ File in Python with laspy
On this page
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 10-million-point tile is 239 MB of coordinates alone as
float64, before any attribute. - Coordinates.
las.Xis the raw stored integer;las.xapplies the header's scale and offset. - Backends. Reading
.lazrequires 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}")9,977,603 points, LAS 1.2, point format 1
x 586500.00..587999.98 z 6.12..154.01
Eastings around 587,000 and northings around 4.5 million are UTM zone 18N — but nothing in las says so in a form laspy can read (step 2). Web-served copies of the same survey arrive in Web Mercator instead, with eastings around −8.2 million; if your coordinates look like that and you intend to measure anything, reproject first — see laspy returns the wrong coordinates.
A LAZ backend, the header first, scaled coordinates
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 implementationWithout 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()}")9,977,603 points
bounds [5.865e+05 4.515e+06 6.120e+00] .. [5.87999980e+05 4.51649998e+06 1.54010000e+02]
CRS None
This reads the header only — instant, against the fraction of a second to seconds it takes to decompress the points, which matters when surveying a directory of hundreds of tiles. And here it shows a trap: parse_crs() returned None, although the file does carry a CRS. It is written as user-defined GeoTIFF keys — NAD83(2011) / UTM zone 18N with NAVD88 heights, EPSG:6347 — which laspy cannot turn into a CRS object (laspy returns the wrong coordinates reads it).
3. Use x, not X
print(las.X.min(), las.X.max()) # raw int32
print(las.x.min(), las.x.max())
print(las.header.scales, las.header.offsets)58650000 58799998
586500.0 587999.98
[0.01 0.01 0.01] [0. 0. 0.]
The relationship is x = X * scale + offset, with scale and offset in the header. This file's offsets are zero, so X is exactly x times 100 — coordinates wrong by a factor of 100 that look like a projection problem. With a non-zero offset, as in the Web Mercator file in Example 1, X is not even proportional to x.
4. Stream large files rather than loading them
highest, n = -np.inf, 0
with laspy.open("tile.laz") as reader:
for points in reader.chunk_iterator(1_000_000):
z = np.asarray(points.z)
highest, n = max(highest, float(z.max())), n + len(z)
print(f"{n:,} points streamed, highest {highest:.2f} m")9,977,603 points streamed, highest 154.01 m
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.
A reader that reports and validates
The data used on this page
tile.laz is one 1,500 m tile of the USGS 3DEP survey of New York City after Hurricane Sandy (NY_New_York_CMGP_SANDY_LiDAR_15, tile 18TWL865150, around the Central Park reservoir), downloaded from the USGS rockyweb server. ept-root.laz is the root node of the same city's copy in the USGS Entwine Point Tile archive on AWS — a thin sample of all of New York, stored in Web Mercator.
import os
import subprocess
FILES = {
"tile.laz": "https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/LPC/Projects/"
"NY_New_York_CMGP_SANDY_LiDAR_15/NY_CMPG_2013/LAZ/USGS_LPC_NY_New_York_CMGP_SANDY_LiDAR_15_18TWL865150.laz",
"ept-root.laz": "https://s3-us-west-2.amazonaws.com/usgs-lidar-public/NY_NewYorkCity/ept-data/0-0-0-0.laz",
}
for name, url in FILES.items():
for _ in range(20): # rockyweb is slow and drops transfers: resume
if os.path.exists(name):
break
if subprocess.run(["curl", "-sfL", "-C", "-", "--retry", "3", "-o", name + ".part", url]).returncode == 0:
os.replace(name + ".part", name)
print({name: os.path.getsize(name) for name in FILES}){'tile.laz': 37341652, 'ept-root.laz': 390042}
Example 1 — a reader that reports and validates
import laspy
import numpy as np
# ASPRS codes common to LAS 1.2 and 1.4; 13 and above mean different things in
# different versions and surveys, so they are left to the survey's own report
CLASS_NAMES = {0: "never classified", 1: "unclassified", 2: "ground",
3: "low vegetation", 4: "medium vegetation",
5: "high vegetation", 6: "building", 7: "low noise",
8: "model key point", 9: "water"}
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 'none that laspy can read'}")
if crs is None:
print(" ! check the GeoTIFF keys and the survey metadata before measuring anything")
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, 'see the survey'):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
for path in ("tile.laz", "ept-root.laz"):
print(path)
data, header = read_cloud(path)tile.laz
LAS 1.2, point format 1, 9,977,603 points
scales [0.01 0.01 0.01], offsets [0. 0. 0.]
CRS none that laspy can read
! check the GeoTIFF keys and the survey metadata before measuring anything
1 unclassified 3,770,404 37.79%
17 see the survey 3,271,051 32.78%
2 ground 1,754,533 17.58%
18 see the survey 932,512 9.35%
9 water 240,767 2.41%
10 see the survey 6,407 0.06%
7 low noise 1,929 0.02%
389 MB in memory
ept-root.laz
LAS 1.2, point format 3, 33,699 points
scales [0.01 0.01 0.01], offsets [-8.23561e+06 4.96936e+06 3.37000e+02]
CRS EPSG:3857
! Web Mercator: distances are inflated by 1/cos(latitude)
2 ground 14,327 42.51%
18 see the survey 6,994 20.75%
9 water 6,136 18.21%
7 low noise 3,447 10.23%
1 unclassified 1,432 4.25%
17 see the survey 1,270 3.77%
10 see the survey 93 0.28%
1 MB in memory
The two warnings are the point of the function. The survey tile carries its CRS in a form parse_crs cannot read, so the code says so instead of guessing; the Entwine node is in Web Mercator, the usual CRS of point clouds served for the web and the one that silently breaks every metric result — at New York's latitude by a factor of 1.32. The class codes 17 and 18 are this survey's own: LAS 1.2 reserves them, and the per-tile metadata USGS publishes does not say what they mean, so the function does not guess.
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))
reservoir_box = (587_187, 4_515_077, 587_957, 4_515_851) # Jacqueline Kennedy Onassis Reservoir, EPSG:6347
x, y, z, c = read_bbox("*.laz", reservoir_box) 1 tiles read, 1 skipped by header bounds
1,629,344 points inside
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}")
ground = np.asarray(laspy.read("tile.laz").classification) == 2
write_subset("tile.laz", "ground.las", ground)
with laspy.open("ground.las") as check:
print(" ", check.header.point_count, check.header.scales, check.header.offsets,
[v.description for v in check.header.vlrs][:3]) 1,754,533 of 9,977,603 points (17.6%) -> ground.las
1754533 [0.01 0.01 0.01] [0. 0. 0.] ['GeoTiff Projection Keys', 'GeoTiff double parameters', 'GeoTiff ASCII parameters']
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.
X against x, backends and memory
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 9,977,603-point tile used here is 239 MB before a single attribute — classification, returns and intensity add more (Example 1 reports the total).
Two mitigations. Cast to float32 where 7 significant digits suffice, which halves the coordinate cost but is not safe for large projected coordinates: float32 steps are 0.0625 m at this tile's eastings, 0.5 m at its northings of 4.5 million, and 0.5 m at New York's Web Mercator eastings — against the 1 cm the file stores. 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 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 — 1.6 ms for nine in the clip-and-tile guide.
For a single large file, the equivalent does not exist, which is what COPC and EPT were designed to fix.
Missing CRSs and header checks
las.xis scaled;las.Xis raw. Off by 100 or 1,000.- Install a LAZ backend:
pip install "laspy[lazrs]". parse_crs()can returnNone. Older files often carry no CRS at all.- Check for Web Mercator. It is common on web-served clouds and breaks every distance.
chunk_iteratorfor large files. 10 M points is 239 MB of coordinates.float32is unsafe for large projected coordinates — 0.5 m steps at a northing of 4.5 million.- Preserve the source header when writing, or the scale changes silently.
- Filter on header bounds first when reading many tiles.
Internal links
- LAS, LAZ and COPC explained — the file layout behind these attributes
- A LAZ file will not open in Python — backend and version errors
- laspy returns the wrong coordinates: scale, offset and CRS — the scale and CRS traps
- LiDAR point clouds explained: returns, classes and intensity — what the attributes mean
- How to filter a point cloud by class and return number — the next step
- How to clip and tile a point cloud in Python — bounding-box workflows
- How to create a DTM from LiDAR ground points in Python — the first product
- How to open any spatial file in Python when you do not know the format — the general problem
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 10-million-point tile is 239 MB of coordinates alone.
How do I get the CRS?
las.header.parse_crs(). It can return None — some files carry no CRS, and some, like the New York tile here, carry one as user-defined GeoTIFF keys laspy cannot parse. Read the keys, or take the CRS 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 northing of 4.5 million, float32 steps are half a metre.