Fixing LAS Scale and Offset Precision Loss

A survey delivered at 8 mm accuracy comes back from the client’s software with every coordinate ending in a round centimetre. Or worse: the cloud opens with points spread across thousands of kilometres, a few of them near the site and the rest in the sea. Both symptoms come from the same two header fields, and neither produces an error when it is written.

The underlying cause is that LAS does not store coordinates as floating-point numbers at all. It stores three signed 32-bit integers per point, and the real coordinate is reconstructed from them through a per-axis scale and offset in the header. Get those two wrong and the file is arithmetically incapable of holding the data — a fact that is invisible until somebody measures something.

This page covers detecting both failures before writing, choosing header values from the data, and what can and cannot be recovered afterwards. It is the detailed version of the header discussion in point cloud formats and interchange in Python.

# The arithmetic, exactly

For each axis, the stored integer is round((real − offset) / scale). Reconstruction is the inverse. Two constraints follow.

The quantisation is the scale. A scale of 0.01 rounds every coordinate to the nearest centimetre. This is not an approximation that averages out across a cloud — it is a hard lattice, and a distance measured between two points on it can be out by a full scale unit.

The range is the signed 32-bit limit, about ±2.147 billion. The largest representable displacement from the offset is therefore 2.147e9 × scale. At a scale of 0.001 that is ±2,147 km, which is ample — but only measured from the offset. With the offset left at zero and a projected northing of 5.6 million metres, the required integer is 5.6 billion, which overflows and wraps to a negative value. The point lands on the other side of the world.

Figure 1 — The lattice, at the three scales in common use.

# Minimal reproducible solution

Compute both header fields from the data, and refuse rather than degrade.

import numpy as np

INT32_SAFE = 2_147_483_000          # a little under the true limit, for rounding


def safe_header_geometry(xyz: np.ndarray, *, precision: float = 0.001,
                         snap_offset: float = 1.0) -> dict:
    """Header scale and offset that hold `precision` without overflow.

    The offset is snapped to a whole metre so it is exactly representable in
    the header's double and reads cleanly to a human. The check is done per
    axis, because a corridor survey can be fine in easting and overflow in
    northing.
    """
    lo, hi = np.asarray(xyz.min(axis=0)), np.asarray(xyz.max(axis=0))
    offset = np.round((lo + hi) / 2.0 / snap_offset) * snap_offset

    reach = np.maximum(np.abs(lo - offset), np.abs(hi - offset))
    needed_scale = reach / INT32_SAFE

    if (needed_scale > precision).any():
        axis = "xyz"[int(np.argmax(needed_scale))]
        raise ValueError(
            f"axis {axis} spans {2 * reach.max():.0f} m about its centre, which needs "
            f"scale {needed_scale.max():.6f} m — coarser than the requested "
            f"{precision} m. Tile the cloud instead of coarsening it.")

    return {"scale": (precision,) * 3,
            "offset": tuple(float(v) for v in offset),
            "headroom": float(precision / needed_scale.max())}

The headroom figure — how many times more range is available than needed — is worth logging. A value of 1.2 means the survey is close to the limit and a slightly larger extent next month will fail; a value of 400 means there is room to go finer if the survey warrants it.

# Detecting the failure in a file somebody else wrote

Incoming clouds are frequently the problem. Two checks identify both faults from the header alone, without reading a point.

import laspy
import numpy as np


def audit_header(path: str, *, required_precision: float = 0.005) -> dict:
    """Judge an existing file's header against the precision it should carry."""
    with laspy.open(path) as fh:
        h = fh.header
        scales = np.asarray(h.scales)
        offsets = np.asarray(h.offsets)
        mins, maxs = np.asarray(h.mins), np.asarray(h.maxs)

    problems = []
    if (scales > required_precision).any():
        problems.append(
            f"scale {scales.tolist()} is coarser than {required_precision} m — "
            "coordinates in this file are quantised and cannot be improved")

    reach = np.maximum(np.abs(mins - offsets), np.abs(maxs - offsets))
    counts = reach / scales
    if (counts > 2_147_483_647).any():
        problems.append("stored integers exceed int32 — coordinates have wrapped")
    elif (counts > 2.0e9).any():
        problems.append("within 7 % of the int32 limit — a slightly larger "
                        "survey will overflow")

    if (np.abs(offsets) < 1.0).all() and (np.abs(mins) > 1.0e5).any():
        problems.append("offset is at zero for projected coordinates — "
                        "this file is one scale change away from wrapping")

    return {"scales": scales.tolist(), "offsets": offsets.tolist(),
            "max_integer": float(counts.max()), "problems": problems}

