🔧 Fix LiDAR & point clouds CRS & projections LiDAR & Point Clouds · step 5 of 14 Updated

laspy Returns the Wrong Coordinates: Scale, Offset and CRS

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 877 × 876 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("tile.laz")
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(las.x.mean(), las.y.mean())[1]
    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.
The first two look obviously wrong. The third looks fine and is the expensive one.

Step-by-step solution

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:

                  Web Mercator     UTM 5N
extent               1,613 m       877 m
area                 260.1 ha     70.7 ha
point density     5.0 per m²   16.9 per m²

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.
The distortion grows with latitude and applies to both axes, so area error is the square of distance error.

Code examples

Example 1 — a load-time guard 

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


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 = (x.max() - x.min()) * (y.max() - y.min())
    print(f"  {len(x):,} points over {x.max() - x.min():.0f} x "
          f"{y.max() - y.min():.0f} m = {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

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
  projected     1613.0 m
  geodesic       841.2 m
  scale factor 1.9175 (areas 3.6768)

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 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(target_epsg)

    out = laspy.LasData(header)
    out.points = las.points.copy()
    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

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.

Explanation

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,500, 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.

Five load-time checks on a point cloud, ending with whether the extent matches what you expect.
The extent check takes ten seconds and catches every one of the failures above.

Edge cases or notes

  • 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.