Reprojecting Point Clouds from UTM to WGS84 with PDAL and Python

TL;DR: Pass filters.reprojection with in_srs: "EPSG:<utm_zone>" and out_srs: "EPSG:4326" in a PDAL pipeline, then write with writers.las using a_srs: "EPSG:4326" and forward: "all" — the header, all attributes, and Z values in their original vertical datum are preserved automatically.

# Context and Motivation

This guide is part of Spatial Reprojection in PDAL, which covers the full range of coordinate system transformations available in a PDAL pipeline.

UTM coordinates are ubiquitous in survey-grade LiDAR deliverables — they are metre-based, zone-specific, and optimised for regional accuracy. WGS84 (EPSG:4326) is the geographic coordinate system that web maps, cloud platforms, and most interoperability formats expect. The mismatch creates a routine but error-prone handoff: naive reprojection either corrupts the LAS/LAZ file header CRS, silently applies approximate datum shifts, or produces unit confusion when Z values in metres are mixed with X/Y in decimal degrees.

Getting this right matters beyond aesthetics. A mismatched coordinate reference system in a production point cloud can cascade into broken ground classification, incorrect hillshade rasters, and failed API ingest. PDAL’s filters.reprojection stage, backed by the PROJ engine, handles the inverse map projection and datum grid application in a single streaming pass — no in-memory array juggling required.

# Prerequisites and Assumptions

Confirm your input file’s embedded CRS before writing the pipeline:

pdal info --metadata input.laz | python -m json.tool | grep -i srs

If the header reports EPSG:0 or empty WKT, supply in_srs explicitly rather than relying on the file — an untagged file will silently produce incorrect output coordinates. See fixing CRS mismatches in point clouds for a full repair workflow. You can also run pipeline validation on the source pipeline JSON before executing it against real data.

# Step-by-Step Implementation

# Step 1 — Identify the source EPSG code

UTM zones follow a deterministic pattern. For a file covering the US East Coast (e.g. UTM Zone 18N, WGS84 datum): EPSG:32618. For NAD83 datum equivalents, the EPSG range is EPSG:26901–26923. Run pdal info --metadata or open the file in QGIS to read the VLR CRS before hardcoding anything.

# Step 2 — Build the reprojection pipeline

The minimal JSON pipeline has three stages: a reader, the reprojection filter, and a writer with explicit CRS tagging. This follows the same PDAL stage chaining pattern used across all PDAL workflows — each stage passes its point buffer to the next in sequence.

{
  "pipeline": [
    "input_utm.laz",
    {
      "type": "filters.reprojection",
      "in_srs": "EPSG:32618",
      "out_srs": "EPSG:4326"
    },
    {
      "type": "writers.las",
      "filename": "output_wgs84.laz",
      "a_srs": "EPSG:4326",
      "forward": "all"
    }
  ]
}

Key decisions:

• in_srs overrides whatever the header says; always supply it explicitly when the source CRS might be absent or wrong.
• out_srs: "EPSG:4326" targets 2D geographic coordinates. X becomes longitude, Y becomes latitude — both in decimal degrees.
• a_srs on writers.las writes the target CRS as a WKT2 VLR into the output header. Without it, the coordinates shift but the header still claims the old UTM CRS.
• forward: "all" propagates every non-geometric VLR and dimension (intensity, return number, classification) into the output unchanged.

# Step 3 — Execute in Python

import pdal
import json