The third check — offset at zero with large projected coordinates — flags files that are currently correct and fragile. They work at 0.01 scale and break the moment anybody re-writes them at 0.001, which is exactly what a precision-conscious pipeline will do.

Figure 3 — Scale is a range-against-precision trade, and the default is rarely right.

# Edge-case matrix

The geographic-coordinates row deserves emphasis. A degree is roughly 111 km, so a scale of 0.001 degrees quantises to about 100 m. Storing unprojected coordinates in LAS at a metre-oriented scale is a mistake that produces a file which is obviously wrong once plotted and entirely plausible in the header.

# Verification snippet

Verification is a round trip: write, read back, compare against the array that was written.

import laspy
import numpy as np


def roundtrip_error(path: str, original_xyz: np.ndarray) -> dict:
    """Largest coordinate error introduced by the write, per axis."""
    las = laspy.read(path)
    back = np.column_stack([las.x, las.y, las.z])
    if back.shape != original_xyz.shape:
        raise ValueError(f"point count changed: {original_xyz.shape[0]} → {back.shape[0]}")

    err = np.abs(back - original_xyz)
    scales = np.asarray(las.header.scales)
    worst = err.max(axis=0)

    return {"worst_error_mm": (worst * 1000).tolist(),
            "scale_mm": (scales * 1000).tolist(),
            # A correct write errs by at most half a scale unit per axis.
            "ok": bool((worst <= scales / 2 + 1e-12).all())}

Comparing the observed error against half the scale unit is what makes this a real test rather than a tolerance guess. A correct encoding cannot do worse than half a lattice step; anything larger means something other than quantisation is happening — most often an overflow that has wrapped a handful of points.

Figure 2 — Why the offset is not an optional nicety. It is what brings the numbers into range at any useful scale.

# Choosing a precision to standardise on

Rather than deciding per job, most operations are better served by standardising on one scale across the fleet and only deviating with a reason. Three considerations point at the same answer.

A drone survey’s own accuracy is rarely better than about 10 mm in the vertical, even with dense ground control. A scale of 1 mm is therefore an order of magnitude finer than the data, which is the right relationship: the encoding should never be the limiting factor, and one order of magnitude of headroom costs nothing. Going to 0.1 mm buys nothing measurable and cuts the representable reach to 215 km, which is still ample but narrows the margin for no benefit.

A scale of 1 mm also keeps the numbers human-readable. Coordinates print as three decimal places, which matches how survey reports quote them, and a technician comparing a file against a coordinate list is not mentally converting anything.

And it is safe for every projected extent a drone survey produces. With the offset at the data centre, 1 mm supports a survey 4,000 km across — which no single flight will ever approach — so the refuse-rather-than-degrade branch effectively never fires for legitimate data. When it does fire, it is telling you something real: usually that a merge combined two sites, or that an outlier at the edge of the solar system survived the noise filter.

Write the chosen scale into the pipeline configuration rather than into each writer call, and audit incoming third-party files against it. That single convention removes an entire class of quiet delivery failure.

# When to escalate

• The delivered file is already quantised and the source is gone. Nothing recovers it. Say so plainly and re-derive from the reconstruction if it still exists; an interpolated “improvement” is fabrication.
• A survey genuinely needs finer than the int32 range allows. That is a tiling decision. A single file cannot hold a 40 km corridor at 0.1 mm precision, and pretending otherwise produces wrapped coordinates.
• A client’s software rewrites the headers on import. Some packages normalise scale and offset on load, quietly coarsening the data. Ask what it writes before delivering at a fine scale, and supply COPC or LAZ if it preserves them better.

# Related

• Point Cloud Formats and Interchange in Python
• Fixing SRS Lost When Writing LAS Files
• Converting LAS to COPC for Cloud Streaming

← Point Cloud Formats and Interchange in Python