This is the final part of a three-part tutorial on how to use LAStools to implement a pipeline that (1) quality checks a newly received set of raw LiDAR flight strips, (2) tiles and prepares the LiDAR for subsequent exploitation, and (3) generates raster and vector derivatives from the LiDAR points such as DTMs, DSMs, contour lines and building footprints with multi-core batch processing and Tin streaming with BLAST.
It is maybe helpful to start with the tutorial on quality checking. But it is mandatory to first complete the tutorial on LiDAR preparation because that is where you generate the cleaned and classified LiDAR tiles in folders ‚tiles_classified‘ and ‚tiles_final‘ that we are using as input in this tutorial.
Generating matching DTM tiles without edge artifacts
Create a new folder called ‚.\lastools\bin\tiles_dtm‘ for storing the bare-earth digital terrain models (DTMs). Then run ‚las2dem.exe‘ in GUI mode and load the 4 tiles from the folder ‚tiles_classified‘. You can either do this by double-clicking ‚las2dem.exe‘ and using the ‚browse …‘ menu or by entering the command below:
C:\lastools\bin>las2dem -i tiles_classified\fitch*.laz -gui
The las2dem tool creates rasters from LiDAR by triangulating the points and then sampling the resulting TIN at the center of each raster pixel. Remember, the tiles in folder ‚tiles_classified‘ have a buffer of 50 meter that was initially created with ‚lastile‘ and that has been carried through ‚lasground‘, ‚lasheight‘, and ‚lasclassify‘. The presence of these buffers allows us now to avoid artifacts along the tile boundaries because our TINs will be 50 meters wider in all directions than the raster that we are sampling. Rasterizing only the extend of the original tile without the surrounding buffer is accomplished with the ‚tile_bb_only‘ command line parameter.
Set the GUI options as shown below to create one DTM per tile. The DTMs will seamlessly fit together. Check the button ‚use tile bounding box‘ for the reason explained above and check the button ‚extra pass‘ which makes two passes over the data. In the first pass it only counts the number of points that pass through the filter so that in the second it can optimize the memory used to triangulate them. Choose ‚bil‘ as the output format. It is much much much more efficient than the ‚asc‘ format many people like to use. Select 4 cores if you have that many and set the output directory to ‚tiles_dtm‘. Use the ‚filter…‘ menu to choose two filters: the first keeps only the ground points with filter ‚keep_classification 2‘ and the second keeps only one point every 0.5 by 0.5 meter area with the filter ‚thin_with_grid 0.5′. The order is really important (excercise: explain why!). Because we create DTMs with the default ’step‘ size of 1 meter there is no sense in constructing a temporary TIN that is more detailed than one ground point per 0.5 by 0.5 meter area. Constructing a large TIN requires a lot of main memory. We should only use as many points as necessary to sample the TIN at our target resolution ’step‘. Adding a ‚thin_with_grid‘ filter with half the ’step‘ size seems always a good choice for ‚las2dem‘.
Press ‚RUN‘ and the ‚RUN‘ window will pop up. You may have to close the ‚output‘ menu before you can actually see the ‚RUN‘ button. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Alternatively, if you prefer to work in the command-line, you can use this equivalent command here.
C:\lastools\bin>las2dem -i tiles_classified\fitch*.laz ^ -keep_class 2 -thin_with_grid 0.5 ^ -extra_pass ^ -use_tile_bb ^ -odir tiles_dtm -obil ^ -cores 4
Generating matching DSM tiles without edge artifacts
Create a new folder called ‚.\lastools\bin\tiles_dsm‘ for storing the first-return digital surface models (DSMs). You could re-use ‚las2dem‘ if you have not closed the application yet. Otherwise re-run it in GUI mode and load the 4 tiles from the folder ‚tiles_classified‘ again. Do this by double-clicking ‚las2dem.exe‘ and using the ‚browse …‘ menu or by entering the command below:
C:\lastools\bin>las2dem -i tiles_classified\fitch*.laz -gui
Set the GUI options as shown below. The only changes are that the output directory name is now ‚tiles_dsm‘ and that the filtering is now keeping first returns instead of ground points. If you re-used the GUI from the last step you need to delete the old filters and enter the new ones in the correct order (exercise: why is the order important?). You can delete old filter entries by double clicking them. Select 4 cores only if you have a computer with that many cores. This will assign tiles to different cores and process them in parallel.
Press ‚RUN‘ and the ‚RUN‘ window will pop up. You may have to close the ‚output‘ menu before you can actually see the ‚RUN‘ button. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Alternatively, if you prefer to work in the command-line, you can use this equivalent command here:
C:\lastools\bin>las2dem -i tiles_classified\fitch*.laz ^ -first_only -thin_with_grid 0.5 ^ -extra_pass ^ -use_tile_bb ^ -odir tiles_dsm -obil ^ -cores 4
Visualizing and processing DEM rasters as if they were points
Here a little trick. All LAStools can read the BIL format „as if“ it was a point format. The other two raster formats that are supported this way are ESRI’s ASC and FUSION’s DTM. But if possible don’t use ASC because it stores the values as ASCII text which is terribly inefficient and much much slower to read than BIL or DTM. Run lasview in GUI mode and load the four BIL rasters from the folder ‚tiles_dsm‘. You need to check the button labelled ‚.bil‘ in the ‚browse …‘ menu to be able to select the BIL rasters as input. Alternatively you can use this command-line:
C:\lastools\bin>lasview -i tiles_dsm\fitch*.bil -gui
Choose the option ’selected file only‘ and pick one of the four tiles. Play with the ‚color_by‘ or ‚render_only‘ options. But these can also be changed once ‚lasview‘ is displaying the point cloud.
Visualize the ‚fitch_273920_4714320.bil‘ from the ‚tiles_dsm‘ folder. After all points have loaded press ‚t‘ on the keyboard to triangulate them. Press ‚-‚ to make the points disappear and press ‚h‘ multiple times to iterate over the possible ways to shade the triangulation.
Using lasgrid to merge DEM tiles into one large DEM raster
Run ‚lasgrid.exe‘ in GUI mode and load the 4 raster DSM tiles in BIL format from the ‚tiles_dsm‘ folder by double-clicking ‚lasgrid.exe‘ and using the ‚browse …‘ menu. Check the button labelled ‚.bil‘ to be able to select the BIL rasters as input. Or use the following in the DOS command-line:
C:\lastools\bin>lasgrid -i tiles_dsm\fitch*.bil -gui
Configure the GUI with the settings shown below. Check the ‚merge files into one‘ button, specify the output file as ‚dsm.tif‘, and specify the northern UTM zone 19 and NAVD88 as the projection information because the BIL files we generated earlier do (currently) not store it.
Press ‚RUN‘. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Or use this in the DOS command-line:
C:\lastools\bin>lasgrid -i tiles_dsm\fitch*.bil -merged ^ -utm 19N -vertical_navd88 ^ -o dsm.tif
Load the merged raster ‚dsm.tif‘ into a GIS package such as QGIS or ArcGIS and repeat the exercise by creating – in the exact same way – a merged ‚dtm.tif‘ from the 4 raster DTM tiles in the ‚tiles_dtm‘ folder.
Using the BLAST extension for seamless generation of rasters
Run ‚blast2dem.exe‘ in GUI mode and load the 4 raster DTM tiles in BIL format by double-clicking ‚blast2dem.exe‘ and using the ‚browse …‘ menu. You need to check the button labelled ‚.bil‘ to be able to select the BIL rasters as input. Or enter this in the DOS command-line window:
C:\lastools\bin>blast2dem -i tiles_dtm\fitch*.bil -gui
Configure the GUI with the settings shown below. Check the ‚merge files into one‘ button, specify the output file as ‚dtm_hillshade_raster.png‘, select ‚hillside shading‘, and choose ‚png‘ as the output format. Use the ‚projection …‘ menu to specify the northern UTM zone 19 and NAVD88 as the projection information because the BIL files we generated earlier do (currently) not store it.
Press ‚RUN‘. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Or use this in the DOS command-line:
C:\lastools\bin>blast2dem -i tiles_dtm\fitch*.bil -merged ^ -hillshade ^ -utm 19N -vertical_navd88 ^ -o dtm_hillshade_raster.png
Whenever ‚blast2dem‘, ‚las2dem‘, ‚lasgrid‘, or ‚lasoverlap‘ generate rasters that are standard color or gray-scale PNG, TIF, or JPG that can be loaded as an image overlay into Google Earth, a KML is automatically generated – assuming there is projection information either available in the input or specified in the command-line. If you double-click the KML file that we just created you will see the hillshading overlaid and nicely aligned with Fitchburg airport.
Repeat the exercise by creating – in the exact same way – a merged and hillshaded ‚dsm_hillshade_raster.png‘ with KML file from the 4 raster DSM tiles in the ‚tiles_dsm‘ folder.
Using the BLAST extension for seamless generation of contours
Run ‚blast2iso.exe‘ in GUI mode and load the 4 BIL raster tiles from the ‚tiles_dtm‘ folder by double-clicking and using the ‚browse …‘ menu or by entering the command below:
C:\lastools\bin>blast2iso -i tiles_dtm\fitch*.bil -gui
Set the GUI as shown below to compute 1 meter contours from on-the-fly merged bare-earth DTM rasters in BIL format. Check the ‚merge files into one‘ button, specify the output file name as ‚dtm_contours_raster_1m.shp‘, set the isovalue extraction to every 1 unit (here: meter), and check the buttons for ’simplify short segments‘, ’simplify small bumps‘ and, ‚clean short lines‘. Read the corresponding README.txt file if you want more details on what these options exactly do. Add the northern UTM zone 19 as the projection information because the BIL files we generated earlier do (currently) not store it. This is less important when generating SHP output but would be required to produce KML instead.
Press ‚RUN‘. If you forgot or need to fix a setting you can directly modify the command-line before you press START. Or use this command-line:
C:\lastools\bin>blast2iso -i tiles_dtm\fitch*.bil -merged ^ -iso_every 1 ^ -simplify_length 0.5 -simplify_area 1 -clean 5 ^ -utm 19N ^ -o dtm_contours_raster_1m.shp
Visualize the result with the convenient ShapeViewer.exe program and check if the 1 meter contours that ‚blast2dem‘ has produced are what you expected:
Repeat the exercise by creating – in the exact same way – a merged and hillshaded ‚dsm_contours_raster_3m.shp‘ with from the 4 raster DSM tiles in the ‚tiles_dsm‘ folder but with a spacing of 3 meters instead of 1 meter between the contour lines.
Extracting building footprints from on-the-fly merged tiles
Run ‚lasboundary.exe‘ in GUI mode and load the 4 classified tiles with:
C:\lastools\bin>lasboundary -i tiles_final\fitch*.laz -gui
Set the GUI as shown below to generate compute polygonal outlines around the buildings points classified in the previous tutorial. Check the ‚merge files into one‘ button to create one seamless output file. Add a ‚keep_classification 6‘ filter that filters out all points except those classified as building. Check the ‚disjoint‘ button – otherwise a funny-looking single polygonal hull will be created connecting all buildings. Set the ‚concavity‘ to 1.5 meters to specify the smallest feature that lasboundary can grow into. This value needs to be 2 to 3 times the average point spacing. If you set it too small, lasboundary will „eat“ its way into the buildings. If you set it too high nearby buildings will be joined together and the outlines will become coarser. Play with a few different settings such as 1.0, 2.0, 5.0, and 10.0 to get an intuitive idea how this parameter affects the result.
Press ‚RUN‘. If you forgot or need to fix a setting you can directly modify the command-line before you press START. Or use the equivalent command-line below directly in the DOS window.
C:\lastools\bin>lasboundary -i tiles_final\fitch*.laz -merged ^ -keep_class 6 ^ -concavity 1.5 -disjoint ^ -o buildings.shp
Visualize the result with the convenient ShapeViewer.exe program and check if the building footprints that lasboundary has produced are what you expected:
By generating KML output instead and overlaying the result on Google Earth imagery you can easily check which and how many buildings were missed and go back to adjust the parameters for ‚lasclassify‘ in the previous tutorial to maybe achieve better results.
C:\lastools\bin>lasboundary -i tiles_final\fitch*.laz -merged ^ -keep_class 6 ^ -concavity 1.5 -disjoint ^ -o buildings.kml
The las2dem GUI specifies the -elevation switch in the generated command line, when the combobox for item is set to „elevation“. But, the help for las2dem does not mention the -elevation switch.
Hello Mike, the ‚-elevation‘ switch is the default of las2dem (and also of lasgrid and blast2dem). It may as well be omitted.
I am using lasboundary to extract building footprints, and output them as features in a shapefile, as per your tutorial. I need those features to have additional attributes, such as height. For my purposes, feature height can be an average of the heights of all the points that were used to create the feature’s footprint, or just the height of any of those points. Is there an existing ability of the lastools to calculate such a feature height?
The height of each point is stored as the Z coordinate in the SHP file. Maybe average all of them or take the 95th percentile as the building height?
Thanks, that is a good idea. I overlooked that the feature type is POLYGONZ.
However, when I crawl through the features of the shapefile using pyshp, I observe that none of the points of those features have a z coordinate.
Can you confirm that lasboundary is really intended to output the z coordinates for the points that it determines are part of a building footprint feature?
I confirm that lasboundary outputs the z coordinates for the points that are part of a building footprint feature. I wrote a different program in C++, using GDAL, to crawl the features of the resulting shapefile, and observed that the points have z coordinates.
It seems to be a bug in pyshp that prevents me from access the z coordinates when I crawl the shapefile with a Python program.
Happy to hear it works as intended … (-: