How to Filter a Point Cloud by Class and Return Number

Problem statement

Almost every point cloud operation starts with a selection, and the two columns that drive it — classification and return_number — are related in ways that are easy to get backwards.

Measured on a real 3DEP survey of 12,968,770 points:

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

Every ground point is a last return, and fewer than two in five last returns are ground. Filtering by "last return" gives you a set two and a half times too large.

And the classes you expect may not exist. This survey is 77.1% class 1 (unclassified) and 22.9% class 2 (ground), with two noise points. There are no vegetation or building classes at all.

Quick answer

Check what is there before filtering:

import numpy as np

codes, counts = np.unique(classification, return_counts=True)
for code, count in zip(codes, counts):
    print(f"  class {code:2d}: {count:11,}  {count / len(classification):6.2%}")

ground = classification == 2
first = return_number == 1
last = return_number == number_of_returns
single = number_of_returns == 1
intermediate = (return_number > 1) & (return_number < number_of_returns)

print(f"  ground {ground.mean():.2%}, first {first.mean():.2%}, "
      f"last {last.mean():.2%}, single {single.mean():.2%}")

  class  1:  10,003,338  77.13%
  class  2:   2,965,430  22.87%
  class  7:           2   0.00%
  ground 22.87%, first 59.47%, last 59.54%, single 23.48%

The class histogram, return predicates, the right filter

1. Print the class histogram, always

The LAS standard defines classes 0–18. It does not require a survey to populate them. Vendors classify what they were paid to classify, and ground-only deliveries are common.

Code that filters classification == 5 on this survey returns an empty array and everything downstream fails with a confusing error about zero-length input.

2. Get the return predicates right

first        = return_number == 1
last         = return_number == number_of_returns
single       = number_of_returns == 1
intermediate = (return_number > 1) & (return_number < number_of_returns)

Two things people get wrong. last is not return_number == number_of_returns.max() — the maximum is a property of the file, not the pulse. And a single return is both first and last, which is correct and occasionally surprising in a Venn diagram.

Measured distribution:

values, counts = np.unique(return_number, return_counts=True)
for value, count in zip(values, counts):
    print(f"return {value}: {count:10,}  {count / len(return_number):6.2%}")
print(f"multi-return pulses: {(number_of_returns > 1).mean():.1%} of points")

return 1:  7,712,998  59.47%
return 2:  4,675,883  36.05%
return 3:    568,150   4.38%
return 4:     11,679   0.09%
return 5:         60   0.00%
multi-return pulses: 76.5% of points

76.5% of points came from a pulse that produced more than one return — a high figure, indicating well-vegetated ground and a sensor that resolves closely spaced echoes.

3. Choose the filter from what the product needs

4. Exclude noise explicitly

Classes 7 (low noise) and 18 (high noise) exist to be dropped. This survey has two class-7 points — the easy case. Unclassified low points are the ones that ruin a minimum-based DTM, and no class flags them.

usable = ~np.isin(classification, [7, 18])

5. Use last returns as a ground proxy only knowingly

If the survey has no ground class, last returns are the available approximation. Measured against the classified ground on this survey, with a 1 m lowest-point DTM and only the cells inside the survey's footprint counted:

import shapely
from rasterio.transform import from_origin
from scipy.spatial import ConvexHull

CELL = 1.0
left, top = x.min(), y.max()
width, height = int(np.ceil((x.max() - left) / CELL)), int(np.ceil((top - y.min()) / CELL))
transform = from_origin(left, top, CELL, CELL)
cols = np.clip(((x - left) / CELL).astype(int), 0, width - 1)
rows = np.clip(((top - y) / CELL).astype(int), 0, height - 1)


def lowest(mask):
    grid = np.full(height * width, np.inf)
    np.minimum.at(grid, rows[mask] * width + cols[mask], z[mask])
    grid[np.isinf(grid)] = np.nan
    return grid.reshape(height, width)


hull = ConvexHull(np.c_[x, y])
cx, cy = np.meshgrid(left + (np.arange(width) + 0.5) * CELL, top - (np.arange(height) + 0.5) * CELL)
surveyed = shapely.contains_xy(shapely.Polygon(np.c_[x, y][hull.vertices]), cx, cy)

dtm = lowest(classification == 2)
from_last = lowest((return_number == number_of_returns) & ~np.isin(classification, [7, 18]))
diff = (from_last - dtm)[surveyed & np.isfinite(from_last) & np.isfinite(dtm)]
print("DTM from last returns vs from classified ground")
print(f"  holes    {np.isnan(from_last[surveyed]).mean():.2%} vs {np.isnan(dtm[surveyed]).mean():.2%}")
print(f"  bias    {diff.mean():+.3f} m mean, p1 {np.percentile(diff, 1):+.2f} m, "
      f"p99 {np.percentile(diff, 99):+.2f} m")

DTM from last returns vs from classified ground
  holes    0.00% vs 3.67%
  bias    -0.040 m mean, p1 -0.33 m, p99 +0.00 m

