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.47%, single 23.48%

Step-by-step solution #

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:

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%

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:

DTM from last returns vs from classified ground
  holes    7.48% vs 10.93%
  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 4 cm downward bias comes from low returns the classifier had excluded.

Code examples #

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

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

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

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.

Explanation #

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.

Edge cases or notes #

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