🔧 Fix LiDAR & point clouds CRS & projections LiDAR & Point Clouds · step 5 Notebook By · Updated

laspy Returns the Wrong Coordinates: Scale, Offset and CRS

On this page

Problem statement

Coordinates from a point cloud come out wrong in three distinct ways, and each has a different fix.

Wrong by a factor. Reading las.X instead of las.x gives the raw stored integer:

las.X.min() = -802212        raw int32
las.x.min() = -17342522.12   scaled

Wrong by a large constant. The header's offset is not the minimum; it is an arbitrary origin. Ignoring it puts the cloud millions of metres from where it belongs.

Right numbers, wrong units. The worst case, because nothing looks wrong. A survey stored in EPSG:3857 reports an extent of 1,613 × 1,613 m for a patch that is 841 × 841 m on the ground — a factor of 1.92 at 58.6° north. Every distance, area, density and slope derived from it is wrong by that factor or its square.

Quick answer

Check all three at load time:

import laspy
import numpy as np

las = laspy.read("nodes/4-3-13-7.laz")       # one node of the survey
h = las.header

print(f"scales  {h.scales}")
print(f"offsets {h.offsets}")
print(f"raw X   {las.X.min()} .. {las.X.max()}")
print(f"x       {las.x.min():.2f} .. {las.x.max():.2f}")

crs = h.parse_crs()
print(f"CRS     {crs.to_string() if crs else 'NOT SET'}")
if crs and crs.to_epsg() == 3857:
    import math
    from pyproj import Transformer
    lat = Transformer.from_crs(3857, 4326, always_xy=True) \
        .transform(np.mean(las.x), np.mean(las.y))[1]      # las.x is a view: no .mean()
    print(f"! Web Mercator at {lat:.1f}° — distances inflated by "
          f"{1 / math.cos(math.radians(lat)):.2f}x, areas by "
          f"{1 / math.cos(math.radians(lat)) ** 2:.2f}x")
scales  [0.01 0.01 0.01]
offsets [-1.73345e+07  8.07826e+06  4.86000e+02]
raw X   -802212 .. -640954
x       -17342522.12 .. -17340909.54
CRS     EPSG:3857
! Web Mercator at 58.6° — distances inflated by 1.92x, areas by 3.68x
Three coordinate failures: raw integers used unscaled, a missing offset, and correct numbers in a distorted CRS. three failures, one of them invisiblesymptomcausecoordinates are plain integersread las.X instead of las.xthe cloud is millions of metresawaythe header offset was ignoredthe numbers look completely fineWeb Mercator — every distanceinflatedThe first two look obviously wrong; the third is the expensive one spatialworkflow.io
The first two look obviously wrong. The third looks fine and is the expensive one.

Scaled attributes, no hand reconstruction, check the CRS

1. Use the scaled attributes

x, y, z = np.asarray(las.x), np.asarray(las.y), np.asarray(las.z)

Never las.X. The uppercase attributes are the stored integers, and the relationship is x = X * scale + offset. With a 0.01 scale that is a factor of 100; millimetre-scaled files give 1,000.

The symptom is coordinates that are plausible integers in the hundreds of thousands with no decimal part.

2. Never reconstruct coordinates by hand

x = las.X * las.header.scales[0] + las.header.offsets[0]   # correct, and pointless
x = las.X * 0.01                                            # wrong: no offset
x = las.X * las.header.scales[0] + las.header.mins[0]       # wrong: mins is not offset

The offset and the minimum are different things. The offset is an arbitrary origin chosen by the writer to keep the integers in range; the minimum is a property of the data. Confusing them displaces the cloud by their difference.

3. Check whether the CRS is metric and undistorted

Three CRS states, in increasing order of subtlety:

  • Missing. parse_crs() returns None. Common in older files; the CRS must come from the survey metadata.
  • Geographic. Coordinates in degrees. Every distance is nonsense, and obviously so once you look at the extent.
  • Conformal but distorted. Web Mercator is projected and metric, so nothing looks wrong — and its scale factor is 1/cos(latitude).

The third is the dangerous one. It is also common, because point clouds served for web viewers are frequently stored in EPSG:3857.

4. Reproject before any metric computation

from pyproj import Transformer

transformer = Transformer.from_crs(3857, 32605, always_xy=True)
x, y = transformer.transform(las.x, las.y)

The effect on this survey, measured on all of tile.laz:

from scipy.spatial import ConvexHull

