LAS, LAZ and COPC Explained: How Point Clouds Are Stored
On this page
Problem statement
A point cloud file is a header plus a long run of fixed-width records. That layout makes it fast to read from start to finish and slow to answer the question everyone actually asks: give me the points in this box.
The formats exist on a ladder that trades compression and indexing against simplicity:
- LAS — the uncompressed base format. Simple, large, no spatial index.
- LAZ — LAS compressed in chunks. Typically 5–10× smaller, still no spatial index.
- COPC — LAZ with an octree index inside it, so a bounding-box query at a chosen resolution is a range request.
- EPT — an octree stored as a directory of many small LAZ files plus a JSON hierarchy.
Measured on a real 3DEP survey: 12,968,770 points over 841 × 841 m arrived as 287 LAZ files totalling 82.8 MB — 6.4 bytes per point on disk, against 32 bytes per point for the same records uncompressed (Example 2).
Quick answer
import laspy
las = laspy.read("nodes/4-3-13-7.laz") # one node of the survey
h = las.header
print(f"LAS {h.version}, point format {h.point_format.id}")
print(f"{h.point_count:,} points")
print(f"scales {h.scales}")
print(f"offsets {h.offsets}")
print(f"raw X range {las.X.min()} .. {las.X.max()} (int32)")
print(f"scaled x {las.x.min():.2f} .. {las.x.max():.2f}")LAS 1.2, point format 1
180,756 points
scales [0.01 0.01 0.01]
offsets [-1.73345e+07 8.07826e+06 4.86000e+02]
raw X range -802212 .. -640954 (int32)
scaled x -17342522.12 .. -17340909.54
las.X and las.x are different attributes: the raw 32-bit integer and the scaled coordinate. Using the wrong one is off by a factor of 100 here, and by 1,000 in files with millimetre scaling.
Integer coordinates, point formats, compression and indexes
1. Understand why coordinates are integers
LAS stores X, Y and Z as signed 32-bit integers, with a scale and offset per axis in the header:
position = raw_integer * scale + offset
That is a deliberate design. Integers of a fixed scale give uniform precision across the file — every coordinate is exact to 1 cm here — whereas floating point gives relative precision, so a coordinate of 17 million metres would lose resolution exactly where a projected CRS puts you.
The cost is a fixed range. With a 0.01 scale, a signed 32-bit integer spans about ±21 million, which is why the offset exists: it moves the origin near the data.
2. Know what the point format determines
The point format ID fixes which attributes exist. Format 1 adds GPS time to the basic format 0; formats 2 and 3 add RGB; formats 6 and above (LAS 1.4) extend the return numbering and the classification range.
print([d.name for d in las.point_format.dimensions])
print(f"{las.point_format.size} bytes per record, "
f"{las.point_format.num_extra_bytes} of them extra bytes")['X', 'Y', 'Z', 'intensity', 'return_number', 'number_of_returns', 'scan_direction_flag', 'edge_of_flight_line', 'classification', 'synthetic', 'key_point', 'withheld', 'scan_angle_rank', 'user_data', 'point_source_id', 'gps_time', 'OriginId']
32 bytes per record, 4 of them extra bytes
OriginId is not part of point format 1: it is an extra-bytes dimension the EPT writer added to record which source file each point came from — four bytes on every record.
An older format cannot hold more than 5 returns per pulse or classification codes above 31. If you need those, the file must be LAS 1.4 with point format 6 or higher.
3. Know what LAZ does and does not add
LAZ compresses the point records in independent chunks, typically 50,000 points each. Compression is lossless and typically gets a factor of 5–10 — here 82.8 MB for records that would occupy 415 MB uncompressed, a factor of five.
What LAZ does not add is a spatial index. Finding the points in a box still means decompressing every chunk and testing. The chunking helps only in that a reader can skip whole chunks whose declared bounds miss the query, if the writer recorded them.
4. Use COPC or EPT when you need spatial queries
COPC stores the points in octree order inside a single LAZ file, with a hierarchy that maps each octree node to a byte range. A viewer can then fetch exactly the nodes overlapping the view, at exactly the depth the zoom level needs, using HTTP range requests.
EPT stores the same idea as a directory: ept.json, a hierarchy of JSON files, and one LAZ per node. The survey measured here is EPT, which is why it arrived as 287 separate files. Its hierarchy says how the whole survey is spread through the octree:
from collections import defaultdict
by_depth = defaultdict(lambda: [0, 0])
for key, count in hierarchy.items():
depth = int(key.split("-")[0])
by_depth[depth][0] += 1
by_depth[depth][1] += count
for depth, (n_nodes, n_points) in sorted(by_depth.items()):
print(f"depth {depth}: {n_nodes:5,} nodes, {n_points:12,} points")
print(f"total {sum(p for _, p in by_depth.values()):,} points")depth 0: 1 nodes, 19,873 points
depth 1: 5 nodes, 85,428 points
depth 2: 11 nodes, 400,774 points
depth 3: 31 nodes, 1,857,900 points
depth 4: 96 nodes, 9,189,283 points
depth 5: 317 nodes, 44,739,745 points
depth 6: 1,413 nodes, 156,264,268 points
depth 7: 4,720 nodes, 281,081,360 points
depth 8: 3,241 nodes, 35,646,686 points
total 529,285,317 points
The shallow depths are thin samples spread over the whole survey — so a low-resolution overview costs one small node rather than the whole file.
5. Check the CRS, and do not assume it is metric
The survey measured in laspy's wrong coordinates is stored in EPSG:3857, Web Mercator. Its extent reads as 1,613 × 1,613 m and the ground truth is 841 × 841 m — a factor of 1.92 at 58.6° north.
Every metric quantity computed in those units is wrong. See laspy returns the wrong coordinates.
A header report before reading points
The data used on this page
The USGS 3DEP lidar survey of Brooks Camp, Katmai National Park, Alaska (2012), from the USGS Entwine Point Tile archive on AWS, AK_BrooksCamp_2012. The setup reads the survey's EPT hierarchy, downloads every node under octree node 4-3-13-7 — 287 LAZ files in nodes/ — and joins them into tile.laz. The quick answer reads one of those files, nodes/4-3-13-7.laz, the node at the top of that branch.
import json
import os
import urllib.request
import laspy
EPT = "https://s3-us-west-2.amazonaws.com/usgs-lidar-public/AK_BrooksCamp_2012"
ROOT = (4, 3, 13, 7) # depth, x, y, z of the octree node
def under(key):
d, x, y, z = map(int, key.split("-"))
step = 2 ** (d - ROOT[0])
return d >= ROOT[0] and (x // step, y // step, z // step) == ROOT[1:]
hierarchy = json.load(urllib.request.urlopen(f"{EPT}/ept-hierarchy/0-0-0-0.json"))
nodes = sorted(k for k in hierarchy if under(k))
os.makedirs("nodes", exist_ok=True)
for key in nodes:
path = f"nodes/{key}.laz"
if not os.path.exists(path):
urllib.request.urlretrieve(f"{EPT}/ept-data/{key}.laz", path + ".part")
os.replace(path + ".part", path)
if not os.path.exists("tile.laz"): # the whole subtree as one file
with laspy.open(f"nodes/{nodes[0]}.laz") as first:
header = first.header
with laspy.open("tile.laz", mode="w", header=header) as writer:
for key in nodes:
writer.write_points(laspy.read(f"nodes/{key}.laz").points)
with laspy.open("tile.laz") as reader:
print(len(nodes), "nodes,", f"{reader.header.point_count:,} points,",
f"{sum(os.path.getsize(f'nodes/{k}.laz') for k in nodes) / 1e6:.1f} MB of nodes")287 nodes, 12,968,770 points, 82.8 MB of nodes
Example 1 — a header report before you read any points
import laspy
import numpy as np
def header_report(path):
"""Everything the header knows, without decompressing the points."""
with laspy.open(path) as reader:
h = reader.header
print(f" LAS {h.version}, point format {h.point_format.id}")
print(f" {h.point_count:,} points")
print(f" scales {h.scales}")
print(f" offsets {h.offsets}")
print(f" bounds x {h.mins[0]:.2f}..{h.maxs[0]:.2f}")
print(f" y {h.mins[1]:.2f}..{h.maxs[1]:.2f}")
print(f" z {h.mins[2]:.2f}..{h.maxs[2]:.2f}")
crs = h.parse_crs()
print(f" CRS {crs.to_string() if crs else 'NOT SET'}")
if crs and not crs.is_projected:
print(" ! geographic CRS — reproject before any metric computation")
extent = (h.maxs[0] - h.mins[0]) * (h.maxs[1] - h.mins[1])
print(f" nominal density {h.point_count / extent:.1f} per square unit"
f" (units of the file CRS)")
return h
header = header_report("nodes/4-3-13-7.laz") LAS 1.2, point format 1
180,756 points
scales [0.01 0.01 0.01]
offsets [-1.73345e+07 8.07826e+06 4.86000e+02]
bounds x -17342522.12..-17340909.54
y 8086328.13..8087940.58
z 15.94..200.65
CRS EPSG:3857
nominal density 0.1 per square unit (units of the file CRS)
Opening with laspy.open reads the header only. On an 80 MB file that is instant, against seconds to decompress every point — worth it when you are surveying a directory of tiles.
Example 2 — reading in chunks, so a large file fits
import laspy
import numpy as np
def chunked_stats(path, chunk_size=1_000_000, ground_class=2):
"""Stream a large LAZ without materialising it."""
total = 0
z_min, z_max = np.inf, -np.inf
class_counts = {}
with laspy.open(path) as reader:
for points in reader.chunk_iterator(chunk_size):
total += len(points)
z = np.asarray(points.z)
z_min = min(z_min, z.min())
z_max = max(z_max, z.max())
codes, counts = np.unique(points.classification, return_counts=True)
for code, count in zip(codes.tolist(), counts.tolist()):
class_counts[code] = class_counts.get(code, 0) + count
print(f" {total:,} points, z {z_min:.2f}..{z_max:.2f}")
for code, count in sorted(class_counts.items(), key=lambda kv: -kv[1]):
print(f" class {code:2d}: {count:12,} {count / total:6.2%}")
return class_counts
import os
class_counts = chunked_stats("tile.laz")
with laspy.open("tile.laz") as reader:
record = reader.header.point_format.size
n = reader.header.point_count
on_disk = sum(os.path.getsize(f"nodes/{k}.laz") for k in nodes)
print(f" {record} bytes per record uncompressed = {n * record / 1e6:.0f} MB; "
f"{on_disk / n:.1f} bytes per point in the 287 LAZ files") 12,968,770 points, z 15.88..203.11
class 1: 10,003,338 77.13%
class 2: 2,965,430 22.87%
class 7: 2 0.00%
32 bytes per record uncompressed = 415 MB; 6.4 bytes per point in the 287 LAZ files
A 13-million-point cloud is 311 MB as float64 coordinates alone. Streaming keeps memory flat, and for anything that reduces — histograms, grids, counts — costs nothing.
Example 3 — merging a directory of tiles safely
import glob
import laspy
import numpy as np
def merge_tiles(pattern, max_points=None):
"""Concatenate tiles, checking that they agree on scale and CRS."""
files = sorted(glob.glob(pattern))
if not files:
raise FileNotFoundError(pattern)
with laspy.open(files[0]) as reader:
reference = reader.header
scales, crs = reference.scales, reference.parse_crs()
xs, ys, zs, cls, ret, nret = [], [], [], [], [], []
total = 0
for path in files:
las = laspy.read(path)
if not np.allclose(las.header.scales, scales):
raise ValueError(f"{path}: scales {las.header.scales} "
f"differ from {scales}")
if las.header.parse_crs() != crs:
raise ValueError(f"{path}: CRS differs from the first tile")
xs.append(las.x); ys.append(las.y); zs.append(las.z)
cls.append(np.asarray(las.classification))
ret.append(np.asarray(las.return_number))
nret.append(np.asarray(las.number_of_returns))
total += len(las.points)
if max_points and total >= max_points:
break
print(f" merged {len(xs)} tiles, {total:,} points")
return (np.concatenate(xs), np.concatenate(ys), np.concatenate(zs),
np.concatenate(cls), np.concatenate(ret), np.concatenate(nret))
x, y, z, cls, ret, nret = merge_tiles("nodes/*.laz") merged 287 tiles, 12,968,770 points
Checking the scales matters more than it looks. Tiles written with different scales still merge without error in x (scaled) space, and produce a cloud that cannot be written back to a single LAS file without re-scaling.
Fixed-point integers, LAZ and octrees
Why fixed-point integers instead of floats
Floating point has relative precision: a float32 has about 7 significant digits, so at a UTM easting of 500,000 m its steps are about 3 cm, and at this survey's Web Mercator easting of 17,000,000 m they are 2 m.
Scaled integers have absolute precision. Every coordinate in the file measured here is exact to 1 cm regardless of position, because the offset moves the origin next to the data and the integer counts centimetres from there.
The trade is a bounded range and a mandatory header lookup. A reader that ignores the scale and offset gets numbers that are wrong by a factor of 100 and look plausible.
Why LAZ does not make queries fast
Compression reduces bytes; it does not reorganise. In a LAZ file the points are still in acquisition order, so the points inside any bounding box are scattered through every chunk.
The consequence is that "read a small area from a big file" costs almost as much as reading the whole file. That is exactly the operation a web viewer, a tiling job or an interactive tool needs, which is why the indexed formats exist.
Why an octree is the right index
Point clouds need two things from an index: spatial selection and level of detail. An octree gives both from one structure.
Each node holds a sample of the points in its cube, and its children hold finer samples of the eight sub-cubes. Reading to depth 3 gives a coarse view of everything; descending only where the viewport is gives detail where it is needed.
That is why the depth counts above rise to a peak and then fall: the deepest level holds only the leftover points in the densest regions, not a uniform layer.
Why merging tiles is not always safe
Adjacent tiles from one survey normally share a CRS, scale and offset. Tiles from different deliveries, or reprocessed at different times, may not.
The failure is silent. Scaled coordinates merge fine, and only the write-back fails — or worse, succeeds with a re-scaled precision that quietly loses a decimal place. Assert on scales and parse_crs() at merge time.
Raw versus scaled values and offsets
las.Xis the raw integer;las.xis scaled. Off by 100 or 1,000.- The offset is not the minimum. It is an arbitrary origin near the data.
- Point format limits the attributes. Formats below 6 cap returns at 5 and classes at 31.
- LAZ has no spatial index. Use COPC or EPT for bounding-box queries.
- CRS may be missing entirely.
parse_crs()returningNoneis common in older files. - Check the CRS is projected before computing anything metric.
- Stream large files with
chunk_iterator; 13 M points is 311 MB of coordinates alone asfloat64. - Assert scale and CRS agreement when merging tiles.
Internal links
- LiDAR point clouds explained: returns, classes and intensity — what the records contain
- How to read a LAS or LAZ file in Python with laspy — the practical reading guide
- laspy returns the wrong coordinates: scale, offset and CRS — the integer trap in detail
- A LAZ file will not open in Python — backend and version errors
- How to clip and tile a point cloud in Python — working around the missing index
- Point density explained — what the density figure needs
- Cloud-optimised GeoTIFF explained — the same indexing idea for rasters
- GIS vector file formats compared — the equivalent ladder for vector data
FAQ
What is the difference between LAS and LAZ?
LAZ is losslessly compressed LAS, typically 5–10× smaller. It adds no spatial index, so bounding-box queries still cost a full read.
What is COPC?
Cloud-Optimised Point Cloud: a LAZ file whose points are stored in octree order with a byte-range index, so a viewer can fetch just the region and resolution it needs over HTTP.
Why are point cloud coordinates stored as integers?
For uniform absolute precision. raw * scale + offset gives every point the same resolution regardless of position, which floating point cannot at large coordinate values.
Why is las.X different from las.x?
X is the raw 32-bit integer as stored; x applies the header's scale and offset. Confusing them is off by a factor of 100 or 1,000.
How do I read a huge LAZ file without running out of memory?
Use laspy.open(...).chunk_iterator(n) and reduce as you go. A 13-million-point cloud is 311 MB of coordinates alone as float64.
Can I merge LAS tiles by concatenating them?
Only after checking that they share a CRS, scale and offset. Mismatched scales merge silently in scaled space and fail or lose precision on write-back.
What is EPT?
Entwine Point Tile: an octree stored as a directory of small LAZ files plus a JSON hierarchy. The 3DEP survey used here arrived as 287 such files.