Converting Point Clouds to 3D Tiles with py3dtiles

This page converts classified LAZ point clouds into a streamable 3D Tiles point cloud tileset with py3dtiles — preparing the input with PDAL so heights are ellipsoidal before the converter sees them, running the conversion from UTM 18N (EPSG:32618) to Earth-centred EPSG:4978, sizing jobs and cache for the machine, and auditing the output .pnts files so every input point is accounted for.

Why you hit this

A browser cannot open a 40 GB LAZ collection, and a digital twin should not ask it to. Point cloud tilesets split the cloud into an octree of small binary tiles that a viewer streams by screen-space error, which is how Cesium, iTwin and deck.gl display city-scale scans. py3dtiles is the most direct open-source route from LAS to that format in a Python pipeline. Its defaults are good; what goes wrong is almost always the input — mixed CRSs across files, orthometric heights treated as ellipsoidal, noise points that stretch the octree’s bounding box a kilometre into the sky. The octree addressing underneath the output is covered in octree indexing point clouds with Morton codes.

Prerequisites

  • py3dtiles>=7.0 with LAZ support (pip install "py3dtiles[las]"), pdal>=3.4 Python bindings on PDAL 2.6+, pyproj>=3.6, numpy>=1.24.
  • Input LAZ in EPSG:32618 with NAVD88 heights (compound EPSG:32618+5703), classified to ASPRS codes.
  • The GEOID18 grid available to PROJ (us_noaa_g2018u0.tif), via PROJ_NETWORK=ON or a local copy.
  • Free disk of roughly the input LAZ size: py3dtiles writes uncompressed intermediate node files while it builds the tree.

Step-by-Step

1. Prepare the input with PDAL: drop noise, make heights ellipsoidal

python
import json
from pathlib import Path
import pdal

RAW = sorted(Path("raw").glob("*.laz"))
PREP = Path("prep"); PREP.mkdir(exist_ok=True)

for src in RAW:
    pipeline = {
        "pipeline": [
            {"type": "readers.las", "filename": str(src), "override_srs": "EPSG:32618+5703"},
            {"type": "filters.range", "limits": "Classification![7:7],Classification![18:18]"},
            {"type": "filters.reprojection",
             "in_srs": "EPSG:32618+5703",
             "out_srs": "+proj=utm +zone=18 +datum=WGS84 +units=m +type=crs"},
            {"type": "writers.las", "filename": str(PREP / src.name),
             "a_srs": "EPSG:32618", "compression": "laszip", "forward": "all",
             "scale_x": 0.001, "scale_y": 0.001, "scale_z": 0.001, "offset_x": "auto",
             "offset_y": "auto", "offset_z": "auto"},
        ]
    }
    count = pdal.Pipeline(json.dumps(pipeline)).execute()
    print(f"{src.name}: {count:,} points written")

Two things happen here that py3dtiles cannot do for you. Classes 7 and 18 — low and high noise — are removed, because a single bird return 900 m above the city expands the root bounding volume and pushes every real point one or two octree levels deeper than necessary. And the NAVD88 heights are converted to heights above the WGS84 ellipsoid: the output CRS is UTM 18N with no vertical component, which PROJ treats as ellipsoidal. The converter’s transformation to EPSG:4978 then starts from the height it assumes. Skip this and the whole tileset floats or sinks by the geoid separation — about 33 m below where it should be in New York.

Why heights must be ellipsoidal before conversion A cross-section showing the geoid about 33 metres below the WGS84 ellipsoid in New York. A LiDAR point with an orthometric NAVD88 height of 12 metres lies 12 metres above the geoid, which is minus 21 metres relative to the ellipsoid. If the converter treats 12 as an ellipsoidal height, the point is placed 33 metres too high. WGS84 ellipsoid, h = 0 geoid, NAVD88 H = 0 N ≈ −33 m true point: H = 12 m, h = −21 m H read as h: 33 m too high not to scale · the geoid lies below the ellipsoid along the US East Coast
py3dtiles transforms whatever z it receives as an ellipsoidal height. Converting in PDAL, where the vertical CRS is explicit, is the only place the geoid is applied.

