Point Cloud Formats and Interchange in Python

A point cloud is the easiest product in a photogrammetry pipeline to get subtly wrong on the way out of the door. The file opens, renders, and looks right. Somebody measures a distance in it and gets an answer two centimetres from the one the pipeline computed, or opens it in a second application and finds it positioned at the origin, or discovers that the classification the pipeline spent an hour producing is not in the file at all.

Every one of those is a format decision. LAS quantises coordinates through a header field most tools set to a default. Coordinate reference systems are carried in several mutually incompatible ways depending on the format version. Extra per-point attributes survive or vanish depending on flags. None of it produces an error.

This page covers writing point clouds that carry their full information to somebody else’s software: the integer encoding underneath LAS, the compression choices, cloud-optimized formats for streaming large clouds to a browser, and a validation routine that reads back what was written.

Audience and prerequisites. Python 3.10+, PDAL 2.5+ or laspy 2.5+, and a processed cloud in a projected CRS. Familiarity with the classification stage in classifying point clouds with PDAL and Python is assumed, because the attributes produced there are the ones most often lost on export.

# Prerequisites

# Conceptual architecture

LAS does not store coordinates as floating-point numbers. It stores three 32-bit integers per point, and the header carries a scale and an offset for each axis. The real coordinate is reconstructed as value × scale + offset.

Two consequences follow, and both bite. The scale is the quantisation step: with the widespread default of 0.01, every coordinate is rounded to the nearest centimetre, which discards precision a survey may have worked hard to achieve. And the signed 32-bit range is about ±2.1 billion, so at a scale of 0.001 the representable span is ±2,147 km from the offset — comfortable, until the offset is left at zero and a UTM northing of 5.6 million metres requires 5.6 billion counts, which overflows.

Figure 1 — Two header fields, three outcomes. Only the third is a survey deliverable.

# Step 1: Derive the header from the data

import numpy as np


def header_geometry(xyz: np.ndarray, target_precision: float = 0.001) -> dict:
    """Scale and offset that preserve `target_precision` without overflow.

    The offset is placed at the data's own centre, rounded to a whole metre so
    it reads cleanly in the header and is exactly representable.
    """
    lo, hi = xyz.min(axis=0), xyz.max(axis=0)
    offset = np.round((lo + hi) / 2.0)

    span = np.abs(np.vstack([lo - offset, hi - offset])).max(axis=0)
    min_scale = float((span / 2_147_483_000).max())
    if min_scale > target_precision:
        raise ValueError(
            f"extent requires scale {min_scale:.6f} m, coarser than the requested "
            f"{target_precision} m — tile the cloud before writing")

    return {"scale": (target_precision,) * 3,
            "offset": tuple(float(v) for v in offset)}

Raising rather than silently coarsening is the design decision that matters. A pipeline that quietly drops to centimetre precision on a large site produces files that open, render and measure wrong — the worst possible failure mode, because nothing signals it.

# Step 2: Write the CRS in a way the reader will find

LAS 1.2 carries a CRS as GeoTIFF keys; LAS 1.4 carries it as a WKT string, and only if the global encoding bit says so. Software differs in which it looks for. Writing 1.4 with WKT is the modern default, but a client on older software may see no CRS at all — and a point cloud without a declared CRS opens at the origin.

import laspy
import numpy as np
from pyproj import CRS


def write_las(xyz: np.ndarray, classification: np.ndarray, rgb: np.ndarray,
              crs_epsg: int, out_path: str, *, precision: float = 0.001) -> None:
    """Write LAS 1.4 with WKT CRS, correct scaling and preserved attributes."""
    geom = header_geometry(xyz, precision)
    header = laspy.LasHeader(point_format=7, version="1.4")  # 7 = RGB + time
    header.scales = np.array(geom["scale"])
    header.offsets = np.array(geom["offset"])
    header.add_crs(CRS.from_epsg(crs_epsg))

    las = laspy.LasData(header)
    las.x, las.y, las.z = xyz[:, 0], xyz[:, 1], xyz[:, 2]
    las.classification = classification.astype(np.uint8)
    # LAS stores colour as 16-bit; 8-bit source values must be scaled, or the
    # cloud renders almost black in any viewer that reads the spec correctly.
    las.red, las.green, las.blue = (rgb[:, i].astype(np.uint16) * 257
                                    for i in range(3))
    las.write(out_path)

The * 257 on the colour channels is the most frequently missed line in point-cloud export. LAS colour is unsigned 16-bit; writing 8-bit values directly produces a cloud whose brightest pixel is at 0.4 % of full scale. Some viewers auto-stretch and hide it, which is why the bug survives so long.

# Step 3: Choose the compression and container

Uncompressed LAS is roughly 34 bytes per point with colour. LAZ compresses losslessly to around a quarter of that, costs a few seconds per hundred million points, and is read natively by every relevant tool. There is no case for shipping uncompressed LAS.