Fewer holes, and never higher — because a per-cell minimum over a superset can only fall. The small downward bias comes from low returns the classifier had excluded from ground.

A filter builder that fails loudly

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: every node under octree node 4-3-13-7 — 287 LAZ files, 12,968,770 points — joined into tile.laz, then reprojected from the archive's Web Mercator to UTM zone 5N and loaded as plain arrays.

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")

import numpy as np
from pyproj import Transformer

las = laspy.read("tile.laz")
x, y = Transformer.from_crs(3857, 32605, always_xy=True).transform(las.x, las.y)   # UTM 5N metres
z = np.asarray(las.z)
classification = np.asarray(las.classification)
return_number = np.asarray(las.return_number)
number_of_returns = np.asarray(las.number_of_returns)

287 nodes, 12,968,770 points, 82.8 MB of nodes

Example 1 — a filter builder that fails loudly

import numpy as np

CLASS_NAMES = {0: "never classified", 1: "unclassified", 2: "ground",
               3: "low vegetation", 4: "medium vegetation",
               5: "high vegetation", 6: "building", 7: "low noise",
               9: "water", 11: "snow", 18: "high noise"}
NOISE = (7, 18)


def build_filter(classification, return_number, number_of_returns,
                 classes=None, returns=None, drop_noise=True,
                 min_points=1000):
    """A boolean mask, with a clear error when a requested class is absent."""
    present = set(np.unique(classification).tolist())

    if classes is not None:
        missing = set(classes) - present
        if missing:
            names = ", ".join(f"{c} ({CLASS_NAMES.get(c, '?')})"
                              for c in sorted(missing))
            raise ValueError(
                f"requested class(es) {names} are not in this file. "
                f"Present: {sorted(present)}. Derive the property from the "
                "returns instead of relying on a classification."
            )
        mask = np.isin(classification, list(classes))
    else:
        mask = np.ones(len(classification), bool)

    if returns == "first":
        mask &= return_number == 1
    elif returns == "last":
        mask &= return_number == number_of_returns
    elif returns == "single":
        mask &= number_of_returns == 1
    elif returns == "intermediate":
        mask &= (return_number > 1) & (return_number < number_of_returns)
    elif returns is not None:
        raise ValueError(f"unknown return filter {returns!r}")

    if drop_noise:
        mask &= ~np.isin(classification, NOISE)

    print(f"  {int(mask.sum()):,} of {len(mask):,} points "
          f"({mask.mean():.2%}) after classes={classes} returns={returns}")
    if mask.sum() < min_points:
        raise ValueError(f"filter left only {int(mask.sum()):,} points — "
                         "check the class histogram before proceeding")
    return mask


ground_last = build_filter(classification, return_number, number_of_returns,
                           classes=[2], returns="last")
try:
    build_filter(classification, return_number, number_of_returns, classes=[5])
except ValueError as err:
    print(f"  ValueError: {err}")

  2,965,430 of 12,968,770 points (22.87%) after classes=[2] returns=last
  ValueError: requested class(es) 5 (high vegetation) are not in this file. Present: [1, 2, 7]. Derive the property from the returns instead of relying on a classification.

The explicit error on a missing class is worth more than any other line here. Silently returning an empty array is how a batch job produces a hundred empty rasters.

Example 2 — the class-versus-return cross-tabulation

import numpy as np
import pandas as pd


def class_return_table(classification, return_number, number_of_returns):
    """Which classes appear in which returns? Usually revealing."""
    kind = np.full(len(classification), "intermediate", dtype=object)
    kind[number_of_returns == 1] = "single"
    kind[(return_number == 1) & (number_of_returns > 1)] = "first of many"
    kind[(return_number == number_of_returns) &
         (number_of_returns > 1)] = "last of many"

    table = pd.crosstab(pd.Series(classification, name="class"),
                        pd.Series(kind, name="return"), normalize="index")
    print(table.round(3).to_string())
    return table


table = class_return_table(classification, return_number, number_of_returns)

return  first of many  intermediate  last of many  single
class                                                    
1               0.467         0.058         0.270   0.206
2               0.000         0.000         0.667   0.333
7               0.000         0.000         0.500   0.500

The ground row is the useful one: class 2 is 66.7% last-of-many and 33.3% single, and zero first-of-many or intermediate. That is exactly what physics predicts and a good check that the classification is sane.

A ground class with first-of-many returns in it would mean the classifier had accepted canopy tops as ground.

Example 3 — deriving height classes when the survey has none

import numpy as np


