Python LiDAR

Building a DSM from First Returns

TL;DR: Add {"type": "filters.range", "limits": "ReturnNumber[1:1]"} to keep only first returns, then rasterize with writers.gdal using output_type=max at a resolution matched to your point spacing — the maximum first-return elevation in each cell becomes the Digital Surface Model, which you confirm by opening the GeoTIFF in rasterio and checking the CRS, elevation range, and void fraction.

# Context and Motivation

This guide is part of DSM Generation from LiDAR with PDAL, and it drills into one specific, high-value technique: constraining the rasterizer to the leading pulse return so the resulting surface is unambiguously the top of the world. Where the parent page surveys the whole surface-modelling workflow, this recipe is the copy-and-run version you reach for when a client needs a canopy-and-structure raster by end of day.

The reason first returns deserve their own page is that they are the physically cleanest definition of a surface. When an airborne pulse descends into a forest, the first return fires off the outermost leaf it touches; subsequent returns leak through gaps and report twigs, branches, and eventually the ground. If you rasterize every return with a maximum statistic you usually get the right answer, but a single anomalously high second return from an adjacent pulse can poke a spike through your surface. Filtering to ReturnNumber == 1 first removes that whole failure mode before the grid is ever populated. This recipe uses EPSG:32611 (WGS84 / UTM Zone 11N) and a 0.4 m grid to keep the examples concrete and distinct from the regional Albers example on the parent page.

Selecting the first return from a multi-return pulse A vertical pulse enters a tree canopy and produces four returns at descending heights: return 1 at the canopy top, returns 2 and 3 in the mid-canopy, and return 4 at the ground. A callout marks return 1 as the point kept by the ReturnNumber filter and fed into the DSM. incoming pulse ground ReturnNumber = 1 2 3 4 (last, ground) filters.range keeps this point → writers.gdal output_type = max

# Prerequisites and Assumptions

Requirement Detail
PDAL 2.5+ with Python bindings (pip install pdal)
rasterio 1.3+ for output verification
Input file LAS/LAZ with populated ReturnNumber and NumberOfReturns
CRS Metre-based projected CRS; this recipe uses EPSG:32611
Point density Known average spacing — see Point Density Metrics

Before writing any pipeline, confirm the tile actually carries multi-return data. A file where every point reports ReturnNumber == 1 came from a single-return sensor, and the filter below will be a harmless no-op:

bash
pdal info --stats --dimensions ReturnNumber,NumberOfReturns first_return_tile.laz

Look for NumberOfReturns maxima above 1. If the maximum is 1 everywhere, skip straight to rasterization — there is no understory to exclude.

# Step-by-Step Implementation

# Step 1 — Keep only the first return

filters.range selects points by a numeric interval on any dimension. The interval ReturnNumber[1:1] means “keep points where ReturnNumber is between 1 and 1 inclusive” — i.e. exactly the first return. This shares the interval syntax used across Pipeline Filtering Logic.

json
{
  "type": "filters.range",
  "limits": "ReturnNumber[1:1]"
}

# Step 2 — Rasterize the maximum elevation

writers.gdal bins the surviving points onto a grid and writes one statistic per cell. output_type set to max records the highest first-return Z in each cell — the top surface.

json
{
  "type": "writers.gdal",
  "filename": "first_return_dsm.tif",
  "gdaldriver": "GTiff",
  "output_type": "max",
  "resolution": 0.4,
  "nodata": -9999.0,
  "override_srs": "EPSG:32611",
  "data_type": "float32"
}

# Step 3 — Assemble the full pipeline

Wired together, the reader, range filter, and GDAL writer form a three-stage pipeline. This is the same reader-filter-writer chaining pattern described in PDAL Stage Chaining.

json
{
  "pipeline": [
    {"type": "readers.las", "filename": "first_return_tile.laz"},
    {"type": "filters.range", "limits": "ReturnNumber[1:1]"},
    {
      "type": "writers.gdal",
      "filename": "first_return_dsm.tif",
      "gdaldriver": "GTiff",
      "output_type": "max",
      "resolution": 0.4,
      "nodata": -9999.0,
      "override_srs": "EPSG:32611",
      "data_type": "float32"
    }
  ]
}

# Step 4 — Verify the raster