COPC — cloud-optimized point cloud — is LAZ with an octree index built into the file, which lets a client fetch only the parts of the cloud they are looking at over HTTP range requests. A 40 GB survey becomes something a browser can open, with no server beyond static hosting.

import json
import subprocess


def write_copc(src: str, dst: str) -> None:
    """Convert to COPC so a client can stream it without downloading it."""
    pipeline = {"pipeline": [
        src,
        {"type": "writers.copc", "filename": dst, "forward": "all"},
    ]}
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)

The forward: all again carries header fields through. Without it the writer applies its own defaults for scale, offset and CRS, and a correctly written LAS becomes an incorrectly written COPC. Converting LAS to COPC for cloud streaming covers the hosting requirements that go with it.

# Step 4: Keep the attributes the pipeline produced

Standard LAS point formats carry a fixed set of fields. Anything else — the excess-green index, a planarity value, a per-point uncertainty — must be declared as an extra dimension, and extra dimensions are dropped by any stage that does not know about them.

import json
import subprocess


def convert_preserving_extras(src: str, dst: str,
                              extras: dict[str, str]) -> None:
    """Convert between formats without losing declared extra dimensions.

    `extras` maps a dimension name to its type, e.g. {"Planarity": "float"}.
    Every stage in the pipeline must be told about them; a writer that is not
    simply omits them, silently and without warning.
    """
    spec = ",".join(f"{name}={dtype}" for name, dtype in extras.items())
    pipeline = {"pipeline": [
        src,
        {"type": "writers.las", "filename": dst, "compression": "laszip",
         "forward": "all", "extra_dims": spec or "all"},
    ]}
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)

extra_dims: "all" is the safe default for internal conversions. For a delivered file, naming the dimensions explicitly is better: it documents what the client is receiving, and it keeps intermediate scratch fields out of a product.

Figure 2 — What a conversion costs, by format. PLY is convenient for meshing tools and is not a delivery format for georeferenced data.

# Step 5: Tile a survey that will not fit in one file

A single LAS file has no hard size limit in the 1.4 specification, but practical limits arrive long before the format’s do. A 300-million-point cloud is around 3 GB compressed, which most desktop software will attempt to load entirely into memory and many will fail on. A survey delivered as one enormous file is a survey the client cannot open.

Tiling solves that, and the tiling scheme is worth choosing deliberately rather than accepting a default. Three properties matter.

Tiles should be square and aligned to a fixed grid. A grid anchored to round coordinates — 500 m tiles starting at a multiple of 500 — means the same tile boundaries recur in every survey of the site, so a monitoring series can compare tile against tile without re-tiling anything.

Tile names should encode position. A name such as site_e502000_n5600500.laz tells a human and a script where the tile is without opening it. Sequential numbering does not, and turns every “which tile covers the weighbridge” question into a search.

Tiles must not overlap in the delivered product. A buffer is essential during processing, as described for classification, and must be removed before delivery — or a client merging the tiles gets duplicate points, and any per-point statistic they compute is subtly wrong.

import json
import subprocess
from pathlib import Path


def tile_for_delivery(src: str, out_dir: str, *, tile: float = 500.0,
                      origin: tuple[float, float] = (0.0, 0.0)) -> None:
    """Split a survey into non-overlapping, grid-aligned, position-named tiles.

    filters.splitter with origin_x/origin_y anchors the grid to fixed
    coordinates rather than to the data's own corner, so the same tile
    boundaries recur across every survey of the site.
    """
    Path(out_dir).mkdir(parents=True, exist_ok=True)
    pipeline = {"pipeline": [
        src,
        {"type": "filters.splitter", "length": tile,
         "origin_x": origin[0], "origin_y": origin[1]},
        {"type": "writers.las", "filename": f"{out_dir}/tile_#.laz",
         "compression": "laszip", "forward": "all"},
    ]}
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)


def rename_tiles_by_position(out_dir: str, prefix: str) -> list[str]:
    """Rename sequentially numbered tiles to position-encoded names."""
    import laspy
    renamed = []
    for path in sorted(Path(out_dir).glob("tile_*.laz")):
        with laspy.open(path) as fh:
            h = fh.header
            name = (f"{prefix}_e{int(h.mins[0]) // 1 :d}"
                    f"_n{int(h.mins[1]) // 1 :d}.laz")
        target = path.with_name(name)
        path.rename(target)
        renamed.append(str(target))
    return renamed

Alongside the tiles, deliver an index: a small GeoJSON with one polygon per tile carrying its filename, point count and bounds. Clients load it in any GIS, see the coverage at a glance, and open only the tiles they need. It costs one function and removes most of the support traffic a tiled delivery otherwise generates.

import json
from pathlib import Path

import laspy


