LiDAR Point Clouds Explained: Returns, Classes and Intensity

Problem statement

A lidar point cloud is not a 3D picture. It is a table of pulse measurements, and almost every question people ask of it — how tall is that tree, where is the ground, how dense is the survey — depends on columns other than X, Y and Z.

Loaded from a real 3DEP survey over Brooks Camp, Alaska:

12,968,770 points over 841 x 841 m
18.3 points per square metre
59.5% first returns, 36.1% second, 4.4% third, 0.1% fourth
76.5% of points came from a pulse that produced more than one return
22.9% classified as ground; 77.1% left unclassified

That last line is the one that ruins most first attempts. This survey has no vegetation classes at all — only ground and unclassified. Any script assuming class 5 means "high vegetation" produces an empty array here.

Quick answer

import laspy
import numpy as np

las = laspy.read("tile.laz")
print(f"{len(las.points):,} points, LAS {las.header.version}, "
      f"point format {las.header.point_format.id}")
print("dimensions:", [d.name for d in las.point_format.dimensions])

classes, counts = np.unique(las.classification, return_counts=True)
for code, count in zip(classes, counts):
    print(f"  class {code}: {count:,} ({count / len(las.points):.1%})")

12,968,770 points, LAS 1.2, point format 1
dimensions: ['X', 'Y', 'Z', 'intensity', 'return_number', 'number_of_returns', 'scan_direction_flag', 'edge_of_flight_line', 'classification', 'synthetic', 'key_point', 'withheld', 'scan_angle_rank', 'user_data', 'point_source_id', 'gps_time', 'OriginId']
  class 1: 10,003,338 (77.1%)
  class 2: 2,965,430 (22.9%)
  class 7: 2 (0.0%)

Always print the class histogram before writing any code that depends on a class.

What a row is, returns, classifications

1. Understand what one row is

A row is one return: a detectable echo from one outgoing laser pulse. A pulse over bare ground gives one return. A pulse through a canopy can give several, as energy reflects off leaves, branches and finally the ground.

That is why return_number and number_of_returns are two separate columns. return_number = 2, number_of_returns = 3 means "the second of three echoes from this pulse".

Measured here, 76.5% of points came from multi-return pulses — an unusually high figure that says this is a well-vegetated landscape scanned by a sensor that resolves closely spaced echoes.

2. Use returns to reason about structure

• First returns are the top of whatever the pulse hit — canopy, roof, ground in the open.
• Last returns are the deepest thing the pulse reached. Over vegetation that is often, but not always, the ground.
• Intermediate returns are canopy interior, which is what makes lidar useful for forest structure.
• Single returns (number_of_returns == 1) came from a surface that stopped the pulse outright: bare ground, water, a roof.

Measured relationship between last returns and classified ground:

returns = np.asarray(las.return_number)
last = returns == np.asarray(las.number_of_returns)
ground = np.asarray(las.classification) == 2
print(f"ground points that are last returns:     {last[ground].mean():6.1%}")
print(f"last returns that are classified ground: {ground[last].mean():6.1%}")

ground points that are last returns:     100.0%
last returns that are classified ground:  38.4%

Every ground point is a last return, and only two in five last returns are ground. "Last return" is a necessary condition for ground, not a sufficient one — which is exactly why classification is a separate, harder step.

3. Check the classification before trusting it

The LAS standard defines class codes: 2 is ground, 3–5 low/medium/high vegetation, 6 building, 9 water, 7 and 18 noise.

Defining them does not mean a survey populated them. This one has 22.9% ground and 77.1% class 1, "unclassified" — meaning the vendor ran a ground filter and stopped.

if 5 not in set(np.unique(las.classification)):
    print("no vegetation class — derive height from returns, not from classes")

no vegetation class — derive height from returns, not from classes

4. Treat intensity as relative, not physical

intensity is the amplitude of the returned echo. It correlates with surface reflectance at the laser's wavelength, which makes it tempting to use like a band of imagery.

It is not calibrated. It varies with range, incidence angle, atmospheric conditions and the sensor's automatic gain — so intensity is not comparable between flight lines, let alone between surveys. Use it within one flight line, for relative contrast, and check point_source_id to see how many flight lines you are mixing.