def reproject_utm_to_wgs84(
    input_path: str,
    output_path: str,
    source_epsg: int,
) -> int:
    """
    Reproject a LAS/LAZ file from UTM to WGS84 geographic coordinates.

    Args:
        input_path:   Path to source LAS/LAZ file in UTM.
        output_path:  Destination path for the WGS84 output.
        source_epsg:  EPSG integer for the source UTM CRS (e.g. 32618 for UTM Zone 18N).

    Returns:
        Number of points written to the output file.

    Raises:
        RuntimeError: If PDAL processes zero points (empty file or bad pipeline).
    """
    pipeline_dict = {
        "pipeline": [
            input_path,
            {
                "type": "filters.reprojection",
                "in_srs": f"EPSG:{source_epsg}",
                "out_srs": "EPSG:4326"
            },
            {
                "type": "writers.las",
                "filename": output_path,
                "a_srs": "EPSG:4326",
                "forward": "all"
            }
        ]
    }

    pipeline = pdal.Pipeline(json.dumps(pipeline_dict))
    count = pipeline.execute()

    if count == 0:
        raise RuntimeError(
            f"PDAL executed but processed 0 points from '{input_path}'. "
            "Check that in_srs matches the actual file CRS."
        )

    return count

PDAL’s streaming model processes the file in configurable chunks — this pipeline handles multi-gigabyte files without loading the full point cloud into RAM, unlike the array-based fallback in the complete example below.

# Step 4 — Lightweight fallback with pyproj + laspy

When PDAL cannot be installed (minimal containers, serverless functions), transform coordinates with laspy and pyproj. This loads the entire cloud into memory, so it is unsuitable for files above a few hundred MB unless your environment has the RAM budget.

import laspy
from pyproj import Transformer, CRS

def reproject_utm_to_wgs84_laspy(
    input_path: str,
    output_path: str,
    source_epsg: int,
) -> None:
    """
    Memory-resident reprojection fallback via laspy + pyproj.
    Loads the entire point cloud; use only for small files (<500 MB).
    Z values are NOT transformed — they remain in the original vertical datum.
    """
    with laspy.open(input_path) as f:
        las = f.read()

    # always_xy=True guarantees (longitude, latitude) ordering regardless of CRS axis order
    transformer = Transformer.from_crs(
        f"EPSG:{source_epsg}", "EPSG:4326", always_xy=True
    )

    lon, lat = transformer.transform(las.x, las.y)

    # laspy 2.x recalculates scale/offset on assignment via .x / .y properties
    las.x = lon
    las.y = lat
    # Z intentionally left unchanged — NAVD88/EGM96 heights remain valid

    # Inject a WKT2 VLR so GIS tools recognise the new CRS
    target_crs = CRS.from_epsg(4326)
    las.header.vlrs = [
        laspy.vlrs.VLR(
            user_id="LASF_Projection",
            record_id=2112,
            description="WKT Coordinate System",
            record_data=target_crs.to_wkt().encode("utf-8"),
        )
    ]
    las.header.global_encoding.wkt = True

    las.write(output_path)
    print(f"Transformed {len(las.x):,} points to EPSG:4326.")

# Complete Working Example

The following script can be copied, saved as reproject_utm_wgs84.py, and run against any LAS/LAZ file in UTM. It wires together the PDAL pipeline approach, verification, and coordinate range assertion in a single executable module.

#!/usr/bin/env python3
"""
reproject_utm_wgs84.py
Reproject a LAS/LAZ file from any UTM zone to WGS84 (EPSG:4326) using PDAL.

Usage:
    python reproject_utm_wgs84.py input_utm18n.laz output_wgs84.laz 32618

Requirements:
    pip install pdal
    PDAL 2.4+ with PROJ 8+ (available via conda-forge: conda install -c conda-forge pdal)
"""

import json
import sys

import pdal


def reproject_utm_to_wgs84(
    input_path: str,
    output_path: str,
    source_epsg: int,
) -> int:
    """Reproject UTM LAS/LAZ to WGS84 geographic coordinates (EPSG:4326)."""
    pipeline_dict = {
        "pipeline": [
            input_path,
            {
                "type": "filters.reprojection",
                "in_srs": f"EPSG:{source_epsg}",
                "out_srs": "EPSG:4326",
            },
            {
                "type": "writers.las",
                "filename": output_path,
                "a_srs": "EPSG:4326",
                "forward": "all",
            },
        ]
    }

    pipeline = pdal.Pipeline(json.dumps(pipeline_dict))
    count = pipeline.execute()

    if count == 0:
        raise RuntimeError(
            f"PDAL processed 0 points from '{input_path}'. "
            "Verify that in_srs matches the actual file CRS."
        )

    return count


