Python LiDAR

Handling Vertical Datum Transforms in PDAL

TL;DR: A vertical transform needs a compound CRS on both sides — EPSG:6339+5703, not EPSG:6339 — plus the geoid grid installed where PROJ can find it. Without both, filters.reprojection moves your coordinates horizontally, leaves Z untouched, and reports success.

# Context and Motivation

This guide is part of Spatial Reprojection in PDAL, which covers the horizontal case in depth. Vertical transformation is the half that fails silently, and it fails in a way that surveys notice months later: a terrain model that fits the control network in plan and sits thirty metres off in height.

The root of the problem is that a plain projected CRS code says nothing about the vertical axis. EPSG:6339 is NAD83(2011) / UTM zone 11N — a horizontal system. Points in a LAS file tagged with it carry a Z value, but the file does not record whether that Z is a height above the ellipsoid, above a geoid model, or above a tide-gauge datum from 1929. PDAL cannot transform what has not been declared, so it does the only safe thing: it transforms X and Y and passes Z through unchanged. That is correct behaviour and it is almost never what the user wanted.

The thirty metres a vertical datum accounts for A cross-section with three surfaces. The ellipsoid is a smooth mathematical reference. The geoid undulates above and below it — the separation is about minus thirty metres in this example. The terrain sits above the geoid. An ellipsoidal height measured to the terrain differs from the orthometric height by exactly the geoid separation at that point, and the separation changes across the tile. ellipsoid — the mathematical reference surface geoid — where mean sea level would sit terrain ellipsoidal height h orthometric height H separation N ≈ −30 m H = h − N, and N varies across a single tile — which is why a constant offset is not a substitute for a geoid grid

# Prerequisites and Assumptions

Requirement Detail
PDAL 2.4+ built against PROJ 8 or later
PROJ data grids proj-data package, or PROJ_NETWORK=ON to fetch grids on demand
A declared source vertical datum from the survey report; the file will rarely tell you
Compound EPSG codes e.g. EPSG:6339+5703 for NAD83(2011) UTM 11N with NAVD88 heights
projinfo to confirm the transformation exists before running the pipeline

If nobody can tell you what the source Z means, stop. Guessing between ellipsoidal and orthometric heights is a thirty-metre coin flip, and no amount of pipeline configuration recovers from choosing wrong.

# Step-by-Step Implementation

# Step 1 — Confirm the grid is installed

bash
projinfo -s "EPSG:6339+5703" -t "EPSG:6318+5703" --spatial-test intersects -o PROJ

The output lists candidate operations with accuracies. If the only candidate is a null transformation, or the output mentions a missing grid file, the vertical shift will not happen.

# Step 2 — Declare both sides as compound CRSs

json
{
  "type": "filters.reprojection",
  "in_srs": "EPSG:6339+5703",
  "out_srs": "EPSG:6318+4979"
}

5703 is NAVD88 height; 4979 is WGS84 ellipsoidal height. Writing EPSG:6339 alone on either side disables the vertical part of the transform entirely.

# Step 3 — Set the writer’s Z scale deliberately

Vertical units may not survive the transform unchanged. If the source was in US survey feet and the target is metres, scale_z on the writer must be reconsidered — the header does not follow automatically, as metadata and header sync explains.

# Step 4 — Check a known point

Take a control point with a published orthometric height, run it through the same transform, and compare. A single point is enough to catch a missing grid, and nothing else will.

# Complete Working Example

python
"""Transform a tile between compound CRSs and verify the vertical shift happened."""
from __future__ import annotations

import json
import logging
from pathlib import Path

import numpy as np
import pdal

LOG = logging.getLogger("vdatum")


def spec(src: Path, dst: Path, in_srs: str, out_srs: str) -> str:
    return json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(src)},
        {"type": "filters.reprojection", "in_srs": in_srs, "out_srs": out_srs},
        {"type": "writers.las", "filename": str(dst), "compression": "laszip",
         "scale_z": 0.001, "offset_z": "auto", "forward": "all"},
    ]})


def z_stats(path: Path) -> tuple[float, float]:
    p = pdal.Pipeline(json.dumps({"pipeline": [{"type": "readers.las", "filename": str(path)}]}))
    p.execute()
    z = p.arrays[0]["Z"]
    return float(np.min(z)), float(np.max(z))