full = laspy.read("tile.laz")
wx, wy = np.asarray(full.x), np.asarray(full.y)
ux, uy = transformer.transform(wx, wy)
for label, px, py in [("Web Mercator", wx, wy), ("UTM 5N", ux, uy)]:
    area = ConvexHull(np.c_[px, py]).volume           # the area the points cover
    print(f"{label:13} extent {px.max() - px.min():6,.0f} x {py.max() - py.min():6,.0f} m   "
          f"area {area / 1e4:6.1f} ha   {len(px) / area:5.1f} points per m²")
Web Mercator  extent  1,613 x  1,613 m   area  260.1 ha     5.0 points per m²
UTM 5N        extent    877 x    876 m   area   70.9 ha    18.3 points per m²

The UTM extent is a bounding box: the survey's square patch is turned 2.4° against the UTM grid there, so its box is 877 × 876 m while the patch itself is 841 m on a side.

Note that z is unaffected. Web Mercator distorts the horizontal plane only, so a cloud read without reprojection has correct heights and wrong horizontal distances — which produces slopes that are too gentle by exactly the scale factor.

5. Sanity-check the extent against something you know

The cheapest check in the whole workflow: does the reported extent match the area you believe you have?

A 1 km tile reporting 1.9 km, or a 500 ha survey reporting 1,800 ha, is a scale-factor problem. It takes ten seconds and catches the failure that no exception will.

The Web Mercator scale factor rising from 1 at the equator to 1.92 at 58.6 degrees and 2.92 at 70 degrees, with area inflated by its square. Web Mercator distance inflation, by latitude0° (equator)1.00×40°1.31×58.6° (this survey)1.92×70°2.92×03.2Areas inflate by the square: 3.68× at 58.6°, so density falls by the same spatialworkflow.io
The distortion grows with latitude and applies to both axes, so area error is the square of distance error.

A load-time guard and the distortion measured

The data used on this page

The USGS 3DEP lidar survey of Brooks Camp, Katmai National Park, Alaska (2012), from the USGS Entwine Point Tile archive on AWS, AK_BrooksCamp_2012, which stores every survey in Web Mercator. The setup downloads the 287 LAZ files under octree node 4-3-13-7 into nodes/ and joins them into tile.laz; the quick answer reads the top node of that branch, nodes/4-3-13-7.laz, on its own.

import json
import os
import urllib.request

import laspy

EPT = "https://s3-us-west-2.amazonaws.com/usgs-lidar-public/AK_BrooksCamp_2012"
ROOT = (4, 3, 13, 7)                                     # depth, x, y, z of the octree node


