Inspecting PROJ Transformations Before Reprojecting

TL;DR: Run projinfo -s <src> -t <dst> --spatial-test intersects -o PROJ for your area to list every candidate operation with its accuracy and the grids it needs. If the best one needs a grid you do not have, either enable PROJ_NETWORK=ON or install the grid with projsync. For reproducible production runs, compare PDAL’s output against the expected operation with pyproj, and pin the environment so the same grids exist everywhere.

# Context and Motivation

This guide is part of Spatial Reprojection. filters.reprojection delegates to PROJ, and PROJ chooses among several possible coordinate operations between two CRSs — a grid-based datum shift accurate to centimetres, a Helmert transform accurate to a metre, or a “ballpark” transformation that ignores the datum difference altogether. It picks the most accurate one whose grids are available on that machine. The same pipeline can therefore produce outputs that differ by a metre between a laptop with grids installed and a container without them, and nothing in PDAL’s output tells you which happened.

Inspecting the candidates before a production run, and checking which one was used, turns that silent variability into a documented choice.

# Prerequisites and Assumptions

• PROJ 7+ with the projinfo and projsync utilities (bundled with PROJ in conda-forge and most Linux packages).
• pyproj 3.x in Python, built against the same PROJ as PDAL where possible.
• Source and target CRSs identified, ideally as EPSG codes, with a bounding box of the project in geographic coordinates.

# Step-by-Step Implementation

# Step 1 — List candidate operations for your area

projinfo with --spatial-test intersects and an area of use lists only operations valid where your data is.

projinfo -s EPSG:4269 -t EPSG:6318 --spatial-test intersects \
  --bbox -76.0,39.5,-75.0,40.5 --summary

# Step 2 — Read accuracy and grid requirements

Each candidate shows its accuracy and whether it is usable with installed grids. projinfo ... -o PROJ prints the full pipeline, including +proj=hgridshift +grids=... steps.

# Step 3 — Make the grids available

Either allow PROJ to fetch grids on demand with PROJ_NETWORK=ON (it caches them locally), or install them ahead of time with projsync --bbox ... so containers work offline.

# Step 4 — Check what PDAL will use

Build the same transformation with pyproj’s TransformerGroup, confirm the best candidate is available, and compare a test point against PDAL’s output.

# Step 5 — Pin and record

Record the chosen operation’s description in your processing log or output metadata, and bake the grids into the container image so every worker uses the same one.

# Complete Working Example

"""Report PROJ candidates for a reprojection and verify PDAL used the best one."""
from __future__ import annotations

import json
import os

import numpy as np
import pdal
from pyproj.transformer import TransformerGroup

SRC, DST = "EPSG:6318+5703", "EPSG:6347+5703"      # NAD83(2011) geographic -> UTM 18N
AREA = (-76.0, 39.5, -75.0, 40.5)                  # west, south, east, north


def candidates() -> TransformerGroup:
    from pyproj.aoi import AreaOfInterest
    group = TransformerGroup(SRC, DST, always_xy=True,
                             area_of_interest=AreaOfInterest(*AREA))
    for i, t in enumerate(group.transformers):
        print(f"[{i}] acc={t.accuracy} m  {t.description}")
    for op in group.unavailable_operations:
        grids = [g.short_name for g in op.grids if not g.available]
        print(f"[unavailable] acc={op.accuracy} m  {op.name}  missing grids: {grids}")
    if not group.best_available:
        print("WARNING: the most accurate operation is not available on this machine")
    return group


def pdal_point(x: float, y: float, z: float) -> np.ndarray:
    arr = np.array([(x, y, z)], dtype=[("X", "f8"), ("Y", "f8"), ("Z", "f8")])
    p = pdal.Pipeline(json.dumps({"pipeline": [
        {"type": "filters.reprojection", "in_srs": SRC, "out_srs": DST}]}), arrays=[arr])
    p.execute()
    out = p.arrays[0][0]
    return np.array([out["X"], out["Y"], out["Z"]])


if __name__ == "__main__":
    print("PROJ_NETWORK =", os.environ.get("PROJ_NETWORK", "OFF"))
    group = candidates()
    x, y, z = -75.5, 40.0, 120.0
    best = np.array(group.transformers[0].transform(x, y, z))
    got = pdal_point(x, y, z)
    diff = np.abs(best - got)
    print("pyproj best:", best.round(3), " PDAL:", got.round(3), " diff:", diff.round(3))
    assert diff.max() < 0.01, "PDAL did not use the best available operation"

Run it once with grids and once in your production container. If the outputs or the warning differ, the container is missing grids.

# Key Parameter Table

# Verification

• Best available equals best overall. group.best_available is true on every machine that runs production.
• PDAL matches pyproj. The test point transformed by PDAL equals pyproj’s first candidate to within a millimetre.
• Environment recorded. The output metadata or run log records the PROJ version and the operation description.

# Gotchas and Edge Cases

Ballpark transformations are silent. When no real datum transformation is known between two CRSs, PROJ may use a ballpark operation that treats different datums as identical. projinfo labels it clearly; PDAL does not. Always inspect when the source and target datums differ.

Network grids in batch jobs. PROJ_NETWORK=ON downloads grids on first use, which works on a laptop but can fail or be slow in a fleet of containers without internet egress. Bake grids into the image with projsync during the build instead.

Different PROJ builds for PDAL and pyproj. pip-installed pyproj wheels bundle their own PROJ. If PDAL uses a different PROJ with a different grid set, the check above can disagree for reasons unrelated to your data. Install both from conda-forge to share one PROJ.

EPSG database updates. New PROJ releases add operations and change preferences. Pin the PROJ version in production and re-run the inspection when you upgrade it.

# Frequently Asked Questions

How do I know which transformation filters.reprojection used?

PDAL does not report it directly. Reproduce the transformation with pyproj using the same PROJ installation, list the candidates, and compare a test point: if PDAL matches the first candidate, it used the best available operation.

What does PROJ_NETWORK=ON do?

It lets PROJ download missing grid files from its content delivery network when an operation needs them, caching them locally. It is convenient interactively but should be replaced by pre-installed grids in reproducible batch environments.

Why do I get different results on different machines?

Almost always because different grids are installed, so PROJ selects different operations. Inspect the candidates on both machines and install the same grids everywhere.

Is a Helmert transformation good enough for LiDAR?

It depends on the accuracy you need. Helmert transformations between NAD83 and WGS84 realizations are typically accurate to about a metre, which is far worse than modern LiDAR’s vertical accuracy. Use grid-based operations when they exist.

# Related

• Spatial Reprojection — the reprojection stage in general
• Handling Vertical Datum Transforms in PDAL — geoid grids for heights
• Reprojecting State Plane Feet to Metres — units and compound CRSs
• Building a Slim PDAL Docker Image — baking grids into containers
• Choosing a Projected CRS for a LiDAR Project — picking the target in the first place