Reading the CRS from LAS WKT and GeoTIFF Keys

TL;DR: laspy.open(path).header.parse_crs() returns a pyproj CRS from either the WKT record (LAS 1.4, user ID LASF_Projection, record ID 2112) or the GeoTIFF key records (34735 key directory, 34736 doubles, 34737 ASCII). pdal info --metadata reports the same as metadata.srs with horizontal, vertical and compoundwkt fields. Always check that a vertical component exists and that the declared units match the coordinate magnitudes.

# Context and Motivation

This guide is part of Coordinate Reference Systems. A LAS file does not have to declare its CRS, and when it does, it can do so in two different ways. LAS 1.4 stores an OGC Well-Known Text string in a variable-length record; older files use the GeoTIFF key mechanism borrowed from raster formats, with EPSG codes and parameters spread across three records. Both can be incomplete — a horizontal CRS without the vertical datum is common — and a file can even contain both, disagreeing. Reading the CRS correctly is the first step of every reprojection, merge and accuracy assessment, and it takes only a few lines.

# Prerequisites and Assumptions

• laspy 2.x with pyproj, or PDAL 2.x.
• LAS or LAZ files of any version.
• A rough idea of where the data is on Earth, to sanity-check coordinates against the declared CRS.

# Step-by-Step Implementation

# Step 1 — List the VLRs

laspy.open(p).header.vlrs lists every record with user ID and record ID. Look for LASF_Projection records: 2112 (WKT), 34735/34736/34737 (GeoTIFF keys).

# Step 2 — Parse the CRS

header.parse_crs() returns a pyproj CRS or None. In PDAL, pdal info --metadata gives metadata.srs.wkt, horizontal, vertical and compoundwkt.

# Step 3 — Split horizontal and vertical

A compound CRS has sub_crs_list with a projected and a vertical part. A plain projected CRS has no vertical information.

# Step 4 — Check units and magnitudes

Compare the declared axis units with the coordinate magnitudes: UTM eastings in metres range roughly 160,000–840,000; State Plane eastings in feet are often above 1,000,000.

# Step 5 — Resolve conflicts

If WKT and GeoTIFF keys both exist and disagree, LAS 1.4 readers follow the WKT when global encoding bit 4 is set. Decide which is right from the delivery documentation, then rewrite the file with a single correct CRS.

# Complete Working Example

"""Report what CRS each LAS/LAZ file declares, and flag common problems."""
from __future__ import annotations

from pathlib import Path

import laspy
import numpy as np
from pyproj import CRS

WKT_ID, GEOKEY_IDS = 2112, {34735, 34736, 34737}


def crs_report(path: Path) -> dict:
    with laspy.open(path) as f:
        h = f.header
        proj_vlrs = {v.record_id for v in h.vlrs if v.user_id == "LASF_Projection"}
        proj_vlrs |= {v.record_id for v in getattr(h, "evlrs", []) or [] if v.user_id == "LASF_Projection"}
        crs = h.parse_crs()
        wkt_bit = bool(int(h.global_encoding.wkt)) if h.version.minor >= 4 else False
        mins, maxs = h.mins, h.maxs

    rep = {"file": path.name, "las": f"1.{h.version.minor}",
           "has_wkt": WKT_ID in proj_vlrs, "has_geokeys": bool(proj_vlrs & GEOKEY_IDS),
           "wkt_bit": wkt_bit, "crs": None, "horizontal": None, "vertical": None,
           "unit": None, "issues": []}
    if crs is None:
        rep["issues"].append("no CRS declared")
        return rep
    rep["crs"] = crs.name
    subs = crs.sub_crs_list or [crs]
    horiz = next((c for c in subs if c.is_projected or c.is_geographic), None)
    vert = next((c for c in subs if c.is_vertical), None)
    rep["horizontal"] = horiz.to_epsg() if horiz else None
    rep["vertical"] = vert.name if vert else None
    if horiz is not None and horiz.axis_info:
        rep["unit"] = horiz.axis_info[0].unit_name
    if vert is None:
        rep["issues"].append("no vertical CRS")
    if rep["has_wkt"] and rep["has_geokeys"]:
        rep["issues"].append("both WKT and GeoTIFF keys present")
    if h.version.minor >= 4 and rep["has_wkt"] and not wkt_bit:
        rep["issues"].append("WKT present but global encoding bit 4 not set")
    if horiz is not None and horiz.is_projected:
        east = float(np.mean([mins[0], maxs[0]]))
        if rep["unit"] == "metre" and east > 2_000_000:
            rep["issues"].append("eastings look like feet for a metre CRS")
    return rep


if __name__ == "__main__":
    for p in sorted(Path("received").glob("*.la[sz]")):
        r = crs_report(p)
        print(f"{r['file']:<24} LAS {r['las']} EPSG {r['horizontal']} / {r['vertical']} "
              f"({r['unit']})  {'; '.join(r['issues']) or 'ok'}")

From the command line, the PDAL equivalent for one file:

pdal info --metadata received/t_0431.laz | jq '.metadata.srs | {horizontal, vertical, units}'

# Key Parameter Table

# Verification

• Magnitudes match units. Coordinates in the header bounds should be plausible for the declared CRS and unit.
• Location check. Transform the header’s centre point to longitude and latitude with pyproj; it should land where the project is.
• Consistency across tiles. Every tile in a delivery should declare the same CRS; group the report by crs and investigate minorities.

from pyproj import Transformer
with laspy.open("received/t_0431.laz") as f:
    h = f.header
    t = Transformer.from_crs(h.parse_crs(), "EPSG:4326", always_xy=True)
    print(t.transform((h.mins[0] + h.maxs[0]) / 2, (h.mins[1] + h.maxs[1]) / 2))

# Gotchas and Edge Cases

Horizontal only. The single most common gap. Heights exist but their datum is undocumented; see setting a vertical CRS on a point cloud.

User-defined GeoTIFF projections. Older files sometimes encode a projection by parameters rather than an EPSG code. pyproj parses them into a custom CRS without an EPSG number; compare parameters carefully with the expected zone.

CRS in extended VLRs. LAS 1.4 permits the WKT in an extended VLR at the end of the file. Readers that only scan the standard VLR block miss it; laspy and PDAL handle both.

Trusting the CRS over the data. A header can simply be wrong. When magnitudes contradict the declared CRS, believe the numbers and fix the header — fixing CRS mismatches in point clouds covers the repair.

# Frequently Asked Questions

How do I find the coordinate system of a LAS file?

Open it with laspy and call header.parse_crs, which returns a pyproj CRS, or run pdal info with the metadata flag and read the srs section. Both read the WKT or GeoTIFF key records.

What is VLR record 2112?

It is the LAS 1.4 variable-length record, under user ID LASF_Projection, that holds the coordinate reference system as an OGC Well-Known Text string. LAS 1.4 point formats 6 to 10 require the CRS in this form.

Why does my file have no vertical datum?

Many files declare only the horizontal CRS, often because the source software wrote an EPSG code for the projection alone. The heights are still in some vertical datum; find it from the delivery documentation and set it explicitly.

What if WKT and GeoTIFF keys disagree?

LAS 1.4 readers use the WKT when global encoding bit 4 is set. Determine which description is correct from the project documentation, then rewrite the file with a single, correct CRS.

# Related

• Coordinate Reference Systems — CRS handling overview
• Choosing a Projected CRS for a LiDAR Project — picking a target
• Setting a Vertical CRS on a Point Cloud — filling the vertical gap
• Reading and Writing LAS VLRs with PDAL — the VLR mechanism
• Reprojecting State Plane Feet to Metres — when units are feet