def run(src: Path, dst: Path, in_srs: str, out_srs: str, expect_shift: float,
        tolerance: float = 2.0) -> None:
    before_min, before_max = z_stats(src)
    pdal.Pipeline(spec(src, dst, in_srs, out_srs)).execute()
    after_min, after_max = z_stats(dst)

    observed = ((after_min + after_max) / 2) - ((before_min + before_max) / 2)
    LOG.info("mean Z moved by %.3f m (expected about %.3f m)", observed, expect_shift)

    if abs(observed) < 0.01:
        raise AssertionError(
            "Z did not move at all — the vertical transform was not applied. "
            "Check that both CRSs are compound and that the geoid grid is installed."
        )
    if abs(observed - expect_shift) > tolerance:
        raise AssertionError(
            f"vertical shift {observed:.3f} m is not the expected {expect_shift:.3f} m"
        )


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    run(Path("tile_navd88.laz"), Path("tile_ellipsoidal.laz"),
        in_srs="EPSG:6339+5703", out_srs="EPSG:6339+4979", expect_shift=31.4)
What projinfo shows before you run anything Three candidate operations between the same pair of compound CRSs. The grid-based operation is accurate to two centimetres but only if the grid files are installed. The Helmert fit is always available at half a metre. The null transform applies no vertical shift at all, and PROJ will silently use it when the grid is missing. NADCON5 + GEOID18 grids 0.02 m installed Helmert 7-parameter 0.50 m always available null vertical transform no shift at all the silent fallback projinfo -s EPSG:6339+5703 -t EPSG:6339+4979 — candidates, best first PROJ picks the best candidate it can actually execute, and reports success either way

# Key Parameter Table

Setting Where Purpose
in_srs filters.reprojection Overrides whatever the file declares; required when the file omits the vertical part
out_srs filters.reprojection Must include a vertical code for Z to move at all
PROJ_NETWORK environment ON lets PROJ fetch missing grids over the network at run time
PROJ_DATA environment Directory holding the installed grids; wrong value means silent null transforms
scale_z / offset_z writers.las Vertical precision after the transform; auto offset avoids integer range problems

# Verification

Z actually moved. The assertion in the example: if the mean elevation shifted by less than a centimetre, no vertical transform was applied.

The shift matches the geoid model. Compare against the published separation for the tile centre. Regional geoid separations range from roughly −100 m to +85 m worldwide, so a shift in the wrong direction is as diagnostic as no shift at all.

Control points agree. The only check that matters to a surveyor. One point with a published orthometric height, transformed, compared.

# Gotchas and Edge Cases

Four configurations, one that works Four combinations of source and target CRS declarations. Horizontal-only codes on both sides move X and Y only. A compound source with a horizontal target still moves nothing vertically. A horizontal source with a compound target fails or assumes ellipsoidal. Only compound codes on both sides, with the geoid grid present, transform the height. what you declare what happens to Z EPSG:6339 → EPSG:6318 unchanged — no vertical axis declared EPSG:6339+5703 → EPSG:6318 unchanged — the target has no vertical EPSG:6339+5703 → EPSG:6318+4979 moves — if the geoid grid is installed the same, with PROJ_DATA unset unchanged — grid missing, null transform

The null transform is not an error. PROJ prefers a lower-accuracy operation over failing, so a missing grid produces a run that succeeds and does nothing. The assertion in the example exists precisely because the exit code will not tell you.

Vertical units and horizontal units differ more often than you expect. US state plane systems in survey feet with NAVD88 heights in metres are common. Check both axes of the compound CRS, not just the horizontal one.

Containers lose grids. An image that works on your laptop may not carry proj-data; the Docker containers guide covers pinning the whole native stack so this cannot drift.

# Frequently Asked Questions

Why did my reprojection leave Z unchanged?

Almost certainly because one of the two CRSs had no vertical component. A code like EPSG:6318 describes only the horizontal system, so PDAL transforms X and Y and passes Z through. Use compound codes on both sides, such as EPSG:6339+5703 to EPSG:6318+4979.

How do I know the geoid grid is installed?

Run projinfo between the source and target codes and read the candidate operations. If the only candidate is a null transformation, or the output names a missing grid file, the shift will not happen. Setting PROJ_NETWORK=ON lets PROJ fetch the grid at run time instead.

Can I just add a constant offset instead?

Only if the tile is small and the accuracy requirement is loose. Geoid separation varies across the ground, often by more than a metre within a kilometre, so a constant offset reproduces the mean and leaves a tilted error surface behind.

What if nobody knows what the source Z means?

Stop and find out. The difference between ellipsoidal and orthometric heights is tens of metres, and no pipeline configuration recovers from guessing wrong. The survey report, the acquisition contract or the data provider will say; the LAS header usually will not.

Related in Spatial Reprojection in PDAL: Coordinate Transformations for LiDAR Pipelines

↑ Spatial Reprojection in PDAL: Coordinate Transformations for LiDAR Pipelines overview --- title: "Spatial Reprojection in PDAL: Coordinate Transformations for LiDAR Pipelines" description: "How to reproject LAS/LAZ point... Reprojecting Point Clouds from UTM to WGS84 with PDAL and Python --- title: "Reprojecting Point Clouds from UTM to WGS84 with PDAL and Python" description: "Step-by-step guide to reprojecting LAS/LAZ point...