Filter and Classify a Point Cloud in PyQGIS

Almost every LiDAR task begins by throwing points away. A terrain model wants ground returns, a canopy study wants high vegetation, a building extraction wants class 6, and every one of them wants the noise gone first. The mechanics are simple; what trips people up is that QGIS gives you two entirely separate filtering mechanisms with two different expression languages, and that a great many deliveries are not classified at all.

This recipe belongs to Point Cloud & LiDAR Workflows in PyQGIS. It covers choosing between a subset string and a written file, the noise removal that should come first, checking what classes actually exist, and the honest options when the answer is "none of them".

A view of the cloud, or a smaller cloudSetting a subset string leaves the file untouched and changes only what the layer renders and what layer-based queries see; it is instant and reversible. Running pdal:filter reads the whole file and writes a new one containing only the matching points, which costs time and disk but produces something every other tool can open.Filter the view, or filter the datathe cloud on disksetSubsetStringinstant, reversibleaffects display and layer readsinvisible to pdal: algorithmspdal:filterreads and rewrites every pointproduces a real, smaller filewhat you archive and hand onexplore with the left one, then commit to the right one once

Prerequisites

  • QGIS 3.34 LTR or newer with the PDAL provider.
  • An indexed cloud — see loading a point cloud layer.
  • Somewhere to write output, if you are producing a filtered file rather than a view.

Check what is actually in the file first

Writing a filter before checking the classes is how you end up with an empty output and no idea why.

from qgis.core import QgsPointCloudLayer

cloud = QgsPointCloudLayer("/data/indexed/tile_0345.copc.laz", "tile", "copc")
stats = cloud.statistics()

present = sorted(stats.classesOf("Classification"))
print("classes present:", present)
print("points:", cloud.pointCount())
print("Z range:", stats.minimum("Z"), "to", stats.maximum("Z"))

Breakdown: classesOf() returns the distinct values of a categorical attribute as they were found during indexing, so it is a cheap read rather than a scan. Three results are common and each means something different: [2, 5, 6, 9] is a properly classified survey; [1] alone means the delivery was never classified and no class filter will ever work; [1, 2] means somebody ran a ground filter and nothing else, which is enough for a terrain model and useless for anything about buildings. The Z range is a second sanity check — a range of thousands of metres over flat ground means noise points are present.

Remove noise before anything else

Class 7 (low noise) and class 18 (high noise) exist because LiDAR sees birds, aircraft, cloud, multipath reflections and sensor artefacts. A single class-18 point at 4,000 m does not sound like much until it stretches the Z ramp of every rendering and pushes a surface model's maximum by four kilometres.

import processing

clean = processing.run("pdal:filter", {
    "INPUT": "/data/indexed/tile_0345.copc.laz",
    "FILTER_EXPRESSION": "Classification != 7 && Classification != 18",
    "OUTPUT": "/data/clean/tile_0345.copc.laz",
})["OUTPUT"]

Breakdown: PDAL expression syntax again — != and &&, not <> and AND. Writing to a .copc.laz extension keeps the output indexed, so the filtered file is immediately usable without another conversion step; writing to .laz would give you an unindexed file and a re-indexing cost on every later open. If the survey has no noise classes at all, this step still costs a full read and write, so gate it on 7 in present or 18 in present.

Where noise was never classified, a statistical outlier filter is the alternative, and QGIS does not wrap one — that is a case for a PDAL pipeline run outside QGIS, or for clamping with a Z range you can defend from the survey's own metadata.

Filter to the returns you need

Four filters, four different productsA terrain filter keeps class two. A surface filter keeps everything except noise and takes the maximum per cell. A canopy filter keeps classes three to five. A building filter keeps class six. Each one produces a completely different picture of the same place, and choosing wrongly is the most common cause of a plausible-looking but incorrect result.The filter is the analysisterrainClassification == 2surfaceClassification != 7canopyClassification >= 3 && <= 5buildingsClassification == 6the same tile, four ways — and only one of them answers your questionflood modelling wants terrain; solar potential wants surface; forestry wants canopy minus terraina surface model used as a terrain model puts every roof into the catchment

For exploring, set a subset string and look at the canvas. It is instant and it undoes cleanly.

cloud.setSubsetString("Classification IN (3, 4, 5)")
cloud.triggerRepaint()

print(cloud.subsetString())
cloud.setSubsetString("")

Breakdown: QGIS expression syntax here — IN with a parenthesised list, single =, AND/OR — because this is the layer's own filter and not PDAL's. It applies as the index is walked, so a heavily filtered view actually renders faster than the unfiltered one. It is not persisted in a project file for point cloud layers in every release, so treat it as a session-level tool rather than as configuration.

For committing, write the file:

ground = processing.run("pdal:filter", {
    "INPUT": clean,
    "FILTER_EXPRESSION": "Classification == 2",
    "OUTPUT": "/data/clean/tile_0345_ground.copc.laz",
})["OUTPUT"]

Breakdown: pdal:filter also accepts an EXTENT parameter, which crops in the same pass — worth using when you know you want both, because it halves the number of full reads. The output is a genuine point cloud file with a new header whose point count reflects the filter, which matters if anything downstream reads the count rather than the points.

Filtering on return number and intensity

Classification is not the only useful attribute, and on an unclassified delivery it is not even an available one. Two others carry real information.

