Fixing SRS Lost When Writing LAS Files

The client reports that the point cloud opens “in the wrong place” — sometimes at the origin, sometimes overlaid on the right site but a few hundred metres out, sometimes fine. On the machine that produced it everything looks correct, which is the detail that makes this frustrating: the processing environment has enough context to place the cloud even when the file does not carry it.

The cause is almost always that the spatial reference was never written into the file, or was written in a form the reading software does not look for. Neither produces an error at write time, and the producing machine’s own tools often fall back to a project setting that masks the problem entirely.

# Where a LAS file keeps its coordinate system

LAS has carried spatial reference information three different ways over its versions, and readers vary in which they consult.

GeoTIFF keys, in variable-length records, are the LAS 1.0–1.3 mechanism. They encode a projected coordinate system code and a handful of parameters, and they cannot express everything modern CRS definitions contain — notably compound horizontal-plus-vertical systems and time-dependent transformations.

OGC WKT, in a variable-length record, is the LAS 1.4 mechanism. It is a full text definition and can express compound systems. Crucially, LAS 1.4 has a global encoding bit that declares whether WKT is present, and a file with WKT but without that bit set is one many readers will not look in.

Nothing at all, which is what a writer produces when it was given no CRS and asked no questions.

A fourth mechanism exists outside the file: a .prj sidecar. It is not part of the LAS specification, several tools read it, and it costs nothing to write.

Figure 1 — Why “the CRS is in the file” is not the same as “readers will find it”.

# Minimal reproducible solution

import laspy
import numpy as np
from pyproj import CRS


def write_with_crs(xyz: np.ndarray, out_path: str, epsg: int,
                   *, vertical_epsg: int | None = None,
                   write_sidecar: bool = True) -> dict:
    """Write LAS 1.4 with an unambiguous, discoverable coordinate system.

    add_crs sets both the WKT record and the global encoding bit that tells a
    reader to look for it. Setting one without the other is the single most
    common way a correctly intentioned write produces an unreadable CRS.
    """
    crs = CRS.from_epsg(epsg)
    if vertical_epsg is not None:
        crs = CRS.from_user_input(
            f"EPSG:{epsg}+{vertical_epsg}")     # compound horizontal + vertical

    header = laspy.LasHeader(point_format=6, version="1.4")
    header.offsets = np.round(xyz.mean(axis=0))
    header.scales = np.array([0.001, 0.001, 0.001])
    header.add_crs(crs)

    las = laspy.LasData(header)
    las.x, las.y, las.z = xyz[:, 0], xyz[:, 1], xyz[:, 2]
    las.write(out_path)

    if write_sidecar:
        from pathlib import Path
        Path(out_path).with_suffix(".prj").write_text(crs.to_wkt())

    return {"epsg": epsg, "compound": vertical_epsg is not None,
            "sidecar": write_sidecar}

The compound-CRS branch matters more than it looks. A cloud whose heights are orthometric but whose declared CRS is horizontal-only carries no statement about the vertical datum at all, which is how a survey ends up differenced against another one on a different datum. The vertical side is covered in geoid models and vertical datum automation.

For PDAL-written files, the equivalent is to forward the source CRS explicitly, and to set it when the source has none:

import json
import subprocess


def write_with_pdal(src: str, dst: str, epsg: int,
                    *, assume_source: str | None = None) -> None:
    """Write through PDAL with the CRS explicitly assigned, not inherited.

    `a_srs` on the reader assigns a CRS to data that has none — it does not
    reproject. Confusing it with `out_srs`, which does, silently produces
    coordinates in the wrong system with a correct-looking declaration.
    """
    reader: dict = {"filename": src, "type": "readers.las"}
    if assume_source:
        reader["spatialreference"] = assume_source

    pipeline = {"pipeline": [
        reader,
        {"type": "writers.las", "filename": dst, "compression": "laszip",
         "a_srs": f"EPSG:{epsg}", "forward": "all", "minor_version": 4},
    ]}
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)

