Point Cloud & LiDAR Workflows in PyQGIS

A LiDAR survey arrives as a folder of .laz files and a specification that says "ground and buildings, 1 m DEM, tiled". Nothing in that sentence is hard. What makes point clouds different from every other dataset in QGIS is scale: a single tile holds more points than the whole rest of your project holds features, and the API is built around that fact rather than around iterating features one at a time.

This guide sits inside Spatial Data Processing & Automation and covers the whole path — getting a cloud into a layer, reading what is actually in it, filtering and clipping it, and converting it into the raster or vector products that the rest of a workflow can consume. The recurring theme is that you almost never touch individual points from Python. You describe what you want and let an indexed reader or a PDAL algorithm do the walking.

How a LiDAR file becomes something you can useRaw LAS or LAZ files are indexed into a COPC or EPT structure. QGIS opens the indexed cloud as a point cloud layer, which then feeds three routes: a renderer for display, the PDAL processing algorithms for filtering and conversion, and exported raster or vector products for downstream analysis.Indexing happens once; everything downstream depends on itraw LAS / LAZone file per tilespatial indexCOPC or EPT octreepoint cloud layerextent · attributes · statsrenderer2D canvas and 3D viewfilter and clipsubset string · pdal:clipexport rasterDEM · DSM · densityexport vectorboundary · thinned pointsan unindexed LAZ opens, then re-indexes itself on every single run

What a point cloud layer actually is

QgsPointCloudLayer is not a vector layer with a lot of features. It has no getFeatures(), no attribute table you can edit, and no editing buffer. What it has is an index: a hierarchical octree that stores points in nodes, coarse near the root and dense at the leaves, so that a renderer or an algorithm can ask for "enough points to fill this screen at this zoom" and read a few megabytes instead of a few gigabytes.

That is why the format matters more here than anywhere else in QGIS. A plain .las or .laz file has no index. QGIS will still open it, but it has to build an index first, and it does that in a background task that writes into a cache directory. Do that inside a batch script and you pay the indexing cost on every run, on every machine, for every file.

The practical consequence is that the first decision in any LiDAR project is a storage decision, not a code decision. Convert the delivery once, keep the converted copies, and point every script at those. A survey delivered as three hundred .laz tiles becomes three hundred .copc.laz tiles occupying roughly the same disk space, after which every subsequent read is a seek into an octree rather than a linear decompression of the whole file. The conversion is one algorithm call per file and it parallelises perfectly, so it is worth doing even when the immediate task only needs two tiles.

It also changes what "open the layer" costs. Adding an indexed cloud to a project reads a header and a handful of octree nodes — a few hundred kilobytes regardless of whether the file is 40 MB or 4 GB. Adding an unindexed one blocks on a background task whose duration is proportional to the file size, and a script that calls isValid() immediately afterwards may see a layer that is valid but has not finished building its index, so the first render is empty and the second is not.

from qgis.core import QgsPointCloudLayer

cloud = QgsPointCloudLayer("/data/lidar/tile_0345.copc.laz", "tile_0345", "copc")
print(cloud.isValid(), cloud.pointCount())
print(cloud.crs().authid(), cloud.extent().toString(0))

Breakdown: The third argument is the provider key and it is not optional the way it is for a file-based vector layer — copc for a Cloud Optimised Point Cloud, ept for an Entwine index (you pass the ept.json), and pdal for a raw .las/.laz that has to be indexed on the fly. Passing the wrong key gives you an invalid layer with an unhelpful message, so it is worth branching on the file extension rather than guessing. pointCount() reads the header and is instant even on a billion-point file.

Reading the header before you read the points

Almost every point cloud bug is a mismatch between what you assumed the file contains and what it actually contains: the classification scheme is non-standard, the intensity is 8-bit rather than 16-bit, the vertical datum is not what the metadata claims, or the file is simply empty over your area of interest. All of that is answerable from the header and the index statistics without touching a point.

attributes = cloud.attributes()
for attribute in attributes.attributes():
    print(attribute.name(), attribute.size(), "bytes")

stats = cloud.statistics()
print("Z range", stats.minimum("Z"), stats.maximum("Z"))
print("classes present", sorted(stats.classesOf("Classification")))

Breakdown: attributes() returns a QgsPointCloudAttributeCollection; its own attributes() method gives the list, and every LAS file will have at least X, Y, Z, Classification, Intensity and ReturnNumber. statistics() is computed during indexing and cached, so it is cheap for COPC and EPT and expensive the first time for an unindexed file. classesOf() is the single most useful call in this whole guide: it tells you whether the survey was actually classified, or whether every point is class 1 (unassigned) and your ground filter is about to return nothing.

