It is easy to generate Digital Surface Models (DSMs) from the point clouds generated by dense-matching photogrammetry. The blast2dem tool that is part of the BLAST extension of LAStools can seamlessly triangulate and raster up to 2 billion points. In an earlier blog article we showed how to also create Digital Terrain Models (DTMs) from photogrammetric points by combining lassort, lasground, and las2dem. We were since asked how to do this for larger inputs. Here the answer.
In the following we describe a LAStools pipeline that turns the point cloud shown above into a DTM, a DSM, and 1 meter elevation contours using tile-based multi-core batch processing. Below you see the result.
The ‘macmahanfarm.laz’ file has 140,616,508 RGB-colored points and was provided by Pix4D as a LAZ file of 1.2 GB that decompresses to 4.5 GB of LAS. The unusually low LASzip compression ratio of “only” 1:4 comes from the incoherent order of the points in the file and their excessive coordinate resolution.
First we use lastile to tile the points of ‘macmahanfarm.laz’ into 500 meter by 500 meter tiles with 25 buffer. This buffer will help to avoid edge artifacts along tile boundaries. A lasinfo report tells us that point coordinates are stored with millimeter resolution in the original file. Because the data is not that precise we rescale it on-the-fly to centimeters, which improves compression and eases I/O bottlenecks.
lastile -i macmahanfarm.laz ^ -rescale 0.01 0.01 0.01 ^ -tile_size 500 -buffer 25 ^ -o poop\mmfarm.laz -olaz
One oddity in the data obtained from Pix4D is the high spatial incoherence in the ordering of points in the original file. Here an illustration of the point order in ‘mmfarm_266500_6175500.laz’, one of the 65 tiles generated in the last step.
Such incoherence in point order slows down subsequent processing with lasground and las2dem (which benefit from spatial locality). We use lassort to reorder the tiles on 4 cores in parallel and add appendix ‘_s’ to each file name:
lassort -i poop\mmfarm*.laz ^ -odix _s -olaz ^ -cores 4
This sort the points along a standard Hilbert or Morton space filling curve. Coherence in point order also improves compression: the sorted tile files are about 30 percent.smaller than the unsorted ones.
Next we use lasground to classify the 65 sorted tiles on 4 cores in parallel. A step size of 20 was chosen after neither 10 nor 15 were able to remove a few dense patches of vegetation. Alternatively, use a smaller step size and reclassify points manually with lasview.
lasground -i poop\mmfarm*_s.laz ^ -step 20 -extra_fine ^ -odix g -olaz ^ -cores 4
Below a closer look at the result for tile ‘mmfarm_266500_6175000_sg.laz’. Note that our data is mostly open terrain. Your milage will vary (a lot!!!) depending on the amount of vegetation in the scene. There will generally not be any elevation data in vegetated areas. For larger forested patches you will need LiDAR to get a DTM. The same is true if you are particularly interested in bare-earth terrain features beneath trees and shrubs.
Then we use las2dem to raster the ground points of each tile as a DTM with 0.5 meter resolution in the BIL format, again on 4 cores in parallel. The ‘-use_tile_bb’ limits rasterization to the original 500 meter by 500 meters of the tile, clipping the 25 meter buffer along the tile boundaries from the output.
las2dem -i poop\mmfarm*_sg.laz ^ -keep_class 2 ^ -step 0.5 -use_tile_bb ^ -ocut 3 -odix _dtm -obil ^ -cores 4
lasgrid -i poop\mmfarm*_dtm.bil -merged ^ -step 0.5 ^ -o macmahanfarm_dtm.tif
The DSM is created the same way but without filtering out points. Here a trick how to quickly generate contours that are not too complex (but also not of catographic quality). We first use lasgrid to merge and average the DTMs into a smoother and coarser DTM from which we extract 1 meter elevation contours with blast2iso.
lasgrid -i poop\mmfarm*_dtm.bil -merged ^ -step 1.0 -average ^ -o poop\dtm_coarse.bil
blast2iso -i poop\dtm_coarse.bil ^ -iso_every 1.0 ^ -smooth 2 -simplify_length 0.5 -simplify_area 0.5 -clean 15 ^ -o macmahanfarm_iso.shp
You can find an example batch script that summarizes all the above steps into one workflow in the latest LAStools distribution.
Finally, a word on compression: lowering the coordinate resolution to centimeters, sorting the points of every tile into a coherent spatial order, and merging the tiles back into one, results in a ‘macmahanfarm_cm.laz’ file that is only 0.56 GB: less than half the size of the ‘macmahanfarm.laz’ file we obtained from Pix4D. Now the LASzip compression ratio is 1:8.