Python LiDAR

Syncing Metadata Between LAS and Shapefiles: A Python Workflow

TL;DR: Extract the bounding box and CRS from a LAS header with laspy, build a GeoDataFrame bounding polygon with geopandas, enforce DBF field-name (10-char) and string-value (254-char) limits, then call to_file() — geopandas writes the synchronized .prj, .dbf, .shp, and .shx in one step.

# Context and Motivation

This guide is part of Metadata & Header Sync, which covers the full lifecycle of validating and reconciling LAS/LAZ header fields with point data and external spatial files.

Point clouds and vector boundaries rarely share identical metadata pipelines. LAS files embed spatial context inside a structured binary header: generation date, software ID, bounding box extents, and a coordinate reference system stored in variable-length records (VLRs). Shapefiles split equivalent metadata across three sidecars — .prj (CRS as WKT), .dbf (tabular attributes), and optional .xml (extended metadata) — each with its own format constraints and legacy size limits.

When engineering teams generate delivery boundaries by deriving a bounding polygon from a LiDAR survey tile, metadata typically degrades at the format handoff. CRS definitions get dropped, software ID strings exceed the 254-character DBF limit and are silently truncated, and generation timestamps disappear entirely. The result is an audit gap: the shapefile boundary can no longer be traced back unambiguously to its source scan. For infrastructure and urban-planning deliverables subject to survey-grade QA, that gap is a liability.

Understanding how a LAS header is structured is covered in LAS/LAZ File Structure. Resolving CRS discrepancies before synchronization is covered in Coordinate Reference Systems.

LAS to Shapefile Metadata Synchronization Flow Five sequential stages: (1) LAS header parsed for VLRs and extents, (2) CRS resolved via WKT2 VLR or fallback EPSG, (3) bounding box polygon built in GeoDataFrame, (4) DBF attributes mapped with 10-char name and 254-char value limits enforced, (5) shapefile written with .prj, .dbf, .shp sidecars. LAS Header VLRs + extents CRS Resolve WKT2 / EPSG fallback Bounding Box shapely box + GDF DBF Mapping 10-char / 254-char limits Shapefile Export .prj / .dbf / .shp .shx written enforce DBF limits before export — GDAL truncates silently at stage 4

# Prerequisites and Assumptions

  • Python 3.10+
  • laspy >= 2.4.0 — for the modern laspy.open() API and VLR access
  • geopandas >= 0.14 — for GeoDataFrame.to_file() with the ESRI Shapefile driver
  • pyproj >= 3.3 — for CRS parsing and optional reprojection
  • shapely >= 2.0 — for box() geometry construction
  • Input LAS/LAZ file with valid x_min, x_max, y_min, y_max header fields
  • Target CRS known in advance (an EPSG code or WKT string) if the source file may lack an embedded CRS

Install with:

bash
pip install "laspy[lazrs]>=2.4.0" geopandas pyproj shapely

# Step-by-Step Implementation

# Step 1 — Open the LAS File and Read the Header

Use laspy.open() in read mode. This loads only the header and VLRs; the full point array stays on disk until explicitly requested, which is efficient for files larger than a few hundred megabytes.

python
import laspy

with laspy.open("survey_tile.las") as las:
    header = las.header
    print(header.version.major, header.version.minor)  # e.g., 1 4
    print(header.point_count)                          # e.g., 4_832_011
    print(header.x_min, header.x_max)                 # UTM eastings

The header’s VLRs are accessible as header.vlrs — a list of laspy.VLR objects. Each VLR carries a user_id, record_id, and raw record_data bytes. The CRS is stored in VLR record_id 2112 (WKT2, preferred) or record_id 34735 (legacy GeoKey, LAS 1.2).

# Step 2 — Extract the CRS from the VLR

python
from pyproj import CRS
import warnings

def extract_las_crs(header: laspy.LasHeader) -> CRS | None:
    """
    Return a pyproj.CRS from the WKT2 VLR (record_id 2112).
    Returns None if no recognized CRS VLR is found.
    """
    for vlr in header.vlrs:
        if vlr.record_id == 2112:
            try:
                wkt_str = vlr.record_data.decode("utf-8").rstrip("\x00")
                return CRS.from_wkt(wkt_str)
            except Exception as exc:
                warnings.warn(f"Malformed WKT2 VLR: {exc}")
    return None

LAS 1.2 files store CRS only in GeoKey VLRs (record IDs 34735–34737). For those, parsing requires a GeoTIFF key decoder — or simply accepting the fallback EPSG code described in step 3. See how to parse LAS headers with Python for a full GeoKey walkthrough.