ReturnNumber and NumberOfReturns describe how a single laser pulse fragmented on its way down. A pulse that hits bare tarmac returns once: return 1 of 1. A pulse that clips a branch, then a lower branch, then the ground returns three times, and the last of those is far more likely to be ground than the first. That gives a usable rough ground proxy with no classification at all.

last_returns = processing.run("pdal:filter", {
    "INPUT": clean,
    "FILTER_EXPRESSION": "ReturnNumber == NumberOfReturns",
    "OUTPUT": "/data/clean/tile_0345_last.copc.laz",
})["OUTPUT"]

Breakdown: Comparing two attributes rather than an attribute and a constant is legal in PDAL expressions and is the whole trick here. The result is not a ground classification — a last return off a flat roof is still a roof — but over vegetated terrain it removes most of the canopy, and it is often enough to produce a usable first-pass surface while you chase the supplier for the classified delivery. Combining it with a classification test where classes do exist (ReturnNumber == NumberOfReturns && Classification != 6) tightens it considerably.

Intensity records how much energy came back, and it separates materials that sit at the same height: water and wet asphalt return very little, dry sand and painted road markings return a lot. It is not calibrated between surveys, or reliably between flight lines within one survey, so a threshold that works on one tile may not work on the next. Read the range from the statistics before choosing one.

low = stats.minimum("Intensity")
high = stats.maximum("Intensity")
print(f"intensity spans {low} to {high}")

bright = processing.run("pdal:filter", {
    "INPUT": clean,
    "FILTER_EXPRESSION": f"Intensity > {low + 0.8 * (high - low)}",
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

Breakdown: Deriving the threshold from the file's own range rather than hard-coding a number is what makes this survive being run on a second tile. Note the f-string: the filter expression is just a string, so building it from computed values is normal and there is no parameter binding to worry about. Because intensity is uncalibrated, treat anything you derive from it as a candidate for review rather than as an answer — it is a good way to find the water bodies in a tile and a poor way to map them.

When the survey was never classified

If classesOf("Classification") returns [1], no class filter will help, and no amount of PyQGIS will retroactively classify the cloud — QGIS ships no ground-classification algorithm. There are three honest routes.

The first is to go back to the supplier: classification is normally part of the deliverable and an unclassified delivery is often simply a mistake in what was sent. The second is to run PDAL's own ground filter outside QGIS, through the pdal executable with a pipeline that applies an SMRF or PMF stage, and bring the result back as a normal file. The third, viable when you only need a surface model and not a terrain model, is to skip classification entirely and grid the maximum return per cell — no class filter required.

What to do with an unclassified deliveryWhen every point is class one, three routes remain. Ask the supplier for the classified delivery. Run a ground classification stage in PDAL outside QGIS and bring the result back. Or accept that a surface model, built from the highest return per cell, needs no classification at all.classesOf("Classification") returned [1]nothing is classifiedask the supplierclassification is usuallypart of the contractcheapest by farclassify outside QGISa PDAL pipeline with anSMRF or PMF ground stageneeds tuning per landscapebuild a surface modelmaximum return per cell,no classes neededworks today, answers lesswhat does not work: pretending class 1 is class 2

QGIS version compatibility

Subset strings on point cloud layers require 3.26 or newer; pdal:filter requires 3.32. Both are present in 3.34 LTR and unchanged through 3.44. The EXTENT parameter on pdal:filter was added alongside the algorithm and takes a rectangle string in the layer's CRS, not a layer — passing a layer there fails validation rather than silently ignoring it.

Troubleshooting

  • The filtered output has zero points. The expression matched nothing. Print classesOf("Classification") and check for = where PDAL wants ==.
  • The subset string had no effect on an exported raster. Expected — pdal: algorithms read the file and ignore the layer's view. Repeat the condition in FILTER_EXPRESSION.
  • Classification IN (2,9) is rejected by the algorithm. IN is QGIS syntax. PDAL wants Classification == 2 || Classification == 9.
  • The output file is bigger than the input. The output was written uncompressed, or to a different point format. Keep the .copc.laz extension.
  • Filtering is slower than expected on a virtual point cloud. Every member tile is read. Add an EXTENT so tiles outside it are skipped entirely.
  • Class 2 exists but the terrain model still has roofs in it. The survey's ground classification is poor, which happens on dense urban blocks. Compare against the surface model before trusting it.

Conclusion

Read the classes before writing the filter, drop noise in its own pass, explore with a subset string and commit with pdal:filter, and remember which expression language you are in. When a delivery is unclassified, say so out loud rather than working around it — the workarounds all quietly change what the output means.

Frequently Asked Questions

Can I edit the classification values in place? No. Point cloud layers are read-only in QGIS, and there is no wrapped algorithm that assigns classifications. Producing a reclassified cloud means writing a new file from a PDAL pipeline run outside QGIS.

Does filtering speed up rendering? A subset string does, because fewer points are read from the index. A written filtered file does too, and additionally makes every subsequent algorithm faster, which is why it is worth doing once for a cloud you will use repeatedly.

How do I filter by something other than classification? Any attribute in the file works in both languages — Intensity, ReturnNumber, NumberOfReturns, ScanAngleRank, GpsTime. Filtering ReturnNumber == NumberOfReturns keeps last returns, which is a useful rough proxy for ground when nothing is classified.

Is there a way to preview the effect before writing a file? Yes, and it is the whole reason subset strings are worth knowing: set the equivalent QGIS expression, look at the canvas, then translate it to PDAL syntax once you are happy.