Python LiDAR

Writing a LAS File from NumPy Arrays

TL;DR: Create laspy.LasHeader(point_format=6, version="1.4"), set header.scales to your precision (0.001 for millimetres) and header.offsets to a round value near the data minimum, add any extra dimensions, attach a CRS with header.add_crs(pyproj.CRS(...)), then build laspy.LasData(header), assign las.x, las.y, las.z and attributes, and las.write("out.laz"). Get the offset wrong and coordinates overflow the 32-bit integers LAS stores; get the scale wrong and you lose precision silently.

# Context and Motivation

This guide is part of laspy and NumPy Workflows for LAS Data. Points created in Python — simulated scenes for tests, photogrammetry exports, model outputs, results of a custom filter — usually need to end up in a LAS or LAZ file so that GIS software, viewers and PDAL can use them. laspy makes the mechanics easy, but LAS has two properties that punish carelessness. Coordinates are stored as 32-bit integers scaled and offset per axis, so the header’s scale and offset decide both precision and the representable range. And the point format fixes which fields exist; assigning colour to a format without colour simply is not possible.

Scale sets precision; offset sets range A number line of stored integers from minus 2.1 billion to plus 2.1 billion. With scale 0.001 and offset zero, a UTM northing of 4,471,380 metres needs 4.47 billion units, overflowing the range. With the same scale and an offset of 4,400,000 metres, the same northing needs only 71 million units and fits comfortably. −2.1 B 0 +2.1 B int32 range available for each stored coordinate offset 0: 4.47 B units — overflow offset 4,400,000: 71 M units — fits northing 4,471,380 m at scale 0.001

# Prerequisites and Assumptions

  • laspy 2.x (with lazrs for LAZ output) and pyproj for the CRS.
  • Coordinates in a projected CRS in metres (or a known unit), as float64.
  • A decision on point format: 6 for plain LAS 1.4, 7 with RGB, 8 with RGB and NIR.

# Step-by-Step Implementation

# Step 1 — Create a header with version and point format

laspy.LasHeader(point_format=6, version="1.4").

# Step 2 — Set scales from required precision

Millimetre precision is 0.001; centimetre is 0.01. Finer than needed costs nothing in LAS but makes LAZ compress less well.

# Step 3 — Set offsets near the data

Use a round number just below the minimum of each axis, for example np.floor(x.min() / 1000) * 1000. That keeps stored integers small, positive and well within range.

# Step 4 — Register extra dimensions and the CRS

header.add_extra_dim(laspy.ExtraBytesParams(name=..., type=...)) for each custom field, and header.add_crs(pyproj.CRS("EPSG:6347+5703")) for the spatial reference.

# Step 5 — Fill and write

Create laspy.LasData(header), assign coordinates and attributes, and write. laspy computes header counts and bounds on write.

# Complete Working Example

python
"""Write a LAS 1.4 / PDRF 7 LAZ from NumPy arrays with an extra dimension and a CRS."""
from __future__ import annotations

import laspy
import numpy as np
from pyproj import CRS


def write_las(path: str, x: np.ndarray, y: np.ndarray, z: np.ndarray, *,
              rgb: np.ndarray | None = None, classification: np.ndarray | None = None,
              confidence: np.ndarray | None = None, crs: str = "EPSG:6347+5703",
              scale: float = 0.001) -> None:
    fmt = 7 if rgb is not None else 6
    header = laspy.LasHeader(point_format=fmt, version="1.4")
    header.scales = np.array([scale, scale, scale])
    header.offsets = np.array([np.floor(x.min() / 1000) * 1000,
                               np.floor(y.min() / 1000) * 1000,
                               np.floor(z.min() / 100) * 100])
    if confidence is not None:
        header.add_extra_dim(laspy.ExtraBytesParams(
            name="confidence", type=np.float32, description="model score 0-1"))
    header.add_crs(CRS.from_user_input(crs))

    las = laspy.LasData(header)
    las.x, las.y, las.z = x, y, z
    n = x.size
    las.return_number = np.ones(n, dtype=np.uint8)
    las.number_of_returns = np.ones(n, dtype=np.uint8)
    las.classification = classification if classification is not None else np.ones(n, np.uint8)
    if rgb is not None:
        las.red, las.green, las.blue = (rgb.astype(np.uint16) * 257).T   # 8-bit to 16-bit
    if confidence is not None:
        las.confidence = confidence.astype(np.float32)
    las.write(path)


if __name__ == "__main__":
    rng = np.random.default_rng(11)
    n = 250_000
    x = 431_000 + rng.uniform(0, 500, n)
    y = 4_471_000 + rng.uniform(0, 500, n)
    z = 210 + 0.02 * (x - 431_000) + rng.normal(0, 0.05, n)
    rgb = rng.integers(0, 256, size=(n, 3))
    cls = np.where(rng.random(n) < 0.7, 2, 1).astype(np.uint8)
    conf = rng.random(n).astype(np.float32)
    write_las("out/synthetic_site.laz", x, y, z, rgb=rgb, classification=cls, confidence=conf)

    back = laspy.read("out/synthetic_site.laz")
    print(back.header.point_format.id, back.header.point_count, back.header.parse_crs().name)
    print("max xyz error:", float(np.abs(back.x - x).max()), float(np.abs(back.z - z).max()))