def write_tile_index(out_dir: str, index_path: str, crs_epsg: int) -> None:
    """A GeoJSON footprint per tile, so a client can see coverage before downloading."""
    features = []
    for path in sorted(Path(out_dir).glob("*.laz")):
        with laspy.open(path) as fh:
            h = fh.header
            x0, y0 = float(h.mins[0]), float(h.mins[1])
            x1, y1 = float(h.maxs[0]), float(h.maxs[1])
            features.append({
                "type": "Feature",
                "properties": {"file": path.name, "points": int(h.point_count),
                               "z_min": float(h.mins[2]), "z_max": float(h.maxs[2])},
                "geometry": {"type": "Polygon", "coordinates": [[
                    [x0, y0], [x1, y0], [x1, y1], [x0, y1], [x0, y0]]]},
            })
    Path(index_path).write_text(json.dumps({
        "type": "FeatureCollection",
        "crs": {"type": "name", "properties": {"name": f"EPSG:{crs_epsg}"}},
        "features": features}, indent=2))

Where the client’s tooling can stream, a single COPC is better than any tiling scheme — one file, no index, and the viewer fetches only what it needs. Tiling remains the right answer for desktop GIS workflows and for anyone who wants the data on a disk rather than behind a URL, which is still most survey clients.

# Parameter deep-dive

# Verification and output inspection

Everything above is checkable by reading the file back and comparing it against the arrays that were written.

import laspy
import numpy as np
from pyproj import CRS


def verify_written(path: str, xyz: np.ndarray, expect_epsg: int,
                   *, tol: float = 0.0015) -> dict:
    """Round-trip check: does the file contain what we meant to write?"""
    las = laspy.read(path)
    back = np.column_stack([las.x, las.y, las.z])

    problems = []
    if back.shape != xyz.shape:
        problems.append(f"point count changed: {xyz.shape[0]} → {back.shape[0]}")
    else:
        worst = float(np.abs(back - xyz).max())
        if worst > tol:
            problems.append(f"coordinates differ by up to {worst*1000:.1f} mm "
                            "— the header scale is too coarse")

    crs = las.header.parse_crs()
    if crs is None:
        problems.append("no CRS in the file — it will open at the origin")
    elif CRS.from_user_input(crs).to_epsg() != expect_epsg:
        problems.append(f"CRS is {CRS.from_user_input(crs).to_epsg()}, "
                        f"expected {expect_epsg}")

    if hasattr(las, "red") and int(np.asarray(las.red).max()) < 4096:
        problems.append("colour maximum is very low — 8-bit values written "
                        "into a 16-bit field")

    return {"points": int(len(back)), "scales": las.header.scales.tolist(),
            "offsets": las.header.offsets.tolist(), "problems": problems}

Running this on every delivered file, in the pipeline rather than by hand, is the single highest-value check in this section. It takes seconds and catches all four of the failures described above.

The tolerance deserves a word. Comparing read-back coordinates against the source at a tolerance slightly above the intended scale — 1.5 mm for a 1 mm scale — allows for the rounding the encoding legitimately performs while still catching a header that is an order of magnitude too coarse. Setting the tolerance below the scale makes the check fail on every correctly written file; setting it far above makes it pass on files that have thrown away most of the survey’s precision.

Two further checks are worth adding for delivered products. Compare the point count in the header against the number of points actually readable, because a truncated write produces a file whose header claims more than it contains and which many readers accept silently. And confirm that every extra dimension the pipeline declared is present in the output, since a dimension dropped by an intermediate stage is invisible until somebody tries to use it months later.

Figure 3 — The requirement picks the format; flexibility is a liability in a deliverable.

# Troubleshooting

The cloud opens at the origin in one application and correctly in another. The CRS is present in one of the two encodings the reader does not support. Write LAS 1.4 with WKT, and if the client is on older software, also supply a .prj sidecar.

Coordinates are out by a few millimetres against the source. The header scale is coarser than the data. Set it from the survey precision and re-write; there is no way to recover the lost digits from the file.

The cloud renders almost black. Eight-bit colour written into LAS’s sixteen-bit fields. Multiply by 257 on write.

Points appear scattered across the globe. Integer overflow from an offset of zero with a fine scale. Set the offset near the data’s centre.

The classification is gone after a conversion. The target point format does not carry it, or an intermediate stage dropped it. Check the point format and use forward: all.

A COPC file streams the wrong area. The file’s own bounds disagree with its index, usually because it was written from a source whose header was wrong. Re-write from the corrected LAS rather than patching the COPC.

# Related

• Converting LAS to COPC for Cloud Streaming
• Fixing LAS Scale and Offset Precision Loss
• Reading E57 and PLY into a Python Pipeline
• Fixing SRS Lost When Writing LAS Files
• Exporting Cloud-Optimized GeoTIFFs with Rasterio

← Point Cloud Processing & 3D Deliverables