2. Confirm every prepared file shares one CRS and sensible bounds

python
import laspy
import numpy as np

mins, maxs, crss = [], [], set()
for f in sorted(PREP.glob("*.laz")):
    with laspy.open(f) as r:
        h = r.header
        crss.add(h.parse_crs().to_epsg())
        mins.append(h.mins); maxs.append(h.maxs)
mins, maxs = np.min(mins, axis=0), np.max(maxs, axis=0)
print("CRS:", crss, "| extent (m):", (maxs - mins).round(1), "| z range:", mins[2].round(1), maxs[2].round(1))
assert crss == {32618}, "all inputs must share one horizontal CRS"
assert maxs[2] - mins[2] < 600, "z range too large: noise survived filtering"

The z-range assertion is tuned to the site. A coastal city with buildings under 450 m has no business spanning 600 m vertically; a mountain valley would need a larger limit. The point is to have one, because the bounding volume of the root tile is exactly this extent.

3. Run the conversion

python
import os
import subprocess

jobs = max(1, os.cpu_count() - 2)
cmd = [
    "py3dtiles", "convert", *[str(p) for p in sorted(PREP.glob("*.laz"))],
    "--out", "tiles/harbour_pc",
    "--overwrite",
    "--srs_in", "32618",
    "--srs_out", "4978",
    "--jobs", str(jobs),
    "--cache_size", "6000",
]
print(" ".join(cmd[:4]), "…")
subprocess.run(cmd, check=True)

--srs_out 4978 makes py3dtiles transform points to Earth-centred coordinates and store them relative to per-tile centres, so the viewer needs no further georeferencing. --jobs sets the number of worker processes; leave a core or two for the main process that reads and dispatches points, or it becomes the bottleneck. --cache_size is in megabytes and bounds how many points are held in memory before nodes are flushed to disk — raise it on a machine with plenty of RAM and the conversion does far less intermediate I/O.

From raw LAZ to a validated point cloud tileset Raw LAZ files in EPSG:32618 plus NAVD88 go through PDAL, which removes noise classes and converts heights to ellipsoidal. A CRS and bounds check follows. py3dtiles converts to 3D Tiles in EPSG:4978 with a pnts file per octree node. A point-count audit compares the sum of points in all pnts files with the prepared input before publishing. raw LAZ32618+5703 PDAL: noise out,h ellipsoidal CRS + boundsassertions py3dtiles→ EPSG:4978 point-countaudit Only the fourth box is the converter; the other four decide whether its output is right.
The conversion itself is one command. Preparation before it and an audit after it are what make the tileset trustworthy.

4. Audit the output: every point accounted for

python
import struct

def pnts_points(path):
    with open(path, "rb") as f:
        header = f.read(28)
        magic, version, byte_len, ft_json_len, ft_bin_len, bt_json_len, bt_bin_len = struct.unpack("<4s6I", header)
        assert magic == b"pnts", f"{path}: not a pnts file"
        feature_table = json.loads(f.read(ft_json_len).rstrip(b" \x00"))
    return feature_table["POINTS_LENGTH"], "RTC_CENTER" in feature_table

out = Path("tiles/harbour_pc")
tiles = list(out.rglob("*.pnts"))
counts = [pnts_points(p) for p in tiles]
total = sum(c for c, _ in counts)

prepared = 0
for f in PREP.glob("*.laz"):
    with laspy.open(f) as r:
        prepared += r.header.point_count

tileset = json.loads((out / "tileset.json").read_text())
print(f"{len(tiles):,} pnts files, {total:,} points in tiles, {prepared:,} prepared, "
      f"refine={tileset['root'].get('refine')}, all RTC: {all(r for _, r in counts)}")
assert total == prepared, "points were lost or duplicated during conversion"

