This is part one 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 such as DTMs, DSMs, building footprints, and iso-contours with multi-core batch processing.
To get started, download the latest version of LAStools and these 7 example flight strips: strip1, strip2, strip3, strip4, strip5, strip6, strip7 and put them into a folder called ‘.\lastools\bin\strips_raw’. For simplicity we will work inside the ‘.\lastools\bin’ directory, so open a DOS command line window and change directory to where this folder is on your computer, for example, with the command ‘cd c:\software\lastools\bin’ or with ‘cd d:\lastools\bin’ followed by ‘d:’.
I usually first run ‘lasview’ on newly received LiDAR to see if the data makes any sense. The command below will on-the-fly down-sample the data and display only around 5 million points:
C:\lastools\bin>lasview -i strips_raw\LDR*.laz -points 5000000
I usually inspect a few area-of-interests (AOIs) to get a few close-up looks at the data. If we do this several times for different areas it is best to first create a spatial indexing of the strips with ‘lasindex’ which it significantly speeds-up our spatial queries.
C:\lastools\bin>lasindex -i strips_raw\LDR*.laz -cores 3
Only add the ‘-cores 3’ option if you have a computer with at least 4 cores. If you have an 8 core machine you may even add ‘-cores 7’ and so on. Now start ‘lasview’ in the GUI mode either by double-clicking the executable and using the ‘browse …’ roll-out on the left side to load the seven flight strips or simply by adding ‘-gui’ to the command line.
C:\lastools\bin>lasview -i strips_raw\LDR*.laz -gui
Now you can maneuver up and down the file list to highlight the different strips (close the ‘browse …’ roll-out to see the header contents listed on the bottom left). You can again look at all the data by pressing the ‘VIEW’ button or you can inspect a smaller area by first picking a rectangular region as shown in the screenshot below.
In the GUI above you can see that the files have no projection information. We will add this later. The x, y, and z scaling factors are set to 0.001 which means that points are stored with millimeter resolution. Because airborne LiDAR is far from being that accurate we will later change the scaling factor to a more appropriate centimeter resolution. Finally, the creation day and year is missing, which – if we know when the data was flown – will be an easy fix.
After starting ‘lasview’ you can use the ‘c’ hot-key to toggle through the different ways to color the points. Another way to change the point coloring is by presetting it in the GUI as shown above or via the ‘color by …’ menu that opens when you right click.
The coloring by flight-line which is supposed to be stored in the ‘point_source_ID’ field of each point seems rather odd. This field was probably used for some other purpose and not cleaned up before the data was distributed. Not nice and we will fix that later. There are a few other visualization shown below. Coloring points by return turns single returns yellow, first of many returns red, and last of many returns blue. Usually there are also a few green points corresponding to intermediate (neither first nor last) of many returns but there are none in this data set as the ‘lasinfo’ report will confirm later. Toggle through the colorings with the ‘c’ key until it displays the intensities. You may need to press the ‘=’ key a few times to increase the point size. Next, press the ‘t’ key to triangulate the points and then decrease the point size with the ‘-‘ key until they have disappeared and you can see a hillside shading of the TIN. Use the ‘h’ key to toggle through different shadings for the TIN triangles.
Before attempting to improve these LAZ files and prepare them for processing we want to make sure that our work will be worthwhile by running a quick visualization based of how well the flight strips fit together. We do this with the ‘lasoverlap’ tool.
C:\lastools\bin>lasoverlap -i strips_raw\LDR*.laz ^ -files_are_flightlines ^ -utm 19N ^ -step 1.0 -max_diff 1.0 ^ -o strip_overlap.png
I happen to know that the missing projection in the LAZ files is the UTM zone 19N. I can simply add this to the command line with ‘-utm 19N’. Then ‘lasoverlap’ will produce in addition to the PNG output also a KML file that allows visualizing the overlap and difference rasters within the geo-spatial context of Google Earth.
Serious troubles in the data would be evidenced by void areas in the overlap raster that are obviously caused by poor flight planning (not water bodies) and by deeply blue or deeply red saturated areas in open (=> non-forested) terrain in the difference raster. If either was the case you should send the data back to the vendor to have him fix it … (-:
Now we compare the LiDAR elevations with a set of 29 exactly known control points that were obtained through a ground survey that are listed below. If you want to follow this exercise along you should copy the text below and store it to a file called ‘strips_raw\cps.csv’.
C:\lastools\bin>more strips_raw\cps.csv Type,Easting,Northing,Z Open/Paved,273299.68,4715133.88,74.17 Open/Paved,273477.61,4714979.29,73.85 Open/Paved,274001.2,4714540.29,72.5 Open/Paved,273670.13,4714817.4,73.34 Open/Paved,273677.42,4715018.66,74.26 Open/Paved,273400.1,4714528.98,73.15 Open/Paved,274511.59,4714905.97,95.7 Open/Paved,275074.66,4714841.98,120.18 Open/Paved,275409.65,4714994.76,113.41 Field,273321.18,4714946.83,73.46 Field,273601.49,4715101.78,74.35 Field,273646.76,4714972.94,73.97 Field,273890.5,4714457.59,71.13 Field,274650.24,4714903.44,105.13 Field,274522.36,4714829.74,98.59 Field,275474.47,4714780.03,127.63 Field,275636.39,4714868.85,120.72 Field,274747.37,4714932.57,116.11 Field,272795.36,4714503.86,126.9 Forested,272547.02,4714623.09,127.71 Forested,273205.33,4714900.27,79.36 Forested,272530.52,4715045.46,113.48 Forested,275237.48,4715049.57,120.31 Forested,275268.37,4714543.82,104.99 Forested,274666.09,4714497.49,108.86 Forested,274403.56,4715053.43,93.17 Forested,274901.6,4714493.63,114.64 Forested,274658.37,4715072.74,104.49 Forested,274121.73,4714524.51,72.06
The first entry describes the location of the control point and the following three entries are its x, y, and z coordinate. We now compare the seven LiDAR strips against these 29 control points. You can also do thus via the GUI as shown below but you need to make sure to that the ‘keep ground (2) ‘ and the ‘keep keypoint (8)’ check-buttons are not checked since our LiDAR data is not yet classified.
C:\lastools\bin>lascontrol -i strips_raw\*.laz ^ -cp strips_raw\cps.csv ^ -parse sxyz diff,lidar_z,Type,Easting,Northing,Z 0.1451,74.3151,Open/Paved,273299.68,4715133.88,74.17 0.0190857,73.8691,Open/Paved,273477.61,4714979.29,73.85 0.08872,72.5887,Open/Paved,274001.2,4714540.29,72.5 -0.00419681,73.3358,Open/Paved,273670.13,4714817.4,73.34 0.054868,74.3149,Open/Paved,273677.42,4715018.66,74.26 0.0318439,73.1818,Open/Paved,273400.1,4714528.98,73.15 -0.0720668,95.6279,Open/Paved,274511.59,4714905.97,95.7 -0.0670072,120.113,Open/Paved,275074.66,4714841.98,120.18 0.0756029,113.486,Open/Paved,275409.65,4714994.76,113.41 0.00132683,73.4613,Field,273321.18,4714946.83,73.46 0.0750162,74.425,Field,273601.49,4715101.78,74.35 0.00228472,73.9723,Field,273646.76,4714972.94,73.97 0.166194,71.2962,Field,273890.5,4714457.59,71.13 -0.392266,104.738,Field,274650.24,4714903.44,105.13 -0.0705893,98.5194,Field,274522.36,4714829.74,98.59 -0.423622,127.206,Field,275474.47,4714780.03,127.63 0.217528,120.938,Field,275636.39,4714868.85,120.72 -0.12489,115.985,Field,274747.37,4714932.57,116.11 0.0195411,126.92,Field,272795.36,4714503.86,126.9 17.0342,144.744,Forested,272547.02,4714623.09,127.71 9.20693,88.5669,Forested,273205.33,4714900.27,79.36 26.355,139.835,Forested,272530.52,4715045.46,113.48 20.2143,140.524,Forested,275237.48,4715049.57,120.31 14.5363,119.526,Forested,275268.37,4714543.82,104.99 0.173861,109.034,Forested,274666.09,4714497.49,108.86 20.5742,113.744,Forested,274403.56,4715053.43,93.17 0.00432622,114.644,Forested,274901.6,4714493.63,114.64 23.8816,128.372,Forested,274658.37,4715072.74,104.49 0.3645,72.4245,Forested,274121.73,4714524.51,72.06 sampled TIN at 29 of 29 control points. avgabs/rms/stddev/avg of elevation errors are 4.63438/9.61987/8.62324/4.55475 meter. skew is 1.47593.
The output of ‘lascontrol’ could be directed into a file using the ‘-o control.txt’ option. The tool prefixes two numbers to the beginning of every line: the elevation value that was computed from the LiDAR data at the xy location of the control point and the difference between the computed elevation and the elevation of the control point. At the end a summary reports the average absolute error, the relative mean-square error, the standard deviation, and the average error across all the differences. Of course, for control points of type ‘Forested’ we have terrible results because we use all LiDAR points in this computation and not just those that correspond to ground returns. Therefore we will have to repeat this exercise at a later point once we have ground-classified the LiDAR. Via the GUI of ‘lascontrol’ you can zoom in on different areas-of-interests and visually inspect how well the data matches the control points (see below).
Finally let us run ‘lasinfo’ on one one of the strips and also compute the density:
C:\lastools\bin>lasinfo -i strips_raw\LDR030828_211804_0.laz -cd reporting all LAS header entries: file signature: 'LASF' file source ID: 0 global_encoding: 0 project ID GUID data 1-4: 0 0 0 '' version major.minor: 1.0 system identifier: 'ALSXX' generating software: 'ALSXX_PP V2.56c' file creation day/year: 0/0 header size: 227 offset to point data: 5587 number var. length records: 3 point data format: 0 point data record length: 20 number of point records: 2404613 number of points by return: 2210130 0 194483 0 0 scale factor x y z: 0.001 0.001 0.001 offset x y z: 0 4000000 0 min x y z: 272254.830 4714389.375 65.050 max x y z: 275391.197 4714711.445 169.208 variable length header record 1 of 3: reserved 43707 user ID 'LeicaGeo' record ID 1001 length after header 5120 description '' variable length header record 2 of 3: reserved 43707 user ID 'LeicaGeo' record ID 1002 length after header 22 description 'MissionInfo' variable length header record 3 of 3: reserved 43707 user ID 'LeicaGeo' record ID 1003 length after header 54 description 'UserInputs' the header is followed by 2 user-defined bytes LASzip compression (version 2.0r2 c2 50000): POINT10 2 reporting minimum and maximum for all LAS point record entries ... X 272254812 275391197 Y 714389375 714711445 Z 65050 169208 intensity 0 255 edge_of_flight_line 0 1 scan_direction_flag 0 1 number_of_returns_of_given_pulse 1 2 return_number 1 3 classification 1 1 scan_angle_rank -13 13 user_data 0 0 point_source_ID 0 255 WARNING: 2 points outside of header bounding box number of last returns: 2209098 covered area in square units/kilounits: 728996/0.73 point density: all returns 3.30 last only 3.03 (per square units) spacing: all returns 0.55 last only 0.57 (in units) overview over number of returns of given pulse: 2014615 389998 0 0 0 0 0 histogram of classification of points: 2404613 Unclassified (1) real max x larger than header max x by 0.000422 real min x smaller than header min x by 0.017715