5. Look at point_source_id and gps_time

point_source_id identifies the flight line. Where two lines overlap, density doubles and intensity jumps, which produces stripes in any density or intensity raster.

gps_time gives the acquisition time, which is how you separate passes and how you detect that a "single survey" was actually flown across three weeks with a metre of leaf-on change in between.

A survey report before any analysis

The data used on this page

The USGS 3DEP lidar survey of Brooks Camp in Katmai National Park, Alaska (2012), from the USGS Entwine Point Tile archive on AWS, AK_BrooksCamp_2012. EPT stores a survey as an octree of LAZ files; the setup downloads every node under node 4-3-13-7 — 287 files, 82.8 MB — and joins them into one file, tile.laz. The archive is in Web Mercator (EPSG:3857), so its "metres" are not metres at 58.6° N; see the coordinates guide.

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 survey report before any analysis

import laspy
import numpy as np

CLASSES = {0: "never classified", 1: "unclassified", 2: "ground",
           3: "low vegetation", 4: "medium vegetation", 5: "high vegetation",
           6: "building", 7: "low noise", 9: "water", 11: "road surface",
           17: "bridge deck", 18: "high noise"}


def survey_report(path):
    """Everything you need to know before writing analysis code."""
    las = laspy.read(path)
    n = len(las.points)
    x, y, z = np.asarray(las.x), np.asarray(las.y), np.asarray(las.z)
    area = (x.max() - x.min()) * (y.max() - y.min())

    print(f"  {n:,} points, LAS {las.header.version}, "
          f"point format {las.header.point_format.id}")
    print(f"  extent {x.max() - x.min():.0f} x {y.max() - y.min():.0f} "
          f"(units of the file CRS)")
    print(f"  z {z.min():.1f} to {z.max():.1f}")
    print(f"  nominal density {n / area:.1f} per square unit "
          f"— meaningless unless the CRS is metric and undistorted")

    codes, counts = np.unique(las.classification, return_counts=True)
    for code, count in sorted(zip(codes.tolist(), counts.tolist()),
                              key=lambda t: -t[1]):
        print(f"    {code:2d} {CLASSES.get(code, 'reserved'):18} "
              f"{count:11,}  {count / n:6.2%}")

    returns, rcounts = np.unique(las.return_number, return_counts=True)
    print("  returns: " + ", ".join(f"{r}:{c / n:.1%}"
                                    for r, c in zip(returns, rcounts)))
    multi = np.asarray(las.number_of_returns) > 1
    print(f"  multi-return pulses: {multi.mean():.1%} of points")

    sources = np.unique(las.point_source_id)
    print(f"  {len(sources)} flight line(s)")
    return {"n": n, "classes": dict(zip(codes.tolist(), counts.tolist()))}


report = survey_report("tile.laz")

# Web Mercator stretches both axes by 1/cos(latitude): measure the area on the ground
from pyproj import Transformer

x, y = np.asarray(las.x), np.asarray(las.y)
lon, lat = Transformer.from_crs(3857, 4326, always_xy=True).transform(x.mean(), y.mean())
scale = np.cos(np.radians(lat))
ground_area = (x.max() - x.min()) * scale * (y.max() - y.min()) * scale
print(f"  centre {lat:.3f}°N {lon:.3f}°; {(x.max() - x.min()) * scale:.0f} x "
      f"{(y.max() - y.min()) * scale:.0f} m on the ground, "
      f"{len(x) / ground_area:.1f} points per square metre")

  12,968,770 points, LAS 1.2, point format 1
  extent 1613 x 1613 (units of the file CRS)
  z 15.9 to 203.1
  nominal density 5.0 per square unit — meaningless unless the CRS is metric and undistorted
     1 unclassified        10,003,338  77.13%
     2 ground               2,965,430  22.87%
     7 low noise                    2   0.00%
  returns: 1:59.5%, 2:36.1%, 3:4.4%, 4:0.1%, 5:0.0%
  multi-return pulses: 76.5% of points
  9 flight line(s)
  centre 58.566°N -155.783°; 841 x 841 m on the ground, 18.3 points per square metre

Example 2 — deriving height above ground without vegetation classes

import numpy as np