A .pnts file starts with a 28-byte header — the magic, a version and five byte lengths — followed by a JSON feature table whose POINTS_LENGTH gives the count. py3dtiles builds an additive octree: every input point is written exactly once, in one node, with coarse nodes holding a sparse subset and finer nodes the remainder. That makes the audit exact. A total below the prepared count means a worker failed quietly or the disk filled; above it means an input file was listed twice.

Expected Output & Verification

text
h_2026_0412.laz: 18,402,110 points written
CRS: {32618} | extent (m): [2488.4 2490.1  311.7] | z range: -38.2 273.5
14,862 pnts files, 212,730,556 points in tiles, 212,730,556 prepared, refine=ADD, all RTC: True

The negative minimum z is the check that the datum conversion ran: ground at a few metres above NAVD88 near the harbour becomes around −30 m above the ellipsoid. Then verify placement against something independent. Pick a surveyed control point — a benchmark or a corner of a known building — and confirm its ECEF position, computed with pyproj from EPSG:32618+5703, lies within a few centimetres of the nearest point in the tile that contains it after applying that tile’s RTC_CENTER.

Points per octree level in an additive tileset Stacked horizontal bars show where points end up by octree level. The root and first levels hold a small, evenly spaced subset, and each deeper level holds more, with most points in the deepest few levels. Because refinement is additive, the sum across all levels equals the input count exactly. level 0 · 0.2 M level 2 · 1.4 M level 4 · 9 M level 6 · 44 M level 8 · 96 M level 10+ · 62 M sum of all levels = 212.7 M = input count, because each point is stored once
Additive refinement is what makes an exact count audit possible: coarse levels are samples of the data, not copies of it.

Common Errors

The tileset renders 33 m in the air, or under the terrain. Heights reached py3dtiles as NAVD88 orthometric values. Confirm the PDAL step ran with the compound in_srs and that PROJ found the GEOID18 grid — projinfo -s EPSG:32618+5703 -t EPSG:4979 should list a GEOID18 operation, not a ballpark one.

pyproj.exceptions.CRSError or points placed at the wrong zone. Some inputs carry no CRS in their header and py3dtiles was not told --srs_in. Always pass it explicitly, and assert header CRSs in step 2 so missing ones fail before a two-hour conversion.

Conversion slows to a crawl halfway through. The cache filled and workers are flushing to a slow disk. Raise --cache_size, reduce --jobs so the dispatcher keeps up, or put the output directory on local SSD rather than network storage.

Frequently Asked Questions

Should I use py3dtiles or PDAL’s writers.copc for streaming?

They target different viewers. COPC is a single LAZ file streamed by byte range, read by Potree, QGIS and web COPC viewers; 3D Tiles is what Cesium-family runtimes and many twin platforms load natively. Pipelines serving both often write both from the same prepared input.

How do I keep classification and intensity?

py3dtiles writes classification and intensity into the batch table by default in recent versions; check the flags on your installed version with py3dtiles convert --help, and confirm by reading a .pnts batch table JSON after conversion.

Can I convert tile by tile and merge later?

Yes, and for very large collections it is easier to operate. Convert each shard to its own tileset and combine them as described in merging shard tilesets into a root tileset. Expect slightly more points at shard edges to be visible at coarse levels, because each shard builds its own sample.

How long should a city-scale conversion take?

Throughput is dominated by LAZ decompression and disk writes, not by the tree building. As a planning figure, a 16-core machine with local NVMe converts in the region of 20 to 40 million points per minute from prepared LAZ; a 200-million-point district is a ten-minute job, and a whole city of several billion points is an overnight one. If a run is an order of magnitude slower than that, look at the storage first — network file systems and container overlay filesystems are the usual culprits.

Does the output need Draco or other compression?

Point cloud tiles compress well, and Cesium supports Draco-compressed .pnts through the 3DTILES_draco_point_compression extension. py3dtiles writes uncompressed tiles, which are simple to audit; compress in a separate post-processing step once the count audit has passed, so the audit always runs on tiles you can read with twenty lines of Python.

Back to Automated Tile Generation for 3D Geospatial.