def verify_reprojection(input_path: str, output_path: str, expected_count: int) -> None:
    """Assert point count parity and check output coordinate ranges."""

    def count_points(path: str) -> int:
        p = pdal.Pipeline(json.dumps({"pipeline": [path]}))
        p.execute()
        return p.arrays[0].shape[0]

    n_out = count_points(output_path)
    assert n_out == expected_count, (
        f"Point count mismatch: {expected_count:,} in vs {n_out:,} out. "
        "Check that no filter stage was unintentionally applied."
    )

    # Spot-check coordinate ranges — WGS84 lon is [-180, 180], lat is [-90, 90]
    p = pdal.Pipeline(json.dumps({"pipeline": [output_path]}))
    p.execute()
    pts = p.arrays[0]
    assert pts["X"].min() > -180 and pts["X"].max() < 180, "X outside WGS84 longitude range"
    assert pts["Y"].min() > -90 and pts["Y"].max() < 90, "Y outside WGS84 latitude range"

    print(f"OK: {n_out:,} points verified at EPSG:4326 coordinate ranges.")


def main() -> None:
    if len(sys.argv) != 4:
        print("Usage: python reproject_utm_wgs84.py <input.laz> <output.laz> <source_epsg>")
        sys.exit(1)

    input_path = sys.argv[1]
    output_path = sys.argv[2]
    source_epsg = int(sys.argv[3])

    print(f"Reprojecting {input_path} (EPSG:{source_epsg}) -> {output_path} (EPSG:4326)...")
    count = reproject_utm_to_wgs84(input_path, output_path, source_epsg)
    print(f"Reprojected {count:,} points.")

    verify_reprojection(input_path, output_path, count)


if __name__ == "__main__":
    main()

# Key Parameter Table

# Including a vertical datum transform

If your deliverable requires heights in WGS84 ellipsoidal metres (e.g. for web APIs or volumetric analysis tools), chain the horizontal and vertical transformation with a compound CRS:

{
  "type": "filters.reprojection",
  "in_srs": "EPSG:32618+5703",
  "out_srs": "EPSG:4979"
}

EPSG:32618+5703 combines UTM Zone 18N with NAVD88 heights. EPSG:4979 is 3D WGS84 (lon/lat/ellipsoidal height). PROJ applies the GEOID12B or GEOID18 grid automatically when it is present in $PROJ_DATA.

# Verification

After running either workflow, confirm the transformation succeeded before the file enters any downstream pipeline or is delivered to a client.

Check the output header CRS:

pdal info --metadata output_wgs84.laz | python -m json.tool | grep -i "srs\|epsg\|wkt"

Expect to see EPSG:4326 or a WKT2 string containing "GEOGCRS["WGS 84" in the metadata output.

Assert point count parity:

import pdal, json

def verify_reprojection(input_path: str, output_path: str) -> None:
    def count_points(path: str) -> int:
        p = pdal.Pipeline(json.dumps({"pipeline": [path]}))
        p.execute()
        return p.arrays[0].shape[0]

    n_in  = count_points(input_path)
    n_out = count_points(output_path)

    assert n_in == n_out, (
        f"Point count mismatch: {n_in:,} in vs {n_out:,} out. "
        "A filter stage may have been unintentionally applied."
    )
    print(f"OK: {n_out:,} points, CRS updated to EPSG:4326.")

Spot-check coordinate range:

WGS84 decimal degrees are always in the range [-180, 180] for longitude and [-90, 90] for latitude. If output X values are still in the millions, the header was read but coordinates were not actually transformed — a sign that in_srs was silently ignored because it matched an already-geographic header.