Two header fields deserve a specific check in any script that will run unattended. The first is the CRS. LAS files carry their projection in a variable-length record that older producers wrote inconsistently, and QGIS reports an empty authid() when it cannot make sense of it — at which point the layer lands wherever the project CRS happens to put it, which is usually the Gulf of Guinea. The second is the vertical datum, which the LAS header does not really model at all: a file whose Z values are ellipsoidal heights and a file whose Z values are orthometric heights look identical to the API and differ by tens of metres on the ground.

if not cloud.crs().isValid():
    raise SystemExit("no CRS on the LAS header — run pdal:assignprojection first")

if cloud.pointCount() == 0:
    raise SystemExit("empty tile — the survey boundary probably does not cover this cell")

Breakdown: Failing loudly at the top of a script is worth far more than it looks, because the alternative failure is silent and downstream. An empty tile propagates as a raster full of nodata, which propagates as a hole in a mosaic, which somebody notices three steps later in a hillshade. Both of these checks cost a header read.

Classification is the attribute everything hangs off

The ASPRS classification codes are a small integer vocabulary baked into the LAS specification, and every ground-filter, building-extraction and canopy-height workflow is a filter on that one attribute. Knowing the codes by heart saves an enormous amount of time.

The classification codes you will actually useA schematic cross-section of a landscape shows returns falling into classes: class two on the bare earth surface, classes three to five rising through vegetation, class six on building roofs, class nine on water, and class seven for noise points above the scene. A filter on Classification is what separates a digital terrain model from a surface model.One integer decides whether you get terrain or rooftopsclass 2 — ground, the surface a terrain model is built from4–5 vegetation6 building9 water7 noise — drop these firstwhat you filter forterrain model → 2surface model → 2, 6, 5canopy height → 5 minus 2building footprints → 6always exclude → 7, 181 = unassigned, never usefulcheck the classes are present before you filter — many surveys ship class 1 only

Filtering on that attribute from Python does not need an algorithm at all. A point cloud layer accepts a subset string, which is applied by the reader as the index is walked, so it costs nothing and it affects rendering and every algorithm that takes the layer as input.

cloud.setSubsetString("Classification = 2")
print(cloud.subsetString())
cloud.setSubsetString("")  # back to every point

Breakdown: The syntax is the QGIS expression subset supported by the point cloud providers — comparisons and boolean operators on attribute names, so Classification IN (2, 9) AND Z < 340 is valid and Classification = 'ground' is not. Setting it to an empty string clears it. Unlike a vector subset string this does not persist into an exported file: it is a view, and if you want a filtered file you need pdal:filter, covered in filtering and classifying a point cloud.

The PDAL algorithms are the stable API

QGIS 3.32 added a family of Processing algorithms with the pdal: prefix, and by 3.34 LTR they cover essentially everything a production LiDAR pipeline needs. They matter for a reason beyond convenience: the low-level index classes (QgsPointCloudIndex, the node identifiers, the request/block objects) have changed shape across several releases, while the algorithm identifiers and parameter names have not. Code written against pdal:exportraster in 3.32 still runs on 3.44. Code written against IndexedPointCloudNode in 3.32 does not.

import processing

dem = processing.run("pdal:exportraster", {
    "INPUT": "/data/lidar/tile_0345.copc.laz",
    "RESOLUTION": 1.0,
    "TILE_SIZE": 1000,
    "FILTER_EXPRESSION": "Classification == 2",
    "OUTPUT": "/data/output/dtm_0345.tif",
})["OUTPUT"]

Breakdown: FILTER_EXPRESSION here is PDAL syntax, not QGIS expression syntax — note the doubled ==. That inconsistency between the layer subset string and the algorithm filter is the single most common trip-up in this area, and it fails quietly by matching nothing rather than raising. TILE_SIZE controls how the algorithm chunks its work; leaving it at the default is fine until you hit memory pressure on a very large input.

There is a second reason to prefer the algorithms, which only becomes obvious on a long job: they report progress and honour cancellation through the standard QgsProcessingFeedback object. A PDAL pipeline invoked directly through a subprocess gives you neither, so a plugin that shells out has no way to draw a progress bar or to stop when the user asks. Passing a feedback object into processing.run() gives both for free, and it is the same object the rest of the site uses for handling processing feedback and errors.