# Step 3 — Resolve the CRS with a Safe Fallback

python
import pyproj

crs = extract_las_crs(header)
if crs is None:
    warnings.warn(
        "No embedded CRS detected in LAS VLRs. "
        "Defaulting to EPSG:4326 (WGS84). "
        "Verify this matches your project datum before delivery."
    )
    crs = pyproj.CRS.from_epsg(4326)

Always emit a warning on fallback rather than silently assuming WGS84. A silent default is indistinguishable from a correctly embedded CRS in logs.

# Step 4 — Build the Bounding Polygon

python
import geopandas as gpd
from shapely.geometry import box

with laspy.open("survey_tile.las") as las:
    header = las.header
    crs = extract_las_crs(header) or pyproj.CRS.from_epsg(4326)
    min_x, min_y = header.x_min, header.y_min
    max_x, max_y = header.x_max, header.y_max
    point_count  = header.point_count
    version_str  = f"{header.version.major}.{header.version.minor}"
    system_id    = str(header.system_identifier).strip()
    gen_soft     = str(header.generating_software).strip()

bbox = box(min_x, min_y, max_x, max_y)
gdf  = gpd.GeoDataFrame(geometry=[bbox], crs=crs)

Passing crs to the constructor is what causes geopandas to generate the .prj sidecar on export. Omitting it produces a shapefile with no spatial reference — a common silent mistake.

# Step 5 — Map Header Fields to DBF-Safe Attributes

The dBase III+ format limits field names to 10 ASCII characters and string values to 254 characters. geopandas does not validate these constraints before writing; violations silently corrupt the .dbf or raise a GDAL error at export time.

python
from datetime import datetime

gdf["LAS_VER"]   = version_str           # e.g., "1.4"
gdf["PT_COUNT"]  = int(point_count)      # numeric — no char limit
gdf["GEN_DATE"]  = datetime.now().strftime("%Y-%m-%d")
gdf["SYSTEM_ID"] = system_id[:254]       # clamp to 254 chars
gdf["GEN_SOFT"]  = gen_soft[:254]
gdf["X_MIN"]     = float(min_x)
gdf["Y_MIN"]     = float(min_y)
gdf["X_MAX"]     = float(max_x)
gdf["Y_MAX"]     = float(max_y)

All column names above are 10 characters or fewer. Use explicit float() casting for coordinate values to avoid numpy dtype coercion warnings from GDAL.

# Step 6 — Export the Shapefile

python
gdf.to_file("survey_tile_boundary.shp", driver="ESRI Shapefile")

geopandas generates four sidecars: .shp (geometry), .shx (index), .dbf (attributes), and .prj (CRS as WKT). The .prj content is derived directly from the crs you passed in step 4.

# Complete Working Example

python
import laspy
import geopandas as gpd
import pyproj
from pyproj import CRS
from shapely.geometry import box
from datetime import datetime
import warnings


def extract_las_crs(header: laspy.LasHeader) -> CRS | None:
    """
    Return a pyproj.CRS from the WKT2 VLR (record_id 2112).
    Returns None if no recognized CRS VLR is present.
    """
    for vlr in header.vlrs:
        if vlr.record_id == 2112:
            try:
                wkt_str = vlr.record_data.decode("utf-8").rstrip("\x00")
                return CRS.from_wkt(wkt_str)
            except Exception as exc:
                warnings.warn(f"Malformed WKT2 VLR — skipping: {exc}")
    return None


def sync_las_to_shapefile(las_path: str, out_shp: str) -> str:
    """
    Extract LAS header metadata, build a bounding polygon,
    and export a synchronized shapefile.

    Args:
        las_path: Path to the source LAS or LAZ file.
        out_shp:  Output .shp path (sidecars written alongside).

    Returns:
        Absolute path to the written .shp file.
    """
    # 1. Read header only — no point array loaded
    with laspy.open(las_path) as las:
        header    = las.header
        crs       = extract_las_crs(header)
        min_x, min_y = header.x_min, header.y_min
        max_x, max_y = header.x_max, header.y_max
        point_count  = header.point_count
        version_str  = f"{header.version.major}.{header.version.minor}"
        system_id    = str(header.system_identifier).strip()
        gen_soft     = str(header.generating_software).strip()

    # 2. Resolve CRS with auditable fallback
    if crs is None:
        warnings.warn(
            f"{las_path}: No embedded CRS. Defaulting to EPSG:4326 (WGS84)."
        )
        crs = pyproj.CRS.from_epsg(4326)

    # 3. Build bounding polygon and GeoDataFrame
    bbox = box(min_x, min_y, max_x, max_y)
    gdf  = gpd.GeoDataFrame(geometry=[bbox], crs=crs)

    # 4. Map to DBF-safe attributes
    gdf["LAS_VER"]   = version_str
    gdf["PT_COUNT"]  = int(point_count)
    gdf["GEN_DATE"]  = datetime.now().strftime("%Y-%m-%d")
    gdf["SYSTEM_ID"] = system_id[:254]
    gdf["GEN_SOFT"]  = gen_soft[:254]
    gdf["X_MIN"]     = float(min_x)
    gdf["Y_MIN"]     = float(min_y)
    gdf["X_MAX"]     = float(max_x)
    gdf["Y_MAX"]     = float(max_y)

    # 5. Export — generates .shp, .shx, .dbf, .prj
    gdf.to_file(out_shp, driver="ESRI Shapefile")
    return out_shp

