Fixing CRS Mismatches in Point Clouds with Python

import laspy
from pyproj import CRS

def diagnose_crs(las_path: str) -> dict:
    """
    Inspect a LAS/LAZ file and return a dict describing its CRS state.
    Returns: {'has_wkt2': bool, 'has_geokey': bool, 'crs': CRS | None, 'note': str}
    """
    result = {"has_wkt2": False, "has_geokey": False, "crs": None, "note": ""}
    with laspy.open(las_path) as f:
        header = f.header
        for vlr in header.vlrs:
            if vlr.record_id == 2112:          # WKT2 record
                result["has_wkt2"] = True
                try:
                    wkt = vlr.record_data.decode("utf-8").rstrip("\x00")
                    result["crs"] = CRS.from_wkt(wkt)
                    result["note"] = "WKT2 CRS parsed from VLR 2112"
                except Exception as e:
                    result["note"] = f"Malformed WKT2 VLR: {e}"
            elif vlr.record_id == 34735:        # GeoKey directory (legacy)
                result["has_geokey"] = True
                if result["crs"] is None:
                    result["note"] = "Only legacy GeoKey CRS found — WKT2 preferred"

    if result["crs"] is None and not result["has_geokey"]:
        result["note"] = "No CRS in VLRs — must be supplied manually"
    return result

from pyproj import Transformer

def build_transformer(source_crs: CRS, target_epsg: int) -> Transformer:
    """
    Build a pyproj Transformer that preserves axis order and applies
    vertical grid shifts when the target CRS includes a geoid model.
    """
    target_crs = CRS.from_epsg(target_epsg)
    return Transformer.from_crs(
        source_crs,
        target_crs,
        always_xy=True,   # Treat first axis as easting, second as northing
    )

import numpy as np

def transform_xyz(las: laspy.LasData, transformer: Transformer) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Apply transformer to X/Y/Z arrays. laspy 2.x applies scale+offset
    automatically, so las.x / las.y / las.z are already in CRS units.
    """
    new_x, new_y, new_z = transformer.transform(las.x, las.y, las.z)
    return new_x, new_y, new_z

import laspy

def write_corrected_las(
    source_las: laspy.LasData,
    new_x: np.ndarray,
    new_y: np.ndarray,
    new_z: np.ndarray,
    target_crs: CRS,
    output_path: str,
) -> None:
    """
    Write a new LAS file with transformed coordinates and the target CRS
    embedded as a WKT2 VLR. All non-coordinate dimensions are preserved verbatim.
    """
    new_header = laspy.LasHeader(
        point_format=source_las.header.point_format.id,
        version=source_las.header.version,
    )
    new_header.offsets = source_las.header.offsets
    new_header.scales = source_las.header.scales

    wkt_bytes = target_crs.to_wkt().encode("utf-8")
    new_header.vlrs = [
        laspy.vlrs.VLR(
            user_id="LASF_Projection",
            record_id=2112,
            description="WKT Coordinate System",
            record_data=wkt_bytes,
        )
    ]
    new_header.global_encoding.wkt = True   # Required for LAS 1.4 compliance

    new_las = laspy.LasData(header=new_header)
    new_las.x = new_x
    new_las.y = new_y
    new_las.z = new_z

    # Copy every non-coordinate dimension: intensity, classification,
    # return_number, scan_angle_rank, rgb, waveform fields, etc.
    for dim_name in source_las.point_format.dimension_names:
        if dim_name.lower() not in ("x", "y", "z"):
            setattr(new_las, dim_name, getattr(source_las, dim_name))

    new_las.write(output_path)

import laspy
import numpy as np
from pyproj import CRS, Transformer
from pathlib import Path

# ── Configure these three values ──────────────────────────────────────────────
INPUT_PATH   = "survey_raw.laz"
OUTPUT_PATH  = "survey_utm18n_navd88.laz"
TARGET_EPSG  = 6347           # NAD83(2011) / UTM zone 18N (horizontal only)
# For a compound CRS add the vertical: CRS.from_authority("EPSG", "6347+5703")
# ──────────────────────────────────────────────────────────────────────────────


def extract_crs(header: laspy.LasHeader) -> CRS | None:
    for vlr in header.vlrs:
        if vlr.record_id == 2112:
            try:
                return CRS.from_wkt(vlr.record_data.decode("utf-8").rstrip("\x00"))
            except Exception:
                pass
    return None


def fix_crs_mismatch(input_path: str, output_path: str, target_epsg: int) -> None:
    input_p  = Path(input_path)
    output_p = Path(output_path)

    if not input_p.exists():
        raise FileNotFoundError(f"Input file not found: {input_path}")

    with laspy.open(input_p) as f:
        las = f.read()

    source_crs = extract_crs(las.header)
    if source_crs is None:
        raise ValueError(
            "No CRS detected in LAS header VLRs. "
            "Set source_crs = CRS.from_epsg(<code>) manually and retry."
        )

    target_crs  = CRS.from_epsg(target_epsg)
    transformer = Transformer.from_crs(source_crs, target_crs, always_xy=True)

    new_x, new_y, new_z = transformer.transform(las.x, las.y, las.z)

    new_header = laspy.LasHeader(
        point_format=las.header.point_format.id,
        version=las.header.version,
    )
    new_header.offsets = las.header.offsets
    new_header.scales  = las.header.scales
    new_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"),
        )
    ]
    new_header.global_encoding.wkt = True

    new_las   = laspy.LasData(header=new_header)
    new_las.x = new_x
    new_las.y = new_y
    new_las.z = new_z

    for dim_name in las.point_format.dimension_names:
        if dim_name.lower() not in ("x", "y", "z"):
            setattr(new_las, dim_name, getattr(las, dim_name))

    new_las.write(output_p)
    print(
        f"Transformed {len(new_x):,} points "
        f"from {source_crs.to_epsg()} → EPSG:{target_epsg}"
    )
    print(f"Output: {output_p}")


if __name__ == "__main__":
    fix_crs_mismatch(INPUT_PATH, OUTPUT_PATH, TARGET_EPSG)

import laspy
import numpy as np
from pyproj import CRS

def verify_corrected_file(output_path: str, expected_epsg: int) -> None:
    with laspy.open(output_path) as f:
        las = f.read()

    # 1. Confirm CRS round-trips correctly
    found_crs = None
    for vlr in las.header.vlrs:
        if vlr.record_id == 2112:
            found_crs = CRS.from_wkt(vlr.record_data.decode("utf-8").rstrip("\x00"))
    assert found_crs is not None, "No WKT2 VLR found in output"
    assert found_crs.to_epsg() == expected_epsg, (
        f"CRS mismatch: expected EPSG:{expected_epsg}, got {found_crs.to_epsg()}"
    )

    # 2. Bounding box sanity — for UTM, X should be 100 000–900 000, Y > 0
    assert 100_000 < np.min(las.x) < 900_000, f"X out of UTM range: {np.min(las.x)}"
    assert np.min(las.y) > 0,                  f"Negative northing: {np.min(las.y)}"

    # 3. Point count preserved
    print(f"Point count: {len(las.x):,} — OK")
    print(f"CRS verified: EPSG:{found_crs.to_epsg()}")
    print(f"Bounds X [{np.min(las.x):.1f}, {np.max(las.x):.1f}]")
    print(f"Bounds Y [{np.min(las.y):.1f}, {np.max(las.y):.1f}]")
    print(f"Bounds Z [{np.min(las.z):.3f}, {np.max(las.z):.3f}]")