p = pdal.Pipeline(json.dumps({"pipeline": ["output_wgs84.laz"]}))
p.execute()
pts = p.arrays[0]
assert pts["X"].min() > -180 and pts["X"].max() < 180, "X outside WGS84 longitude range"
assert pts["Y"].min() >  -90 and pts["Y"].max() <  90, "Y outside WGS84 latitude range"

# Gotchas and Edge Cases

1. Header says UTM but coordinates are already geographic. Some LAZ files created by third-party exporters embed the wrong CRS in the header while the stored X/Y are already in decimal degrees. If you apply UTM to WGS84 to these, the output is garbage. Always run pdal info --stats and eyeball the X/Y ranges: UTM easting typically falls between 160,000 and 834,000 m; values below 180 almost certainly indicate geographic coordinates already.

2. Missing PROJ datum grid files produce silent metre-scale errors. When PROJ cannot find a required transformation grid (e.g. us_noaa_conus.tif for NAD27 to WGS84), it silently falls back to an approximate 7-parameter Helmert transformation, introducing errors of 1–10 m. Install the full PROJ grid package before running survey-grade reprojection:

conda install -c conda-forge proj-data
# or
pip install pyproj[grids]

Verify available transformations with:

projinfo -s EPSG:32618 -t EPSG:4326 --summary

3. Scale/offset corruption in the laspy fallback. The laspy fallback assigns decimal-degree values to las.x and las.y, which triggers automatic scale/offset recalculation. However, if the original LAS file had a large UTM northing offset (e.g. 4,500,000 m), the auto-offset may store coordinates with insufficient decimal places. Check las.header.offsets and las.header.scales after assignment and confirm the stored precision is no larger than 1e-7 degrees (about 1 cm at the equator).

4. Mixing horizontal WGS84 degrees with vertical metre values. A common oversight in survey workflows: horizontal units shift to decimal degrees while Z remains in metres. Downstream tools that compute 3D distances will produce wildly incorrect results if they assume a single unit system. Document this in the output file’s VLR description, or apply the compound CRS approach in the parameter table above to convert Z to ellipsoidal heights in the same pass.

# Frequently Asked Questions

Does filters.reprojection require the input file to have a valid CRS in its header?

No. If you supply in_srs explicitly in the pipeline JSON, PDAL uses that value regardless of what the header says. This is the correct approach for files with empty or incorrect CRS metadata — a situation that fixing CRS mismatches in point clouds covers in depth.

Does PDAL automatically update the LAS header CRS after reprojection?

Yes — when you set a_srs on writers.las, PDAL writes a WKT2 VLR into the output header. Without a_srs the header retains the source CRS, which will confuse downstream GIS tools. Parsing and checking that VLR is covered in how to parse LAS headers with Python.

Are Z values transformed when reprojecting from UTM to WGS84?

Only if you include a compound CRS (e.g. EPSG:32618+5703) in out_srs. A plain EPSG:4326 target leaves Z values in their original vertical datum (NAVD88, EGM96, etc.) untouched, as shown in the diagram above.

Can I chain reprojection with other filters in the same pipeline?

Yes, and that is the normal pattern. You can combine reprojection with statistical outlier removal in a single pipeline pass — see applying statistical outlier filters in PDAL. Place filters.reprojection before any filter that operates in geographic space (e.g. spatial bounds filtering) but after filters that work on raw intensity or classification dimensions.

What if I need to reproject to a different UTM zone rather than WGS84?

Change out_srs to the target UTM EPSG code (e.g. "EPSG:32619" for UTM Zone 19N). The same pipeline structure applies; set a_srs on the writer to match the output zone.

# Related

• Spatial Reprojection in PDAL — parent guide covering CRS concepts, reprojection stage parameters, and batch reprojection patterns
• PDAL Pipeline Architecture & Execution — how readers, filters, and writers chain together into a streaming execution model
• Fixing CRS Mismatches in Point Clouds — diagnose and repair embedded CRS errors before reprojection
• Applying Statistical Outlier Filters in PDAL — combine with reprojection in a single pipeline pass
• LAS/LAZ File Structure — VLR layout and header fields that carry the CRS after transformation