Never ship a surface raster you have not opened. rasterio reads the CRS, cell size, elevation range, and void mask in a few lines — enough to catch every common failure.

# Complete Working Example

Save this as first_return_dsm.py and run it against any multi-return LAZ tile. It builds the pipeline in Python, executes it, then verifies the output GeoTIFF and raises loudly if anything looks wrong.

python
#!/usr/bin/env python3
"""
first_return_dsm.py
Build a Digital Surface Model from first-return LiDAR points with PDAL,
then verify the output raster with rasterio.

Usage:
    python first_return_dsm.py input_tile.laz output_dsm.tif 0.4 32611

Requirements:
    pip install pdal rasterio
    PDAL 2.5+ (conda install -c conda-forge pdal python-pdal)
"""

import json
import sys

import pdal
import rasterio


def build_first_return_dsm(
    input_path: str,
    output_path: str,
    resolution: float,
    epsg: int,
    nodata: float = -9999.0,
) -> int:
    """
    Rasterize the first-return top surface of a LiDAR tile into a DSM GeoTIFF.

    Args:
        input_path:  source LAS/LAZ tile with multi-return data.
        output_path: destination GeoTIFF for the surface raster.
        resolution:  cell size in CRS units (metres).
        epsg:        projected EPSG code stamped into the raster.
        nodata:      sentinel written to cells with no first return.

    Returns:
        Number of points rasterized (first returns only).
    """
    pipeline_dict = {
        "pipeline": [
            {"type": "readers.las", "filename": input_path},
            # Keep only the leading pulse return — the visible top surface.
            {"type": "filters.range", "limits": "ReturnNumber[1:1]"},
            {
                "type": "writers.gdal",
                "filename": output_path,
                "gdaldriver": "GTiff",
                "output_type": "max",
                "resolution": resolution,
                "nodata": nodata,
                "override_srs": f"EPSG:{epsg}",
                "data_type": "float32",
            },
        ]
    }

    pipeline = pdal.Pipeline(json.dumps(pipeline_dict))
    count = pipeline.execute()

    if count == 0:
        raise RuntimeError(
            f"No first returns rasterized from '{input_path}'. "
            "Check that ReturnNumber is populated and non-zero."
        )
    return count


def verify_dsm(output_path: str, nodata: float = -9999.0) -> None:
    """Assert the DSM has a CRS, a plausible elevation range, and low voids."""
    with rasterio.open(output_path) as src:
        band = src.read(1, masked=True)
        void_fraction = float(band.mask.mean()) if band.mask.ndim else 0.0

        assert src.crs is not None, "Output has no CRS — check override_srs"
        assert band.min() > nodata, "nodata sentinel leaking into valid elevations"

        print(f"CRS:            {src.crs}")
        print(f"Grid:           {src.width} x {src.height} @ {src.res[0]} m")
        print(f"Elevation:      {band.min():.2f} to {band.max():.2f} m")
        print(f"Void fraction:  {void_fraction:.1%}")

        if void_fraction > 0.10:
            print(
                "WARNING: void fraction above 10% — resolution may be finer "
                "than the first-return spacing; try coarsening the grid."
            )


def main() -> None:
    if len(sys.argv) != 5:
        print("Usage: python first_return_dsm.py <in.laz> <out.tif> <resolution> <epsg>")
        sys.exit(1)

    input_path = sys.argv[1]
    output_path = sys.argv[2]
    resolution = float(sys.argv[3])
    epsg = int(sys.argv[4])

    print(f"Building DSM from first returns: {input_path} -> {output_path}")
    n = build_first_return_dsm(input_path, output_path, resolution, epsg)
    print(f"Rasterized {n:,} first-return points.")
    verify_dsm(output_path)


if __name__ == "__main__":
    main()

# Key Parameter Table

Parameter Stage Type Value here Notes
limits filters.range string ReturnNumber[1:1] Keeps first returns only; the interval is inclusive on both ends
output_type writers.gdal string max Highest Z per cell — the top-surface definition of a DSM
resolution writers.gdal float 0.4 Cell size in metres; coarsen if voids climb above ~10%
nodata writers.gdal float -9999.0 Empty-cell sentinel; keep outside the real elevation range
override_srs writers.gdal string EPSG:32611 Stamps the CRS into the GeoTIFF
data_type writers.gdal string float32 Preserves centimetre-level elevation precision