def height_above_ground(x, y, z, classification, cell=1.0):
    """Normalise heights against a per-cell ground minimum.

    Works when the survey classifies ground and nothing else, which is
    the common case.
    """
    ground = classification == 2
    if not ground.any():
        raise ValueError("no ground points — run a ground filter first")

    left, top = x.min(), y.max()
    width = int(np.ceil((x.max() - left) / cell))
    height = int(np.ceil((top - y.min()) / cell))

    col = np.clip(((x - left) / cell).astype(int), 0, width - 1)
    row = np.clip(((top - y) / cell).astype(int), 0, height - 1)
    flat = row * width + col

    dtm = np.full(width * height, np.inf)
    np.minimum.at(dtm, flat[ground], z[ground])
    dtm[np.isinf(dtm)] = np.nan

    normalised = z - dtm[flat]
    covered = np.isfinite(normalised)
    print(f"  {covered.mean():.1%} of points have a ground cell beneath them")
    print(f"  height above ground: {np.nanmin(normalised):.1f} to "
          f"{np.nanmax(normalised):.1f} m")
    return normalised


normalised = height_above_ground(x, y, np.asarray(las.z), np.asarray(las.classification))

  75.8% of points have a ground cell beneath them
  height above ground: -1.2 to 22.1 m

The 24.2% without a ground cell beneath them are the honest output — and the cells here are one Web Mercator unit, 0.52 m on the ground, so many are simply smaller than the spacing of ground returns. The small negative heights are unclassified points below the lowest ground point in their cell. Filling them by interpolating the DTM is reasonable and should be recorded, because those heights rest on an interpolated ground rather than an observed one.

Example 3 — separating flight lines

import numpy as np


def flight_line_report(point_source_id, gps_time, intensity):
    """Overlapping lines double the density and shift the intensity."""
    lines = np.unique(point_source_id)
    print(f"  {len(lines)} flight lines")
    for line in lines[:10]:
        m = point_source_id == line
        span = (gps_time[m].max() - gps_time[m].min()) / 60
        print(f"    line {line:5d}: {m.sum():10,} points, "
              f"{span:6.1f} min, median intensity {np.median(intensity[m]):6.0f}")

    if len(lines) > 1:
        medians = [np.median(intensity[point_source_id == l]) for l in lines]
        spread = (max(medians) - min(medians)) / np.mean(medians)
        print(f"  intensity varies {spread:.0%} between lines "
              f"— do not compare intensity across them")
    return lines


lines = flight_line_report(np.asarray(las.point_source_id), np.asarray(las.gps_time),
                           np.asarray(las.intensity))

  9 flight lines
    line     5:    946,235 points,    0.2 min, median intensity    127
    line     7:  3,003,561 points,    0.3 min, median intensity    108
    line     8:    844,661 points,    0.2 min, median intensity    118
    line    10:  2,918,818 points,    0.3 min, median intensity    108
    line    14:    150,097 points,    0.1 min, median intensity    130
    line    15:    671,206 points,    0.2 min, median intensity    139
    line    18:    222,052 points,    0.1 min, median intensity    127
    line    25:  2,320,279 points,    0.3 min, median intensity    113
    line    26:  1,891,861 points,    0.2 min, median intensity    119
  intensity varies 26% between lines — do not compare intensity across them

An intensity raster built across several flight lines shows the flight pattern rather than the ground. Rasterise intensity per line, or normalise per line, or do not use intensity.

Density in the wrong CRS, and classification as a product

Why point density is the wrong number in the wrong CRS

The naive density of this survey is 5.0 points per square metre — and the true density is 18.3.

The file is stored in EPSG:3857, Web Mercator, where distances are inflated by 1/cos(latitude). At 58.57° north that is a factor of 1.92, so areas are inflated by 3.68, and any density computed in those units is understated by the same factor.

The measured extent makes it concrete: 1,613 × 1,613 Web Mercator "metres" for a patch that is 841 × 841 m on the ground. Reprojecting helps only if you then measure the area properly: in UTM zone 5N the same points have a bounding box of 877 × 876 m, because the UTM grid is turned 2.4° against Web Mercator's there, and dividing by that box gives 16.9 — the convex hull gives 18.3.

