How to generate a clean Digital Terrain Model (DTM) from point clouds that were generated with the image matching techniques implemented in various photogrammetry software packages like those from Pix4D, AgiSoft, nframes, DroneDeploy and others has become an ever more frequent inquiry. There are many other blog posts on the topic that you can peruse as well [1,2,3,4,5,6,7,8]. In the following we go step by step through the process of removing low noise from a high-density point cloud that was generated with Pix4D software. A composite of the resulting DTM and DSM is shown below.
Final DSM and DTM created with LAStools for a photogrammetric point cloud of a road generated by Pix4D.
After downloading the data it is useful to familiarize oneself with the number of points, the density of points and their geo-location. This can be done with lasview, lasinfo, and lasgrid using the command lines shown below. There are around 19 million points in the file and their density averages around 2300 points per square meter. Because the RGB values have a 16 bit range (as evident in the lasinfo report) we need to add the option ‘-scale_rgb_down’ to the command line when producing the RGB raster with lasgrid.
lasview -i 0_photogrammetry\densified_point_cloud.laz
lasinfo -i 0_photogrammetry\densified_point_cloud.laz ^
-cd ^
-o 1_quality\densified_point_cloud.txt
lasgrid -i 0_photogrammetry\densified_point_cloud.laz ^
-scale_rgb_down ^
-step 0.10 ^
-rgb ^
-fill 1 ^
-o 1_quality\densified_point_cloud.png
The first step is to use lastile and create smaller and buffered tiles for these 19 million photogrammetry points. We use a tile size of 100 meters, request a buffer of 10 meters around every tile, and flag buffer points as withheld so they can be easily be dropped later. We also make sure that all classification codes are reset to 0.
lastile -i 0_photogrammetry\points.laz ^
-set_classification 0 ^
-tile_size 100 -buffer 10 -flag_as_withheld ^
-o 2_tiles_raw\seoul.laz -olaz
We start with lassort as a pre-processing step that rearranges the points into a more coherent spatial order which often accelerates subsequent processing steps.
lassort -i 2_tiles_raw\seoul_*.laz ^
-odir 3_tiles_temp0 -olaz ^
-cores 4
Next we use a sequence of four modules, namely lasthin, lasnoise, lasground, and lasheight with fine-tuned parameters to remove the low noise points that are typical for point clouds generated from imagery by photogrammetry software. A typical example for such noise points are shown in the image below generated with this call to lasview:
lasview -i 3_tiles_temp0\seoul_210400_542900.laz ^
-inside 210406 542894 210421 542921 ^
-points 20000000 ^
-kamera 0 -95 90 0 -0.3 1.6 ^
-point_size 4
Ground surface noise (exaggerated by pressing <]> in lasview which doubles the scale in z).
As always, the idea is to construct a surface that is close to the ground but always above the noise so that it can be used to declare all points beneath it as noise. Below is a processing pipeline whose parameters work well for this data and that you can fine tune for the point density and the noise profile of your own data.
First we use lasthin to give those points the classification code 8 that are closest to the 70th percentile in elevation within every 20 cm by 20 cm cell. As statistics like percentiles are only stable for a sufficient number of points we only do this for cells that contain 25 points or more. Given that we have an average of 2300 points per square meter this should easily be the case for all relevant cells.
lasthin -i 3_tiles_temp0\seoul_*.laz ^
-step 0.20 ^
-percentile 70 25 ^
-classify_as 8 ^
-odir 3_tiles_temp1 -olaz ^
-cores 4
The we run lasnoise only points on the points with classification code 8 and reclassify all “overly isolated” points with code 9. The check for isolation uses cells of size 20 cm by 20 cm by 5 cm and reclassifies the points in the center cell when the surrounding neighborhood of 27 cells has only 3 or fewer points in total. Changing the parameters for ‘-step_xy 0.20 -step_z 0.05 -isolated 3’ will remove isolated points more or less aggressive.
lasnoise -i 3_tiles_temp1\seoul_*.laz ^
-ignore_class 0 ^
-step_xy 0.20 -step_z 0.05 -isolated 3 ^
-classify_as 9 ^
-odir 3_tiles_temp2 -olaz ^
-cores 4
Next we use lasground to ground-classify only the surviving points (that still have classification code 8) by ignoring those with classification codes 0 or 9. This sets their classification code to either ground (2) or non-ground (1). The temporary surface defined by the resulting ground points will be used to classify low points as noise in the next step.
lasground -i 3_tiles_temp2\seoul_*.laz ^
-ignore_class 0 9 ^
-town -ultra_fine -bulge 0.1 ^
-odir 3_tiles_temp3 -olaz ^
-cores 4
Then we use lasheight to classify all points that are 2.5 cm or more below the triangulated surface of temporary ground points as points as noise (7) and all others as unclassified (1).
lasheight -i 3_tiles_temp3\seoul_*.laz ^
-classify_below -0.025 7 ^
-classify_above -0.225 1 ^
-odir 4_tiles_denoised -olaz ^
-cores 4
The progress of each step is illustrated visually in the two image sequences shown below.
all points.
after lasthin: all 70th elevation percentile points per 20 cm by 20 cm cell.
same. all 70th elevation percentile points per 20 cm by 20 cm cell.
after lasnoise with isolated points in blue.
after lasnoise without isolated points.
after lasground.
after lasheight: using surface defined by ground points to classify low noise.
low noise only.
points without low noise.
all points.
temporary ground surface for classifying low noise.
low noise points.
points without low noise.
Now that all noise points are classified we start a standard processing pipeline, but always ignore the low noise points that are now classified with classification code 7.
The processing steps below create a 10 cm DTM raster. We first use lasthin to classify the lowest non-noise point per 10 cm by 10 cm cell. Considering only those lowest points we use lasground with options ‘-town’, ‘-extra_fine’, ‘-bulge 0.05’, and ‘-spike 0.05’. Using las2dem the resulting ground points are interpolated into a TIN and rasterized into a 10 cm DTM cutting out only the center 100 meter by 100 meter tile. We store the DTM raster as a gridded LAZ for maximal compression and finally merge these gridded LAZ files to create a hillshaded raster in PNG format with blast2dem.
lasthin -i 4_tiles_denoised\seoul_*.laz ^
-ignore_class 7 ^
-step 0.10 ^
-lowest ^
-classify_as 8 ^
-odir 5_tiles_thinned_lowest -olaz ^
-cores 4
lasground -i 5_tiles_thinned_lowest\seoul_*.laz ^
-ignore_class 1 7 ^
-town -extra_fine ^
-bulge 0.05 -spike 0.05 ^
-odir 6_tiles_ground -olaz ^
-cores 4
las2dem -i 6_tiles_ground\seoul_*.laz ^
-keep_class 2 ^
-step 0.10 ^
-use_tile_bb ^
-odir 7_tiles_dtm -olaz ^
-cores 4
blast2dem -i 7_tiles_dtm\seoul_*.laz -merged ^
-hillshade ^
-step 0.10 ^
-o dtm_hillshaded.png
The processing steps below create a 10 cm DSM raster. We first use lasthin to classify the highest non-noise point per 10 cm by 10 cm cell. With las2dem the highest points are interpolated into a TIN and rasterized into a 10 cm DSM cutting out only the center 100 meter by 100 meter tile. Again we store the raster as gridded LAZ for maximal compression and merge these files to create a hillshaded raster in PNG format with blast2dem.
lasthin -i 4_tiles_denoised\seoul_*.laz ^
-ignore_class 7 ^
-step 0.10 ^
-highest ^
-classify_as 8 ^
-odir 8_tiles_thinned_highest -olaz ^
-cores 4
las2dem -i 8_tiles_thinned_highest\seoul_*.laz ^
-keep_class 8 ^
-step 0.10 ^
-use_tile_bb ^
-odir 9_tiles_dsm -olaz ^
-cores 4
blast2dem -i 9_tiles_dsm\seoul_*.laz -merged ^
-hillshade ^
-step 0.10 ^
-o dsm_hillshaded.png
The final result is below. The entire script is linked here. Simply download it, modify it as needed, and try it on this data or on your own data.
final DSM.
final DTM.