The algorithms also normalise a genuine awkwardness in PDAL itself: most of them accept either a file path or a loaded QgsPointCloudLayer as INPUT, and they resolve a layer back to its underlying file for you. That means the same call works in a script that has just built a layer and in a batch job that only has paths, without a branch.

The full list is worth skimming once: pdal:info, pdal:convertformat, pdal:reproject, pdal:clip, pdal:filter, pdal:merge, pdal:tile, pdal:thinbyradius, pdal:thinbydecimate, pdal:boundary, pdal:density, pdal:exportraster, pdal:exportrastertin and pdal:exportvector. Running the PDAL algorithms from PyQGIS works through the parameters that are easy to get wrong.

Turning points into products

Nothing downstream in a GIS consumes points directly. A hydrologist wants a DEM, a planner wants building footprints, a forester wants a canopy height model. Every one of those is a conversion, and each conversion is a decision about what to do with the cells that have no points in them.

The same cloud, three productsFiltering to ground returns and gridding produces a digital terrain model. Gridding the highest return in each cell produces a digital surface model. Subtracting terrain from surface produces a normalised height model in which trees and roofs stand at their true height above the ground.Terrain, surface and the difference between themterrain model (DTM)Classification == 2, mean per cellbare earth, buildings removedsurface model (DSM)all returns, maximum per cellroofs and canopy includedheight model (DSM − DTM)raster calculator, one expressionheights above ground, terrain flatthe cells with no points in them decide how the output looksunder a dense canopy a 0.5 m ground grid is mostly holes; 1 m or 2 m is usually the honest resolution

Choosing the output resolution is the decision that matters, and the rule of thumb is arithmetic rather than taste: a grid cell needs at least a few points in it to produce a value you can defend. Divide the survey's stated point density into one and take the square root. A 4 points-per-square-metre survey supports a 0.5 m grid over open ground and, once you filter to ground returns under trees, realistically a 1 m or 2 m one. Asking for 0.25 m gives you a raster full of nodata holes and a slope map full of artefacts, as described in creating a DEM from a point cloud.

Once you have the rasters, the rest is ordinary raster work: the raster calculator subtracts one from the other, zonal statistics summarises the result inside parcels, and slope, aspect and hillshade turn the terrain model into something readable.

Rendering without melting the canvas

A point cloud layer has four renderers, and picking the right one is mostly about which attribute carries the meaning. QgsPointCloudClassifiedRenderer colours by classification and is the default when classes are present. QgsPointCloudRgbRenderer uses the Red, Green and Blue attributes if the survey was colourised from imagery. QgsPointCloudAttributeByRampRenderer maps any numeric attribute — Z, Intensity, NumberOfReturns — through a colour ramp. QgsPointCloudExtentRenderer draws only the tile boundary, which is what you want in an overview map with forty tiles loaded.

from qgis.core import QgsPointCloudClassifiedRenderer, QgsPointCloudCategory
from qgis.PyQt.QtGui import QColor

renderer = QgsPointCloudClassifiedRenderer("Classification")
renderer.setCategories([
    QgsPointCloudCategory(2, QColor("#8d6e3f"), "ground"),
    QgsPointCloudCategory(5, QColor("#2f7a3d"), "high vegetation"),
    QgsPointCloudCategory(6, QColor("#8f4b4b"), "building"),
])
renderer.setPointSize(1.4)
cloud.setRenderer(renderer)
cloud.triggerRepaint()

Breakdown: Categories you do not list are not drawn, which is a fast way to get a classification-filtered view without touching the subset string. QgsPointCloudClassifiedRenderer.defaultCategories() returns the full ASPRS set with the standard colours if you would rather start from that and edit. Point size is in millimetres by default; setPointSizeUnit() switches it to pixels or points, and on a dense cloud the difference between 1 mm and 3 mm is the difference between a readable surface and a solid block of colour.

The other half of performance is the point budget — the maximum number of points QGIS will draw in one canvas refresh. It is a global setting rather than a layer property, and a script that renders a cloud to an image needs to raise it or accept a sparse picture. That, along with the interaction with rendering a layer to an image without the GUI, is covered in styling a point cloud renderer.

Scaling past one tile

A survey is never one file. The moment you have forty of them, three things change: you want a single logical layer rather than forty entries in the tree, you want to process only the tiles that intersect your area of interest, and you want the intermediate products to land somewhere predictable.

The virtual point cloud format — a small .vpc JSON file listing member clouds with their extents — solves the first two at once. QGIS opens it as a single layer, and the reader consults the member extents so that a clipped operation touches only the tiles it must. Building one is a single algorithm call, and it is almost always the right first step of a multi-tile pipeline. Clipping and tiling point clouds covers building and consuming them, and the loop that walks a folder is the same loop as running an algorithm over a folder of files.