The final lines read the file back: the maximum coordinate error should be at most half the scale, 0.0005 m here.

Pick the format for the fields you have Three LAS 1.4 point formats as rows. Format 6 holds the base fields: coordinates, intensity, returns, classification, flags, scan angle, user data, point source ID and GPS time, 30 bytes. Format 7 adds red, green and blue for 36 bytes. Format 8 adds near-infrared for 38 bytes. Extra dimensions can be appended to any of them. PDRF 6 base fields + GPS time · 30 B PDRF 7 base fields + GPS time RGB 36 B PDRF 8 base fields + GPS time RGB NIR 38 B extra dimensions append after the standard fields in any format

# Key Parameter Table

Setting Recommended Why
version "1.4" Supports classes above 31, WKT CRS, extra bytes
point_format 6, 7 or 8 Match the fields you have
scales 0.001 or 0.01 Storage precision per axis
offsets rounded minimum per axis Keeps int32 values small and positive
add_crs compound CRS Horizontal and vertical reference in the WKT VLR
RGB 16-bit (× 257 from 8-bit) LAS colour is 16-bit per channel

# Verification

  • Round-trip error at most half the scale on every axis, as in the example.
  • Header counts and bounds match the arrays: back.header.point_count == n, back.header.mins close to the minima.
  • CRS present. back.header.parse_crs() returns the CRS you set, and pdal info --metadata shows it too.
  • Other readers agree. Open the file with PDAL (pdal info --summary) to confirm it is valid for tools other than laspy.

# Gotchas and Edge Cases

Overflow with zero offsets. At scale 0.001, any coordinate above about 2,147 km overflows. UTM northings and State Plane coordinates exceed that routinely. laspy raises an error in recent versions, but older code or other writers may wrap silently. Always set offsets.

Required fields left at zero. A return number of 0 is invalid in LAS; set return_number and number_of_returns to 1 for single-return data, as the example does. Some readers reject or misinterpret zeros.

8-bit colour. Writing 8-bit colour values into 16-bit fields makes images look almost black in viewers that expect 16-bit. Multiply by 257 to map 0–255 onto 0–65535.

What the scale throws away Three stored versions of the same coordinate 431522.3147 metres. At scale 0.001 it is stored as 431522.315, an error of 0.3 millimetres. At scale 0.01 it becomes 431522.31, an error of 4.7 millimetres. At scale 0.1 it becomes 431522.3, an error of 14.7 millimetres, coarser than many survey-grade deliverables allow. true value 431522.3147 m scale 0.001 431522.315 · error 0.3 mm scale 0.01 431522.31 · error 4.7 mm scale 0.1 14.7 mm

Classification above 31 in old formats. PDRF 0–5 store classification in 5 bits. Writing class 64 there wraps or fails. Use PDRF 6 or later for modern class ranges.

# Frequently Asked Questions

How do I create a LAS file from NumPy arrays in Python?

Create a laspy LasHeader with the version and point format, set scales and offsets, optionally add extra dimensions and a CRS, build a LasData from the header, assign the coordinate and attribute arrays, and call write with a .las or .laz filename.

What scale and offset should I use?

Choose the scale from the precision you need, commonly 0.001 or 0.01 metres, and set each offset to a round value just below that axis’s minimum. That preserves precision and keeps the stored integers well within the 32-bit range.

How do I add a custom attribute to a LAS file?

Register it on the header with add_extra_dim, giving a name, a NumPy type and optionally a description, then assign values to the field of that name before writing.

How do I set the CRS when writing with laspy?

Call add_crs on the header with a pyproj CRS object, ideally a compound CRS with both horizontal and vertical components. laspy writes it as a WKT record, which LAS 1.4 readers such as PDAL recognize.

Related in laspy and NumPy Workflows for LAS Data

↑ laspy and NumPy Workflows for LAS Data overview --- title: "laspy and NumPy Workflows for LAS Data" description: "Work with LAS and LAZ directly in Python using laspy 2.x and NumPy: scaled... Chunked Reading of Large LAS Files with laspy --- title: "Chunked Reading of Large LAS Files with laspy" description: "Process LAS and LAZ files larger than memory with laspy's... laspy vs PDAL: When to Use Which --- title: "laspy vs PDAL: When to Use Which" description: "A practical comparison of laspy and PDAL for Python LiDAR work: installation,... Reading LAS into NumPy with laspy --- title: "Reading LAS into NumPy with laspy" description: "Load LAS and LAZ files into NumPy arrays with laspy: scaled versus raw...