In this article, you’ll find step-by-step instructions on how to create a DTM from an unclassified and classified LIDAR cloud, using LAStools and QGis. We’ll look at the LASground and LASground_new tools, las2dem and lastile.
Tutorial HD LIDAR data processing with QGis
2- Download LIDAR HD data from IGN and load it into QGis
3- Tools for LIDAR data in QGis 3.32
4-Colorize a point cloud from an orthophoto
5-Colorize from an image with LAStools
6- Digital surface model (DSM) with CloudCompare and LAStools
7- digital terrain model (DTM) with CloudCompare
8- Creating a Digital Terrain Model with LAStools
The process for creating a DTM from an unclassified LIDAR point cloud is the same as that described in the previous article:
1- The points in the cloud must be classified into “ground” points and “non-ground” points.
2- Interpolate a raster from “ground” points
If you’ve read the previous article and found DTM generation fairly straightforward, think again! Here, simplicity is not the order of the day.
- You have 2 tools for classifying points that are similar but different
- The number of possible parameters for each tool exceeds a dozen.
- The best results are obtained by running the tools twice, with different parameters.
- Whether you have a LAStools license or not, an IGN LIDAR cloud cannot be processed without splitting it into smaller clouds.
It’s impossible to see all the possibilities here. It’s up to you to take your courage in both hands and explore all the ones we won’t show you.
Classifying ground and off-ground points with LAStools
LASGround and LASGround_New are two tools in the LAStools toolbox that are used to extract ground points from a LiDAR point cloud. However, there are significant differences between these two tools in terms of methods and functionality.
- LASGround:
LASGround is a tool that identifies and classifies ground points in a LiDAR point cloud. It mainly uses algorithms based on morphological operations to detect ground points. LASGround uses a grid of cells and a variable-size window to identify ground points based on the density of neighboring points. It is effective for extracting ground points in flat or gently sloping areas.
This tool works very well in natural environments such as mountains, forests, fields, hills or any other terrain with few man-made objects. - LASGround_New:
LASGround_New is an enhanced version of LASGround, introduced to provide better performance and more accurate extraction of ground points. Unlike LASGround, LASGround_New uses a machine learning approach. It can use methods such as Random Forest to classify ground points based on features such as height, density, etc. LASGround_New is designed to handle more complex environments, including areas with steeper slopes and varied topographical features.
It’s a completely redesigned version of lasground that does a much better job of handling complicated terrain where there are steep mountains close to urban areas with many buildings.
LASground
This is a tool for extracting bare ground: it classifies LIDAR points into ground points (class = 2) and above-ground points (class = 1). This tool works very well in natural environments such as mountains, forests, fields, hills or any other terrain with few man-made objects.
The tool also produces excellent results for cities, but buildings larger than the step size can cause problems.
Step size
- The default step size is 5 meters, which is good for forests or mountains.
- For cities or flat ‘-towns’, the step size is increased to 10 meters.
- For cities or warehouses, ‘-city’ increases the step size to 25 meters. For very large cities, use ‘-metro’ and the step size will be increased to 50 meters.
- You can also set it directly with ‘-step 35’.
Experienced users can refine the algorithm by specifying the threshold in meters above which peaks are suppressed. The ‘-spike 0.5+5’ parameter will suppress ascending peaks greater than 50 centimetres and descending peaks less than 5 metres in the coarsest TIN.
Another interesting parameter is ‘-bulge 1.0’, which specifies how much the TIN is allowed to bulge when including points as it becomes more refined. The default bulge is one tenth of a pitch for pitch sizes greater than 5 meters, and one fifth of a pitch in all other cases.
The maximum standard deviation for planar patches in centimetres can be set with ‘-stddev 10’. The maximum offset in meters up to which points above the current ground estimate are included can be set with ‘-offset 0.1’.
LASground_new
This is a tool for extracting bare ground: it classifies LIDAR points into ground points (class = 2) and non-ground points (class = 1). This is a totally reworked version of lasground, which handles much better complicated terrain where there are steep mountains close to urban areas with many buildings.
The default step size has been increased to 25 meters, allowing most buildings in common cities to be removed (with the exception of very large warehouses or factories, which require a larger step).
For small towns or terrains with few buildings, use the ‘-town’ option, which reduces the step size to 10 meters. For very large towns, use the ‘-metro’ option and the step size will be increased to 50 meters.
You can also experiment by directly setting the step size using ‘-step 35’.
Then there’s the ‘-nature’ option, which uses a step size of 5 meters for terrain without buildings, and the ‘-wilderness’ option, which uses 3 meters if you’re interested in smaller terrain features in high-resolution LiDAR.
Experienced users can refine the algorithm by specifying the threshold in meters above which peaks are suppressed. The ‘-spike 0.5’ parameter suppresses ascending peaks greater than 50 centimetres and the ‘-spike_down 2.5’ parameter suppresses descending peaks less than 2.5 metres in the coarsest TIN.
The maximum offset in meters up to which points above the current ground estimate are included can be set with ‘-offset 0.1’. Alternatively – if you want to try out a large number of parameters – you can play with ‘-bulge 1.5’ which – by default – is set to a tenth of the step size and then tightened in the range 1.0 to 2.0.
This applies to both tools
You can use classifications other than ASRPS with ‘-ground_class 10’ or ‘-non_ground_class 25’. All points not classified as ground can be left unchanged with ‘-non_ground_unchanged’.
If you wish to exclude certain classifications from the bare earth calculation, you can do so with ‘-ignore_class 7’. For very steep hills, you can intensify the search for initial ground points with ‘-fine’, ‘-extra_fine’, ‘-ultra_fine’ or ‘-hyper_fine’. Similarly, for very flat terrain, you can simplify the search with ‘-coarse’ or ‘-extra_coarse’, but try the default setting first.
Finally, you can ask lasground or lasground_new to calculate the height above ground for each point (so that you can use lasclassify afterwards without needing to run lasheight first) with ‘-compute_height’ or even ask for the calculated height to replace the elevation value with the ‘-replace_z’ option. The result is a height-normalized LAS/LAZ file that can be used, for example, with the lascanopy or lasgrid tools or the CHM (canopy height model) algorithm without a pit. If lasground or lasground_new malfunctions, try disabling certain optimizations using the ‘-no_clean’ or ‘-no_bulge’ flags.
If there are too few points to perform reliable terrain classification, the files are simply copied (and in the case of ‘-replace_z’, all elevations are set to zero). These files can also be ignored using the ‘-skip_files’ command line option.
Pre-processing the example cloud with LAStile
We’re going to use the same LIDAR point cloud as in the previous article, a classified cloud to which we’ve removed the classification. This allows us to compare the results of the tools with the DTM made directly with the IGN “ground” classification.
The first problem we have is that this point cloud contains 28 million points. If we use the LAStools tools in this article on this cloud, we get a message:
ERROR processing file 'C:\telechargements\ncLHD_FXX_0722_6403_PTS_C_LAMB93_IGN69.laz'. maybe out of memory?So we need to split our cloud into smaller clouds (less than 1.5 million if you don’t have a license) and use “directory” tools instead of “file” tools. We’ll reconstruct a single output at the end of processing.
To create the tiles, we create an empty directory, then open the LAStile tool
You’ll need to calculate the “tile size” so that the tiles have the desired number of points. Our input cloud is 1000mx1000m and contains 28 million points. If we want less than 1.5 million points per tile, 25 tiles of 200mx200m will have around 1.1 million points. We’ll then enter 200 in tile size.
On the other hand, the Output LAS/LAZ field, which is marked optional, is not optional at all, otherwise you’ll have to search for it in the recesses of /Users/username/AppData/Local/Temp/. You need to enter a directory for the tiles and a name, e.g. “tile.las”. The tiles will then be named “tuile_coordonnéeX_coordonnéeY.las”.
Our treatment of these tiles includes interpolations between points. If we had a classic edge-to-edge tiling, we’d have edge effects: points located on the edges wouldn’t have neighbors located on the adjacent tile. To avoid this, LAStools adds a buffer zone around each tile. When the resulting cloud is reconstituted, these will be removed automatically, but there will be no edge effects in intermediate processing.
Extraction of ground points
We load the point cloud to be processed, in its colored version, as this will be more practical for analyzing the results.
We have an urban area surrounded by natural zones. We’re going to use the LASground tool, configuring it with the parameters -city for terrain type and -ultra_fine for preprocessing. This is the classic LIDAR processing configuration for archaeology and it corresponds very well to our LIDAR cloud. Don’t forget to check the 64bits executable box.
The result is a series of clouds with two classes:
- 1- unclassified (above-ground points)
- 2- ground
This is where we see the biggest difference from processing with CloudCompare. To evaluate the result we have two problems:
- the first is the fact that we have divided our cloud into multiple batches of data, making it difficult to get an overall view.
- the second is the lack of a true 3D interface that lets you see the result from all angles. The 3D view in QGis is too coarse to see the result properly.
The only way left is to use the Elevation Profile view, but this requires us to draw several transects to detect any problems.
Depending on the types of terrain present in our initial cloud, it will be necessary to check several tiles.
In our example, we’ll stop here. But it’s quite common to need a second lasground pass. The classic case is that of archaeological research. Let’s suppose we have an interesting area in the natural part of the image in our example point cloud. We have used the -city option to classify the buildings above ground. This option uses a pitch of 25 meters. If we have small structures, there’s a good chance they won’t be classified as above ground. A second run of lasground, this time with the -wilderness option (natural area without buildings) and applying it exclusively to points classified as “ground” by the first run, could then remove the small structures.
DTM calculation with LAS2dem
This tool reads LIDAR points, temporarily triangulates them into a TIN, then layers the TIN onto a DTM raster. The tool can use ‘-elevation’, ‘-slope’, ‘-intensity’, ‘-rgb’ values, or ‘-hillshade’ or ‘-gray’ or ‘-false’ coloring. Output can be in BIL, ASC, IMG, FLT, XYZ, DTM, TIF, PNG or JPG format . The particular range of values to be displayed can be limited using ‘-set_min_max 10 100’ or defined using ‘-compute_min_max’. When using filters such as ‘-last_only’ or ‘-keep_class 2’, you can use the ‘-extra_pass’ option to first determine the number of points to be triangulated. This saves memory.
By default, the size of the generated raster is based on the extent of the bounding box. If the LAS/LAZ file has been generated using lastile, as in our example, its extent can be set to that of the tile using the ‘-use_tile_bb’ option. Any “border buffer” that the tile may have had is then not rasterized. This avoids border artifacts while creating corresponding tiles in parallel.
A KML file is automatically generated for displaying the DEM in Google Earth (for TIF/PNG/JPG). As the LAS/LAZ file contains projection information, this is used to georeference the KML file.
For our example, here’s the las2dem dialog window:
The keep_class 2 filter means that only points classified as “ground” (code 2) will be processed.
We set the DTM pitch to 0.5 meters.
We check the “use tile bounding box” so that the raster output of each tile is cut to keep only the original tile, without the buffer zone created by lastile.
-extra_pass optimizes memory allocation.
The result is the creation of three files per tile: the tif file, the associated world file and the kml file.
Note that the output tiles are joined (without overlap):
LAStools has the tool to tile the starting point cloud, but no tool to reconstitute the final DTM. We need to use the standard QGis tool: menu-> raster -> various -> merge
We select all the tiles, name the output tif file and execute the command.
We can now see the complete result and compare it with the reference MNT made from the IGN classification.
In the natural part of our LIDAR zone, the two results match:
In the urbanized area, however, there are differences:
Logically, we’d have to go back to the first lasground pass, change the pitch to a value of 10 metres (-town) and do it all over again…