Where point clouds meet the rest of the project

A point cloud is unusual in a QGIS project because it participates in almost nothing that vector and raster layers participate in. It has no attribute table dialog, no joins, no relations, no editing, no labelling, and it is invisible to most of the native: algorithms. Everything it contributes to a project it contributes through a conversion.

There are three bridges worth knowing. The first is elevation: a point cloud layer has elevation properties (cloud.elevationProperties()) carrying a Z scale and Z offset, which is how you correct a file delivered in feet inside a metric project without rewriting it, and which is what the elevation profile tool reads. The second is the 3D view, where QgsPointCloudLayer3DRenderer gives the layer a place in a 3D map scene alongside extruded vector layers and a terrain surface. The third — the one that carries the most weight in practice — is the export back to vector or raster, after which the data behaves like any other layer and the whole rest of this site applies to it.

That framing is useful because it tells you where to stop. If a task involves per-feature attributes, joins, styling by category or anything that expects a feature id, the answer is not to reach deeper into the point cloud API; it is to export the points you care about and do the work on the resulting layer. pdal:exportvector produces a normal point layer with the LAS attributes as fields, at which point spatial joins and attribute management work exactly as they do everywhere else.

Failure modes worth recognising on sight

Point cloud problems tend to look like nothing at all — an empty canvas, a blank raster, a job that takes four hours — rather than like an exception, so it is worth being able to recognise the handful that account for most of them.

A cloud that renders as a flat grey rectangle is usually the extent renderer, which QGIS falls back to when a layer is loaded but its index is not yet available. A cloud that renders as a thin ribbon along one edge of its extent is a CRS problem: some points carry sensible coordinates and the rest are at the origin, typically because a merge combined tiles in two projections. A DEM export that produces a raster of nodata means the filter expression matched nothing, and nine times in ten the reason is the = versus == difference between QGIS and PDAL expression syntax. A job that runs for hours on a small area is nearly always reading unindexed files, or reading a whole tiled survey because no virtual point cloud told it which tiles overlap the area of interest.

The diagnostic that resolves most of these in one step is pdal:info, which reports the header, the CRS, the point count and the per-attribute statistics for a file without loading it into the project. Running it first, and logging what it says, converts a class of silent failures into a line in a log file.

Key takeaways

  • QgsPointCloudLayer needs an explicit provider key: copc, ept or pdal. Index once, up front, rather than paying for it on every run.
  • Read statistics().classesOf("Classification") before writing any filter. A great many surveys ship entirely as class 1.
  • The layer subset string uses QGIS expression syntax (Classification = 2); the PDAL algorithms use PDAL syntax (Classification == 2). Mixing them fails silently.
  • Prefer the pdal: Processing algorithms over the low-level index classes — the algorithm identifiers are stable across releases and the C++ index API is not.
  • Choose output resolution from point density, not from what looks nice. Ground returns under canopy are far sparser than the survey's headline figure.
  • Build a virtual point cloud over a tiled survey so downstream operations read only the tiles they need.

Frequently Asked Questions

Do I need PDAL installed separately? No. The QGIS installers bundle PDAL, and the pdal: algorithms appear in the Processing registry automatically from 3.32 onwards. If they are missing, the QGIS build is older or was compiled without PDAL support — check with QgsApplication.processingRegistry().algorithmById("pdal:info") returning something other than None.

Can I iterate individual points from Python the way I iterate features? Technically yes, through QgsPointCloudIndex and a QgsPointCloudRequest, but the classes involved have been renamed and reshaped several times across releases and there is no supported iterator that survives a version bump. For anything you have to maintain, express the work as a pdal: algorithm or export to a vector layer first.

Why is my cloud invisible even though the layer is valid? Three usual causes: the point budget is exhausted by another dense layer, the renderer has categories that do not match the classes present, or the CRS is missing from the LAS header so the layer sits in a different part of the world from everything else. pdal:assignprojection fixes the third.

How big is the index compared to the source file? A COPC file is a reorganised LAZ and is generally within about ten per cent of the original size, because the compression is the same. An EPT index is a directory of many small files and is typically a little larger. Neither is a copy on top of the original — COPC replaces it.

Is there a way to edit classifications from PyQGIS? Not in place. Point cloud layers are read-only in QGIS. You reclassify by writing a new file with pdal:filter or an assignment expression, which is what filtering and classifying a point cloud walks through.