Every metric quantity — density, slope, canopy area, distance between points — inherits that error. Reproject to a projected CRS suited to the area before computing anything metric. See laspy returns the wrong coordinates.

Why classification is a separate product

The ground/not-ground decision is a filtering problem, not a property of the measurement. Algorithms such as progressive TIN densification or cloth simulation start from the lowest returns and iteratively decide which points are consistent with a plausible ground surface.

They can fail: on steep slopes they clip ridges, under dense canopy they have little to work with, and near buildings they can accept a flat roof as ground. Vendors deliver ground classification because it is hard and worth paying for — and they often deliver only that, because the other classes are harder still.

The 100% / 38.4% relationship measured above is the reason. Filtering by last return gets you a superset of the ground that is nearly three times too large.

Why intensity looks like a band and is not one

Returned amplitude depends on the target's reflectance at the laser wavelength — usually 1,064 nm, in the near-infrared, so vegetation is bright and water is nearly black.

It also depends on the range to target (inverse square), the incidence angle, the atmospheric transmission on the day, and any automatic gain the sensor applied. None of that is recorded in a way that lets you invert it.

The practical rule: intensity is a texture, useful for eyeballing and for local contrast within one flight line. It is not a reflectance measurement and it is not comparable across lines, dates or sensors.

Why the file is a table and behaves like one

A LAS file is a header plus fixed-width records. The coordinates are stored as 32-bit integers with a scale and offset in the header — hence X and x being different attributes in laspy, one raw and one scaled.

That layout is why point clouds are fast to read sequentially and slow to query spatially: there is no index. Formats that add one — COPC, EPT — exist precisely to make "give me this bounding box at this resolution" a cheap operation. See LAS, LAZ and COPC explained.

Unclassified points and raw integers

• Class 1 means "unclassified", not "unknown surface". A cloud that is 77% class 1 has only had a ground filter run.
• las.X is the raw integer, las.x is the scaled coordinate. Mixing them is off by a factor of 100 or 1,000.
• Every ground point is a last return; only 38.4% of last returns are ground in this survey.
• Intensity is uncalibrated. Do not compare it across flight lines.
• Check the CRS before computing density. Web Mercator understated it by a factor of 3.68 here, and a bounding box in a rotated grid by another 8%.
• number_of_returns can be 0 in malformed files; guard division by it.
• GPS time can be either "week seconds" or "adjusted standard time" — the header's global encoding bit says which.
• Withheld and synthetic flags mark points the vendor wants excluded; check them before analysis.

Internal links

• LAS, LAZ and COPC explained — how the file is laid out and why queries are slow
• Point density explained: what resolution a point cloud supports — the density measurement done properly
• How to read a LAS or LAZ file in Python with laspy — the practical reading guide
• How to filter a point cloud by class and return number — using these columns
• DSM, DTM and CHM from LiDAR: how each surface is derived — turning returns into rasters
• laspy returns the wrong coordinates: scale, offset and CRS — the Web Mercator trap
• Digital elevation models explained: DEM, DSM and DTM — the raster products
• How to create a DTM from LiDAR ground points in Python — the first real product

FAQ

What is a lidar return?

One detectable echo from one outgoing laser pulse. A pulse through vegetation can produce several returns, which is why return_number and number_of_returns are separate columns.

What do the classification codes mean?

The LAS standard defines 2 as ground, 3–5 as vegetation by height, 6 as building, 9 as water and 7/18 as noise. A survey only contains the classes it was paid to produce — this one had ground and nothing else.

Can I use last returns instead of ground classification?

Only as a rough filter. Every ground point here was a last return, but only 38.4% of last returns were ground.

Is lidar intensity like a satellite band?

No. It is uncalibrated and varies with range, incidence angle, atmosphere and sensor gain, so it is not comparable between flight lines.

How dense is my point cloud?

Points divided by ground area — but only after reprojecting to a suitable projected CRS. In Web Mercator this survey appeared to be 5.0 points/m² when it is really 18.3.

Why does my cloud have no vegetation points?

Because the survey did not classify vegetation. Derive height above ground from the returns and the ground surface instead of filtering by class.

What is point_source_id for?

It identifies the flight line. Overlapping lines double the local density and shift the intensity, which shows up as stripes in derived rasters.