Choosing a Projected CRS for a LiDAR Project

TL;DR: Process in a metric projected CRS with a current datum realization — for example NAD83(2011) / UTM (EPSG:6330–6348) in the US, ETRS89 / UTM (EPSG:25828–25838) in Europe, or the national grid your client uses — combined with the official vertical datum as a compound CRS. Find the UTM zone from the project centroid with pyproj.database.query_utm_crs_info, and keep a project that straddles two zones in one zone rather than splitting it.

# Context and Motivation

This guide is part of Coordinate Reference Systems. Everything that measures distances or areas in LiDAR processing — point spacing, SMRF windows, raster cell sizes, buffer widths, slope — assumes a projected CRS in metres with small distortion. Geographic coordinates in degrees break all of it, and a projection chosen carelessly introduces scale errors or forces data across zone boundaries. The choice is usually easy, but it should be made once, deliberately, at the start of a project, and written into every file as a compound CRS so that heights are as well defined as positions.

# Prerequisites and Assumptions

• The project boundary, in any CRS, to compute a centroid and extent.
• pyproj 3.x.
• Knowledge of the client’s requirements and of the official vertical datum and geoid model for the region.

# Step-by-Step Implementation

# Step 1 — Check for a mandated CRS

Contracts and national programmes often fix the delivery CRS. If so, process in it too, unless it is geographic, in which case process in a projected CRS and reproject once at delivery.

# Step 2 — Find the UTM zone

Query pyproj for UTM CRSs whose area of use contains the project centroid, filtered by datum name.

# Step 3 — Check the extent against the zone

UTM zones are 6° wide. A project that crosses a zone boundary by a few kilometres should stay in the zone containing most of it; scale error just outside a zone edge is small and far better than splitting the project.

# Step 4 — Choose the vertical datum

Use the official vertical datum for the region — NAVD88 (EPSG:5703) in the conterminous US today, EVRF2007 or a national height system in Europe — and record the geoid model used.

# Step 5 — Write the compound code everywhere

EPSG:6344+5703 style codes in writers.las (a_srs), writers.gdal and every output’s metadata.

# Complete Working Example

"""Pick a UTM CRS for a project boundary and report scale distortion across it."""
from __future__ import annotations

import geopandas as gpd
from pyproj import CRS, Proj
from pyproj.aoi import AreaOfInterest
from pyproj.database import query_utm_crs_info


def pick_utm(boundary: gpd.GeoDataFrame, datum_name: str = "NAD83(2011)") -> CRS:
    geo = boundary.to_crs("EPSG:4326")
    lon, lat = geo.geometry.union_all().centroid.coords[0]
    infos = query_utm_crs_info(datum_name=datum_name,
                               area_of_interest=AreaOfInterest(lon, lat, lon, lat))
    if not infos:
        raise ValueError(f"no {datum_name} UTM zone found at {lon:.3f}, {lat:.3f}")
    return CRS.from_epsg(infos[0].code)


def scale_report(boundary: gpd.GeoDataFrame, crs: CRS) -> dict:
    geo = boundary.to_crs("EPSG:4326").geometry.union_all()
    minx, miny, maxx, maxy = geo.bounds
    proj = Proj(crs)
    factors = []
    for lon in (minx, (minx + maxx) / 2, maxx):
        for lat in (miny, maxy):
            factors.append(proj.get_factors(lon, lat).meridional_scale)
    return {"crs": crs.name, "epsg": crs.to_epsg(),
            "scale_min": round(min(factors), 6), "scale_max": round(max(factors), 6),
            "max_distortion_ppm": round(max(abs(f - 1) for f in factors) * 1e6, 1)}


if __name__ == "__main__":
    area = gpd.read_file("project/boundary.gpkg")
    utm = pick_utm(area)
    print(scale_report(area, utm))
    compound = f"EPSG:{utm.to_epsg()}+5703"
    print("write outputs with a_srs =", compound, "->", CRS.from_user_input(compound).name)

A typical result for a county near the central meridian:

{'crs': 'NAD83(2011) / UTM zone 15N', 'epsg': 6344, 'scale_min': 0.999602, 'scale_max': 0.999701, 'max_distortion_ppm': 398.0}
write outputs with a_srs = EPSG:6344+5703 -> NAD83(2011) / UTM zone 15N + NAVD88 height

A scale factor of 0.9996 at the central meridian is by design; 400 ppm is 4 cm per 100 m, negligible for point spacing and cell sizes but worth knowing for survey-grade distance comparisons.

# Key Parameter Table

Confirm each code in your PROJ database before use; national agencies periodically publish new realizations.

# Verification

• Centroid lands in the zone. Transform the centroid to the chosen CRS; the easting should fall well within 166,000–834,000 m for UTM.
• Round trip. Transform a few boundary vertices to the CRS and back; errors should be sub-millimetre.
• Vertical present. CRS.from_user_input(compound).sub_crs_list has two members, projected and vertical.

# Gotchas and Edge Cases

Old datum realizations. Plain NAD83 (EPSG:26915 etc.) and NAD83(2011) differ by up to a metre or more in places. Choose the realization that matches your control survey, and use it consistently.

Feet in state plane. Many state plane CRSs have both metre and US-survey-foot variants. Using a feet CRS for processing complicates every metric parameter; process in metres and convert at delivery if the client wants feet, as in reprojecting State Plane feet to metres.

Projects spanning two zones. Splitting a project across zones creates a seam in every product. Keep one zone, accept the slightly larger scale factor beyond the edge, or use a custom transverse Mercator centred on the project for very wide areas.

Geographic processing. Some public data is distributed in EPSG:4326. Reproject to a projected CRS before any metric processing; filters that take distances in metres interpret degrees as metres without complaint.

# Frequently Asked Questions

Which UTM zone should I use?

The zone containing the project centroid, which pyproj can find with query_utm_crs_info. If the project crosses a zone boundary, keep it in the zone containing most of the area rather than splitting it.

Should I use UTM or a national grid?

Use what your client and local users work in. National and state grids are designed to minimize distortion locally and match existing mapping; UTM is a sound default when no local grid is standard.

Why use a compound CRS?

A compound CRS declares both the horizontal system and the vertical datum. Without the vertical part, heights are ambiguous, and datum or geoid mismatches become impossible to detect or correct later.

Can I process LiDAR in latitude and longitude?

Not sensibly. Point spacing, window sizes, raster resolution and slopes all assume metric coordinates. Reproject to a projected CRS first and, if needed, back to geographic at delivery.

# Related

• Coordinate Reference Systems — CRS handling overview
• Reading the CRS from LAS WKT and GeoTIFF Keys — what a file declares
• Inspecting PROJ Transformations Before Reprojecting — getting into the chosen CRS accurately
• Setting a Vertical CRS on a Point Cloud — the vertical half
• Reprojecting Point Clouds from UTM to WGS84 — delivering in geographic coordinates