Figure 3 — Three losses, one remedy.

# Edge-case matrix

The a_srs versus out_srs distinction deserves the emphasis. a_srs asserts — it relabels data without moving it. out_srs transforms. Using the first where the second was intended produces a file whose coordinates are in the old system and whose header claims the new one, which is worse than having no CRS at all because it looks authoritative.

# Verification snippet

from pathlib import Path

import laspy
from pyproj import CRS


def verify_crs(path: str, expect_epsg: int) -> dict:
    """Read back what a reader would actually find, and compare."""
    las = laspy.read(path)
    found = las.header.parse_crs()

    problems = []
    if found is None:
        problems.append("no CRS discoverable in the file")
    else:
        got = CRS.from_user_input(found)
        if got.to_epsg() != expect_epsg:
            problems.append(f"declares EPSG:{got.to_epsg()}, expected {expect_epsg}")
        if not got.is_projected:
            problems.append("CRS is geographic — metre-based scales are wrong for it")
        if len(got.sub_crs_list) < 2 and "orthometric" in (got.name or "").lower():
            problems.append("heights look orthometric but no vertical CRS is declared")

    sidecar = Path(path).with_suffix(".prj")
    return {"epsg": (CRS.from_user_input(found).to_epsg() if found else None),
            "version": str(las.header.version),
            "sidecar_present": sidecar.exists(),
            "problems": problems}

The geographic check is worth including because it catches a compound mistake: a file whose CRS is in degrees but whose header scale of 0.001 was chosen for metres. That combination quantises to roughly 100 m and is the single most extreme precision loss available in the format.

Figure 2 — Two options, one letter apart in intent, with opposite effects on the numbers.

# Making the declaration survive the whole pipeline

A CRS is lost somewhere between the reconstruction and the delivery far more often than it is written wrongly at the end. Each intermediate write is an opportunity to drop it, and a pipeline of six stages has six of them.

Two conventions close the gap. Forward the header on every intermediate write — forward: all in PDAL, or explicit add_crs in laspy — so a CRS set once at ingest propagates to the end without being restated. And assert the CRS after every write, not only at delivery, so the stage that lost it is identified immediately rather than at the end of a long chain.

def crs_trace(paths: list[str], expect_epsg: int) -> list[dict]:
    """Walk the intermediate files of a run and report where the CRS was lost."""
    import laspy
    from pyproj import CRS

    trace = []
    for path in paths:
        try:
            with laspy.open(path) as fh:
                found = fh.header.parse_crs()
            epsg = CRS.from_user_input(found).to_epsg() if found else None
        except Exception as exc:
            epsg, found = None, str(exc)[:80]
        trace.append({"file": path, "epsg": epsg,
                      "ok": epsg == expect_epsg})
        if epsg != expect_epsg:
            break                      # the first loss is the one to fix
    return trace

Stopping at the first loss is deliberate: once a CRS is gone, every file after it inherits the absence, and reporting all of them buries the one that matters. The same pattern applies to any header property worth preserving — scale, offset, extra dimensions — and a single trace function covering all of them is a worthwhile twenty lines on a pipeline that runs weekly.

# When to escalate

• The data is in a local engineering grid with no EPSG code. There is nothing correct to write. Document the transformation to a national grid alongside the delivery and do not assign a code that approximates it.
• A client’s software rewrites or discards the CRS on import. Some packages normalise on load. Ask what their workflow preserves before choosing a delivery format, and supply the sidecar as well.
• Two files of the same site declare different systems. Find out which is right before reprojecting either; a wrong declaration propagated through a reprojection is much harder to unpick.

# Related

• Troubleshooting Point Cloud Processing Failures
• Point Cloud Formats and Interchange in Python
• Managing Coordinate Reference Systems in GDAL

← Troubleshooting Point Cloud Processing Failures