# Key Parameter Table

Parameter / Field Type Default / Example Notes
record_id (VLR) int 2112 WKT2 CRS record; LAS 1.2 uses 34735 (GeoKey)
fallback_epsg int 4326 Override with project datum (e.g., 26918 for NAD83 UTM 18N)
DBF field name length int max 10 GDAL truncates silently; always verify manually
DBF string value length int max 254 Apply [:254] slice before assignment
driver (to_file) str "ESRI Shapefile" Use "GPKG" (GeoPackage) to lift DBF constraints for modern workflows
float() cast type numpy.float64 Prevents GDAL dtype coercion warnings for coordinate columns

# Verification

After export, confirm three properties:

1. CRS fidelity — re-read the .prj and compare against the source:

python
import pyproj

written_crs = pyproj.CRS.from_file("survey_tile_boundary.prj")
assert written_crs.equals(crs), f"CRS mismatch: {written_crs} vs {crs}"

2. Bounding box accuracy — load the shapefile and check geometry bounds:

python
result = gpd.read_file("survey_tile_boundary.shp")
bounds = result.geometry.iloc[0].bounds  # (minx, miny, maxx, maxy)
assert abs(bounds[0] - min_x) < 1e-6, "X_MIN drift detected"
assert abs(bounds[2] - max_x) < 1e-6, "X_MAX drift detected"

3. DBF field lengths — inspect with ogrinfo:

bash
ogrinfo -al -so survey_tile_boundary.shp

Look for any field name that has been auto-truncated (shorter than the original column name) and any String(n) where n < len(original_value).

# Gotchas and Edge Cases

GeoPackage over Shapefile for modern deliverables. The 10-character field-name limit and 254-character string limit are dBase III+ artefacts. If your downstream GIS accepts GeoPackage (.gpkg), pass driver="GPKG" to to_file() and all constraints disappear. Shapefiles remain necessary only for workflows locked to legacy ESRI interchange formats.

LAS 1.2 files without a WKT2 VLR. extract_las_crs() above targets record_id 2112 only. Files generated by older scanners or flight management software may embed CRS solely in GeoKey VLRs (record_id 34735). If the function returns None and you know the datum, pass the project EPSG code explicitly rather than defaulting to WGS84. Incorrect CRS assumptions compound across every downstream spatial join and can shift geometry by tens of metres when UTM zones are confused.

Bounding box header vs. actual point extents. The LAS header x_min / x_max values are written by the generating software and may be stale after in-place edits or partial point deletions. For audit-critical deliverables, verify the header extents against actual point-array statistics — the full procedure is in Metadata & Header Sync.

Empty system_identifier or generating_software fields. Some encoders leave these fields as null bytes. str(header.system_identifier).strip() collapses a null-filled 32-byte field to an empty string, which is a valid DBF string value. Log the empty value rather than assigning a placeholder to avoid fabricating metadata.

Reprojection before export. If the LAS CRS does not match the project delivery CRS, reproject before writing:

python
target_crs = pyproj.CRS.from_epsg(26918)  # NAD83 / UTM zone 18N
if not gdf.crs.equals(target_crs):
    gdf = gdf.to_crs(target_crs)

Always reproject the GeoDataFrame, not just the CRS attribute. Changing only gdf.crs reassigns the label without transforming coordinates and is a silent correctness error.

Related in Metadata & Header Sync in Python LiDAR Workflows

↑ Metadata & Header Sync in Python LiDAR Workflows overview --- title: "Metadata & Header Sync in Python LiDAR Workflows" description: "How to validate, correct, and reconcile LAS/LAZ header fields,...