Setting a Vertical CRS on a Point Cloud

TL;DR: Write a compound CRS with a_srs on the writer — "EPSG:6339+5703" — so the file records what its Z values mean. A horizontal-only code leaves the vertical datum undeclared, and every downstream tool is then free to assume something different.

# Context and Motivation

This guide belongs to Coordinate Reference Systems, which covers where a CRS lives in a LAS file and how PDAL resolves it. This page is about the half of the CRS that most deliveries omit entirely.

A projected EPSG code describes a horizontal system and says nothing about height. Points in a file tagged EPSG:6339 carry a Z value that might be an ellipsoidal height, an orthometric height above NAVD88, or a local site datum — and nothing in the file distinguishes them. The consequence is not theoretical: two tiles from the same campaign, one delivered in ellipsoidal heights and one in orthometric, differ by around thirty metres and merge into a terrain model with a cliff down the middle. Declaring the vertical datum costs one option and removes the entire class of failure.

# Prerequisites and Assumptions

# Step-by-Step Implementation

# Step 1 — Confirm the compound code exists

projinfo "EPSG:6339+5703"

The output names both components. An error here means the pairing is not defined and you need a WKT string instead.

# Step 2 — Write it with a_srs

{"type": "writers.las", "filename": "tile.laz", "compression": "laszip",
 "minor_version": 4, "dataformat_id": 6,
 "a_srs": "EPSG:6339+5703", "forward": "all"}

a_srs relabels; it does not transform. Use it when the coordinates are already correct and only the declaration is missing. When the heights must actually move, that is a vertical datum transform instead.

# Step 3 — Prefer WKT when the EPSG pairing does not exist

{"a_srs": "COMPD_CS[\"NAD83(2011) / UTM 11N + NAVD88\", ...]"}

Verbose, but it is the only way to express many national and site datums.

# Step 4 — Read it back and confirm both components

pdal info tile.laz --metadata | grep -i -A6 "compound\|vertical"

# Complete Working Example

"""Stamp a compound CRS onto a directory of tiles and verify both components."""
from __future__ import annotations

import json
import logging
from pathlib import Path

import pdal

LOG = logging.getLogger("vcrs")


def stamp(src: Path, dst: Path, compound: str) -> int:
    spec = json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(src)},
        {"type": "writers.las", "filename": str(dst), "compression": "laszip",
         "minor_version": 4, "dataformat_id": 6,
         "a_srs": compound, "forward": "all"},
    ]})
    return pdal.Pipeline(spec).execute()


def declared_srs(path: Path) -> str:
    p = pdal.Pipeline(json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(path), "count": 1}]}))
    p.execute()
    srs = p.quickinfo["readers.las"].get("srs", {})
    return srs.get("wkt", "") or srs.get("horizontal", "")


def run(in_dir: Path, out_dir: Path, compound: str, vertical_name: str) -> list[dict]:
    out_dir.mkdir(parents=True, exist_ok=True)
    results = []
    for src in sorted(in_dir.glob("*.laz")):
        dst = out_dir / src.name
        n = stamp(src, dst, compound)
        wkt = declared_srs(dst)
        if vertical_name.lower() not in wkt.lower():
            raise AssertionError(
                f"{dst.name}: vertical datum {vertical_name!r} is absent from the written CRS"
            )
        LOG.info("%s — %d points, compound CRS recorded", dst.name, n)
        results.append({"tile": dst.name, "points": n})
    return results


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    print(json.dumps(run(Path("tiles"), Path("tiles_vcrs"),
                         compound="EPSG:6339+5703",
                         vertical_name="NAVD88"), indent=2))

# Key Parameter Table

# Verification

The vertical component is in the written WKT. Asserted in the example by name, which is cruder than parsing and considerably harder to fool.

Coordinates did not move. a_srs must not change a single value. Compare min and max Z before and after; any difference means a transform crept in.

Downstream tools agree. Open the file in GDAL or a desktop GIS and confirm it reports a compound system. A tool that shows only the horizontal part is telling you the record was written as GeoTIFF keys rather than WKT.

# Gotchas and Edge Cases

GeoTIFF keys cannot really express a compound CRS. Write LAS 1.4 with a WKT record. A 1.2 file with GeoTIFF keys will lose the vertical component whatever you pass to a_srs.

a_srs on a file whose heights are actually ellipsoidal makes matters worse. You have now confidently mislabelled the data. Confirm what Z means before declaring it.

Vertical units are separate from horizontal units. US state plane in survey feet with NAVD88 heights in metres is common and legal. Check both axes of the compound definition.

Merging still needs matching datums. Declaring two tiles honestly as NAVD88 and ellipsoidal does not let them merge; it lets you notice that they cannot, which is the entire benefit.

# Frequently Asked Questions

Why does a projected EPSG code not describe my heights?

Because it describes only the horizontal system. A file tagged EPSG:6339 carries Z values that could be ellipsoidal heights, orthometric heights above NAVD88, or a local site datum, and nothing in the file distinguishes them. A compound code such as EPSG:6339+5703 states which.

Does a_srs change my coordinates?

No. It relabels the file, recording a different declaration for the same numbers. That is the right tool when the coordinates are correct and the declaration was missing or wrong. When the heights themselves must move, you need a reprojection with compound CRSs on both sides.

Why did the vertical part disappear from my output?

Almost certainly because the file was written as LAS 1.2, or with GeoTIFF keys rather than an OGC WKT record. Neither can express a compound coordinate system properly, so the vertical component is dropped without an error.

What if I do not know the vertical datum?

Find out before declaring anything. Guessing between ellipsoidal and orthometric heights is a thirty-metre error, and a confidently wrong declaration is worse than an absent one because it stops anyone downstream from asking.

# Related

• Coordinate Reference Systems — the parent guide to where a CRS lives in a LAS file
• Fixing CRS Mismatches in Point Clouds — diagnosing an offset before deciding it is vertical
• Handling Vertical Datum Transforms in PDAL — when the heights have to move rather than be relabelled
• Metadata and Header Sync — keeping the header honest through a pipeline
• Point Cloud Data Standards and Fundamentals — the section overview