🧭 Concept LiDAR & point clouds LiDAR & Point Clouds · step 1 of 14 Updated

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 877 x 876 m
16.9 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%})")
180,756 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.

One laser pulse producing several returns as it passes through canopy, with the first return at the treetop and the last on the ground.
One pulse, several returns. Which return you keep decides whether you are mapping the canopy or the ground.

Step-by-step solution

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:

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

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.

All ground points being last returns while only 38 percent of last returns are ground, shown as nested sets.
Last return is necessary for ground and nowhere near sufficient. That gap is what a ground classifier exists to close.

Code examples

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

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
  89.1% of points have a ground cell beneath them
  height above ground: 0.0 to 22.7 m

The 10.9% without a ground cell beneath them are the honest output. 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

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.

Explanation

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

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 877 × 876 m on the ground.

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.

A real survey classified 77 percent unclassified, 23 percent ground, and two noise points, with no vegetation or building classes.
The LAS standard defines twelve classes. A survey contains the ones it was paid to produce.

Edge cases or notes

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

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

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.