def height_classes(x, y, z, classification, dtm, transform, cell,
                   breaks=(0.5, 2.0, 5.0)):
    """Assign vegetation classes from height above ground."""
    left, top = transform.c, transform.f
    height, width = dtm.shape
    col = np.clip(((x - left) / cell).astype(int), 0, width - 1)
    row = np.clip(((top - y) / cell).astype(int), 0, height - 1)
    above = z - dtm[row, col]

    derived = np.array(classification, copy=True)
    unclassified = classification == 1
    ok = unclassified & np.isfinite(above)

    derived[ok & (above < breaks[0])] = 2                    # ground-ish
    derived[ok & (above >= breaks[0]) & (above < breaks[1])] = 3
    derived[ok & (above >= breaks[1]) & (above < breaks[2])] = 4
    derived[ok & (above >= breaks[2])] = 5

    codes, counts = np.unique(derived[ok], return_counts=True)
    for code, count in zip(codes.tolist(), counts.tolist()):
        print(f"    -> {code} {CLASS_NAMES.get(code, '?'):18} "
              f"{count:11,}  {count / ok.sum():6.2%}")
    print(f"  {int((~np.isfinite(above)).sum()):,} points had no ground "
          f"beneath them and keep their original class")
    return derived


derived = height_classes(x, y, z, classification, dtm, transform, CELL)

    -> 2 ground               3,180,076  33.11%
    -> 3 low vegetation         336,872   3.51%
    -> 4 medium vegetation    1,225,922  12.77%
    -> 5 high vegetation      4,860,813  50.61%
  399,655 points had no ground beneath them and keep their original class

This is height stratification, not classification: it says nothing about whether a point is vegetation or a wall. It is enough for canopy metrics and not enough for anything that needs to distinguish a tree from a building.

Every ground point is a last return

Why every ground point is a last return

A pulse stops when it hits something opaque. Ground is opaque, so nothing beyond it returns, so a ground echo is necessarily the final one from that pulse.

Measured: 100.0% of class-2 points satisfy return_number == number_of_returns. That is not a coincidence, it is geometry, and it is a useful invariant to assert as a data-quality check.

Why so few last returns are ground

The converse fails badly — 38.4% here — because a pulse can stop on anything opaque. A roof, a dense branch cluster, a wall and a car all produce a last return well above the ground.

That gap is the whole reason ground classification exists as a separate, difficult step. Filtering by last return gives you a superset that includes every hard surface in the scene.

Under closed canopy the gap widens: many pulses never reach the ground at all, so their last return is a branch, and no filter on return number can tell.

Why "unclassified" does not mean "not ground"

Class 1 is "unclassified" — the classifier looked and did not assign a class, or was never asked. It is not a statement that the point is vegetation.

In a ground-only delivery, class 1 contains everything that is not ground: canopy, buildings, wires, vehicles and any ground point the filter missed. Treating class 1 as vegetation will include roofs, and treating it as noise will discard the entire above-ground scene.

Why single returns matter

A single return means the pulse hit something opaque with nothing in front of it — open ground, a roof, water, a road.

Measured, 23.5% of points overall are single returns and 33.3% of ground points are. The other two thirds of ground points are the last of several returns — the pulse passed through vegetation and still reached the ground.

That ratio is a useful summary of how obstructed the site is. Over genuinely open terrain the ground would be mostly single returns; here two thirds of ground echoes came through something first.

Last returns and single-return pulses

• Print the class histogram before filtering. Many surveys classify only ground.
• last is return_number == number_of_returns, per pulse — not the file maximum.
• A single return is both first and last.
• Class 1 is "unclassified", not "vegetation".
• Drop classes 7 and 18 before any minimum-based reduction.
• Ground points must all be last returns. If they are not, the classification is suspect.
• Point formats below 6 cap returns at 5 and classes at 31.
• Check withheld and synthetic flags; they mark points the vendor excluded.

Internal links

• LiDAR point clouds explained: returns, classes and intensity — what these columns mean
• How to read a LAS or LAZ file in Python with laspy — getting the arrays
• How to create a DTM from LiDAR ground points in Python — the main consumer of a class filter
• DSM, DTM and CHM from LiDAR: how each surface is derived — which filter each surface needs
• How to create a canopy height model from LiDAR — using normalised heights instead of classes
• Point density explained — density per class
• My LiDAR DTM has holes, spikes or terraces — what noise classes cause
• How to measure building heights from LiDAR in Python — when the building class is missing

FAQ

How do I select ground points from a point cloud?

classification == 2, after checking the class histogram — many surveys classify only ground, and some classify nothing.

Can I use last returns instead of ground classification?

Only knowingly. Every ground point is a last return, but only 38.4% of last returns are ground, so you get a superset that includes roofs and dense canopy.

What does classification 1 mean?

Unclassified — the classifier did not assign a class. It is not a statement that the point is vegetation, and in a ground-only delivery it contains the entire above-ground scene.

How do I select first and last returns?

return_number == 1 for first, return_number == number_of_returns for last. A single return satisfies both.

Why does my class filter return nothing?

The class is not in the file. Print np.unique(classification) first; vegetation and building classes are frequently absent.

Should I drop noise classes?

Yes, before any minimum- or maximum-based reduction. Classes 7 and 18 exist for exactly that purpose.

What if I need vegetation classes and the survey has none?

Derive height above ground from a DTM and stratify by height. That is height stratification, not classification — it cannot tell a tree from a wall.