def under(key):
    d, x, y, z = map(int, key.split("-"))
    step = 2 ** (d - ROOT[0])
    return d >= ROOT[0] and (x // step, y // step, z // step) == ROOT[1:]


hierarchy = json.load(urllib.request.urlopen(f"{EPT}/ept-hierarchy/0-0-0-0.json"))
nodes = sorted(k for k in hierarchy if under(k))
os.makedirs("nodes", exist_ok=True)
for key in nodes:
    path = f"nodes/{key}.laz"
    if not os.path.exists(path):
        urllib.request.urlretrieve(f"{EPT}/ept-data/{key}.laz", path + ".part")
        os.replace(path + ".part", path)

if not os.path.exists("tile.laz"):                       # the whole subtree as one file
    with laspy.open(f"nodes/{nodes[0]}.laz") as first:
        header = first.header
    with laspy.open("tile.laz", mode="w", header=header) as writer:
        for key in nodes:
            writer.write_points(laspy.read(f"nodes/{key}.laz").points)
with laspy.open("tile.laz") as reader:
    print(len(nodes), "nodes,", f"{reader.header.point_count:,} points,",
          f"{sum(os.path.getsize(f'nodes/{k}.laz') for k in nodes) / 1e6:.1f} MB of nodes")
287 nodes, 12,968,770 points, 82.8 MB of nodes

Example 1 — a load-time guard

import math
import laspy
import numpy as np
from pyproj import CRS, Transformer
from scipy.spatial import ConvexHull


def load_metric(path, target_crs=None, max_scale_factor=1.02):
    """Load a cloud in a CRS where a metre is a metre, or say why not."""
    las = laspy.read(path)
    x, y, z = np.asarray(las.x), np.asarray(las.y), np.asarray(las.z)
    crs = las.header.parse_crs()

    if crs is None:
        raise ValueError(f"{path} has no CRS. Supply it from the survey "
                         "metadata; nothing metric is safe without it.")

    if not crs.is_projected:
        raise ValueError(f"{path} is in {crs.to_string()}, a geographic CRS. "
                         "Coordinates are degrees, not metres.")

    if crs.to_epsg() == 3857:
        lat = Transformer.from_crs(crs, 4326, always_xy=True) \
            .transform(float(x.mean()), float(y.mean()))[1]
        k = 1 / math.cos(math.radians(lat))
        if k > max_scale_factor:
            if target_crs is None:
                raise ValueError(
                    f"Web Mercator at {lat:.1f}°: distances inflated {k:.2f}x, "
                    f"areas {k ** 2:.2f}x. Pass target_crs to reproject.")
            transformer = Transformer.from_crs(crs, target_crs, always_xy=True)
            x, y = transformer.transform(x, y)
            print(f"  reprojected {crs.to_string()} -> {target_crs} "
                  f"(scale factor was {k:.2f})")
            crs = CRS.from_user_input(target_crs)

    area = ConvexHull(np.c_[x, y]).volume                  # not the bounding box
    print(f"  {len(x):,} points over a {x.max() - x.min():.0f} x "
          f"{y.max() - y.min():.0f} m box, covering {area / 1e4:.1f} ha")
    print(f"  density {len(x) / area:.1f} per m², z {z.min():.1f}..{z.max():.1f} m")
    return x, y, z, np.asarray(las.classification), crs


try:
    load_metric("tile.laz")
except ValueError as err:
    print(f"  ValueError: {err}")
x, y, z, classification, crs = load_metric("tile.laz", target_crs="EPSG:32605")
  ValueError: Web Mercator at 58.6°: distances inflated 1.92x, areas 3.68x. Pass target_crs to reproject.
  reprojected EPSG:3857 -> EPSG:32605 (scale factor was 1.92)
  12,968,770 points over a 877 x 876 m box, covering 70.9 ha
  density 18.3 per m², z 15.9..203.1 m

Raising rather than warning is deliberate. A warning in a batch job is a line in a log nobody reads; an exception stops the run while the mistake is still cheap.

Example 2 — measuring the distortion instead of assuming it

import numpy as np
from pyproj import Geod, Transformer

GEOD = Geod(ellps="WGS84")


def measure_distortion(x, y, source_crs):
    """Compare a projected distance with the true geodesic distance."""
    to_wgs = Transformer.from_crs(source_crs, 4326, always_xy=True)

    x0, x1 = float(x.min()), float(x.max())
    y0 = float(y.mean())
    lon0, lat0 = to_wgs.transform(x0, y0)
    lon1, lat1 = to_wgs.transform(x1, y0)

    projected = x1 - x0
    _, _, true = GEOD.inv(lon0, lat0, lon1, lat1)

    print(f"  projected {projected:10.1f} m")
    print(f"  geodesic  {true:10.1f} m")
    print(f"  scale factor {projected / true:.4f} "
          f"(areas {(projected / true) ** 2:.4f})")
    return projected / true


k = measure_distortion(wx, wy, "EPSG:3857")
k_utm = measure_distortion(np.asarray(x), np.asarray(y), "EPSG:32605")
  projected     1612.6 m
  geodesic       843.1 m
  scale factor 1.9128 (areas 3.6588)
  projected      877.0 m
  geodesic       877.1 m
  scale factor 0.9999 (areas 0.9998)

This works for any CRS, not just Web Mercator, and it is the check to run when you are unsure. A scale factor within a few parts per thousand of 1 is fine; 1.92 is not.

Example 3 — reprojecting and writing a corrected file

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


def reproject_cloud(source_path, out_path, target_epsg, scale=0.01):
    """Write a copy in a projected CRS, with a header that matches."""
    las = laspy.read(source_path)
    source_crs = las.header.parse_crs()
    if source_crs is None:
        raise ValueError("source has no CRS")

    transformer = Transformer.from_crs(source_crs, target_epsg, always_xy=True)
    x, y = transformer.transform(np.asarray(las.x), np.asarray(las.y))
    z = np.asarray(las.z)

    header = laspy.LasHeader(version=las.header.version,
                             point_format=las.header.point_format)
    header.scales = np.array([scale, scale, scale])
    header.offsets = np.array([np.floor(x.min()), np.floor(y.min()),
                               np.floor(z.min())])
    header.add_crs(CRS.from_epsg(target_epsg))           # a pyproj CRS, not a bare code

    out = laspy.LasData(header)
    # las.points.copy() is a plain packed record with no scales, and assigning
    # out.x to it fails; wrap the raw array in a record that knows the new header
    out.points = laspy.ScaleAwarePointRecord(las.points.array.copy(), header.point_format,
                                             header.scales, header.offsets)
    out.x, out.y, out.z = x, y, z
    out.write(out_path)

    print(f"  {len(x):,} points -> {out_path} in EPSG:{target_epsg}")
    print(f"  extent {x.max() - x.min():.0f} x {y.max() - y.min():.0f} m")
    return out_path


reproject_cloud("tile.laz", "tile_utm5n.laz", 32605)
with laspy.open("tile_utm5n.laz") as check:
    print(" ", check.header.parse_crs().to_epsg(), check.header.offsets, check.header.point_count)
moved = laspy.read("tile_utm5n.laz")
print(f"  largest difference from the reprojected coordinates: "
      f"{np.abs(np.asarray(moved.x) - x).max() * 1000:.1f} mm (1 cm scale)")
  12,968,770 points -> tile_utm5n.laz in EPSG:32605
  extent 877 x 876 m
  32605 [3.376760e+05 6.494627e+06 1.500000e+01] 12968770
  largest difference from the reprojected coordinates: 2.7 mm (1 cm scale)

Setting the offsets from the reprojected minima matters. Keeping the old offsets, which were near 17 million, would push the new UTM coordinates far outside the range a scaled 32-bit integer can hold, and the write would either fail or wrap.

Why the offset exists

Why the offset exists at all

A signed 32-bit integer spans about ±2.1 billion. At a scale of 0.01 that is ±21 million units of the coordinate system — enough for any single CRS, but only if the origin is near the data.

The offset moves that window. Here it is −17,334,459, placing the representable range around the Web Mercator easting of the survey. Ignore it and every coordinate is displaced by 17 million metres.

The offset is chosen by whoever wrote the file and is not the data minimum, though writers often pick something close to it.

Why Web Mercator is worse than a geographic CRS

A geographic CRS announces itself. Coordinates are between −180 and 90, distances come out absurd, and the mistake is caught in the first plot.

Web Mercator produces metre-like numbers in a projected CRS with a correct EPSG code. The file is valid, is_projected is True, and the numbers are the right order of magnitude. Only the scale is wrong, by 1/cos(latitude) — 1.92 at 58.6° north, 1.31 at 40°, 2.92 at 70°.

Because the distortion is a pure scale, everything remains internally consistent: shapes are right, angles are right, and only measurements against the real world are wrong.

Why heights are unaffected and slope is not

Web Mercator is a horizontal projection. It does nothing to z, so elevations are correct.

That means slope — a vertical difference over a horizontal distance — is understated by exactly the scale factor. A true 30° slope reads as 17° at 58.6° north. Because the elevations are right and the shape is right, nothing about the hillshade looks wrong.

Why the extent check is the best defence

Every one of these failures changes the reported extent, and the extent is the one number you can usually check against outside knowledge: the tile size in the survey specification, the width of a known feature, a measurement from a map.

Ten seconds spent asking "is this patch really 1.6 km across?" catches all three failure modes, and nothing else does.

Checks on a point cloud at load time. the cheapest test in the whole workflowCRS is presentparse_crs()CRS is projectednot degreesCRS is not Web Mercatoror reprojectextent matches what you expect841 m, not 1,613scaled coordinates usedlas.x, not las.XAsking "is this patch really 1.6 km across?" catches all three failures spatialworkflow.io
The extent check takes ten seconds and catches every one of the failures above.

Raw values, offsets and missing CRSs

  • las.x is scaled; las.X is raw. Factor of 100 or 1,000.
  • The offset is not the minimum. It is an arbitrary origin near the data.
  • parse_crs() can return None. Get the CRS from the survey metadata.
  • Web Mercator is projected and distorted. is_projected being True is not enough.
  • Heights are unaffected by a horizontal projection — so slope is wrong while elevation is right.
  • Reset the offsets after reprojecting, or the new coordinates overflow the integer range.
  • Vertical datums differ too. A cloud in ellipsoidal heights is not comparable with one in orthometric heights.
  • Check the extent against something you know. It is the cheapest test available.

FAQ

Why are my laspy coordinates so large?

Either you read las.X (the raw integer) instead of las.x, or the file is in a CRS with large coordinate values such as Web Mercator. Print both and compare.

What is the difference between scale and offset?

coordinate = raw_integer * scale + offset. The scale sets the precision; the offset moves the representable range near the data. The offset is not the data minimum.

My point cloud is in EPSG:3857 — is that a problem?

For anything metric, yes. At 58.6° north distances are inflated by 1.92 and areas by 3.68. Reproject to a suitable projected CRS first.

Why does my point density look too low?

Because it was computed in Web Mercator. The same survey reads as 5.0 points per square metre there and 16.9 in UTM.

Are my elevations affected by the projection?

No. Web Mercator distorts the horizontal plane only. But slope, being vertical over horizontal, is understated by the scale factor.

What if the file has no CRS?

parse_crs() returns None, and you must supply the CRS from the survey metadata. Nothing metric is safe until you do.

Why did my reprojected file fail to write?

Probably the header offsets were still set for the old CRS, pushing the new coordinates outside the range a scaled 32-bit integer can represent. Reset the offsets from the new minima.