# Verification

Beyond the assertions baked into the script, three checks confirm the surface is trustworthy.

Confirm the point count dropped as expected. On multi-return vegetated data, filtering to first returns typically retains 55–75% of points. If the count barely changed, the data is largely single-return; if it collapsed to a tiny fraction, ReturnNumber may be mis-encoded.

python
import pdal, json

def count(path, first_only):
    stages = [{"type": "readers.las", "filename": path}]
    if first_only:
        stages.append({"type": "filters.range", "limits": "ReturnNumber[1:1]"})
    p = pdal.Pipeline(json.dumps({"pipeline": stages}))
    return p.execute()

raw = count("first_return_tile.laz", False)
first = count("first_return_tile.laz", True)
print(f"first returns: {first:,} / {raw:,}  ({first/raw:.0%})")

Eyeball the elevation range. z_max should equal the tallest feature in the scene. A DSM over a forested UTM tile might read 612 m of bare valley floor up to 648 m of canopy crown; a max that matches the ground minimum means first returns never reached anything tall.

Render a quick hillshade. Loading the DSM into a hillshade — see Hillshade, Slope and Aspect — instantly reveals whether buildings and tree crowns cast crisp shadows (good) or the surface is mushy and full of holes (resolution too fine).

# Gotchas and Edge Cases

1. ReturnNumber is zero-based or unpopulated. A handful of exporters write ReturnNumber starting at 0, or leave it at 0 entirely. The interval [1:1] then keeps nothing and the pipeline rasterizes zero points. Run pdal info --stats on ReturnNumber first; if the minimum is 0 and the maximum is 0, the field is empty — rasterize all points instead.

2. Overlapping flight lines double-count first returns. Where two flight lines overlap, a cell may receive first returns from both passes. max still picks the highest, so the surface stays correct, but the point density in the overlap is inflated, which can mask a resolution that is otherwise too fine elsewhere. Judge void fraction from a non-overlap region.

3. Birds, dust, and low noise become spikes. A first return off a bird sits well above the canopy and will win the max for its cell. If your hillshade shows isolated pinnacle spikes, insert an outlier removal stage before the range filter — the statistical outlier filter handles exactly this.

4. The CRS is geographic, not projected. If the tile is in EPSG:4326, resolution: 0.4 means 0.4 degrees — tens of kilometres per cell. Reproject to a metre-based CRS first with Spatial Reprojection, or the raster will be a single meaningless pixel.

# Frequently Asked Questions

Why filter to ReturnNumber==1 instead of using all points?

First returns are the outermost surface the pulse struck — canopy top, roof edge, wire. Keeping only ReturnNumber == 1 drops understory and last-return ground hits before rasterization, so the max statistic reflects the true visible skin of the scene rather than being pulled by stray high points from other pulses in the same cell. On clean data the difference is small, but it eliminates a whole class of spike artifacts for free.

What if my LiDAR only has single returns?

Single-return sensors record one point per pulse, so every point has ReturnNumber 1 and NumberOfReturns 1. The filter is then a no-op and the DSM is built from all points — which is correct, because there is no understory signal to exclude. The recipe still works unchanged.

How do I choose the resolution for a first-return DSM?

Filtering to first returns roughly halves the point count on vegetated scenes, so the effective spacing between surviving points is coarser than the raw density suggests. Start at the raw average spacing and coarsen one step if the void fraction exceeds about 10 percent. The parent DSM Generation page tabulates the resolution-versus-void tradeoff.

Can I keep last returns instead to approximate the ground?

You can filter ReturnNumber against NumberOfReturns to keep last returns, but last returns are not the same as classified ground — they include understory and building interiors. For a true bare-earth surface use ground classification and DTM Raster Generation rather than a last-return proxy.

Related in DSM Generation from LiDAR with PDAL

↑ DSM Generation from LiDAR with PDAL overview --- title: "DSM Generation from LiDAR with PDAL" description: "Building a Digital Surface Model from LiDAR first returns with PDAL... DTM vs DSM: Which Surface Model to Generate --- title: "DTM vs DSM: Which Surface Model to Generate" description: "When to build a bare-earth DTM versus a first-return DSM from LiDAR —...