Tutorial: LiDAR preparation

This is part two 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 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:‘. It may be helpful (not mandatory) to follow the tutorial on quality checking first.

Create a new folder called ‚.\LAStools\bin\tiles_raw‘ for storing the LiDAR tiles. Then run ‚lastile.exe‘ in GUI mode and load the 7 example flight strips. You can either do this by double-clicking ‚lastile.exe‘ and using the ‚browse …‘ menu or by entering the command below:

C:\LAStools\bin>lastile -i strips_raw\LDR*.laz -gui

Set the GUI options as shown below to create a buffered tiling from the original flight strips. Check the button ‚files are flightlines‘ so that points from different flight lines get a unique flight line ID stored in their ‚point source ID‘ attribute. This makes it possible to identify later which points of a tile came from the same flight strip. Use the default 1000 to specify the tile size and request a buffer of 50 meters around every tile. This buffer helps to reduce edge artifacts at the tile boundaries in a tile-based processing pipeline.  Shift the tiling grid by 920 in x direction and by 320 in y direction so the flight strips fit in exactly 4 tiles. This is not mandatory but results in fewer tiles and therefore fewer cuts (experimentally discovered). Requests compressed output to overcome the I/O bottleneck. Set the projection to UTM zone 19 north and the vertical datum to NAVD88. When projection information is missing in the original flight strips it is a good idea to add this as early as possible so it is available in subsequent processing steps.

Press ‚RUN‘ and the ‚RUN‘ window will pop up. 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 almost equivalent command here. It has been enhanced by one additional parameter shown in red. It quantizes the input on-the-fly to a more appropriate centimeter [cm] resolution because airborne LiDAR does not have the millimeter [mm] resolution that the raw flight strips are stored with. See the lasinfo report shown at the end of the tutorial on quality checking to see the excessive millimeter resolution that the additional parameter in red is fixing.

C:\LAStools\bin>lastile -i strips_raw\LDR*.laz ^
                        -files_are_flightlines ^
                        -rescale 0.01 0.01 0.01 ^
                        -utm 19N -vertical_navd88 ^
                        -tile_size 1000 -buffer 50 ^
                        -tile_ll 920 320 ^
                        -odir tiles_raw -o fitch.laz

Create a new folder called ‚.\LAStools\bin\tiles_ground‘ for storing the ground-classified LiDAR tiles. Then run ‚lasground.exe‘ in GUI mode and load the 4 tiles. You can do this by double-clicking ‚lasground.exe‘ and using the ‚browse …‘ menu or by entering the command below:

C:\LAStools\bin>lasground -i tiles_raw\fitch*.laz -gui

Set the GUI options as shown below to classify the bare-earth points in all tiles. Select the ‚metro‘ setting (which uses a step size of 50m) because there are very large hangars in the point cloud and the ‚ultra fine‘ setting for the initial ground estimate. For more details on these parameters read the corresponding README.txt file. Specify the output directory as ‚tiles_ground‘ and select compressed LAZ output. Select 4 cores only if you have a computer with that many cores. This will assign tiles to different cores to process multiple tiles in parallel.

Press ‚RUN‘ and the ‚RUN‘ window will pop up. 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>lasground -i tiles_raw\fitch*.laz ^
                          -metro -ultra_fine ^
                          -odir tiles_ground -olaz ^
                          -cores 4

Run ‚lasview.exe‘ on ‚fitch_273920_4714320.laz‘ in the ‚tiles_ground‘ folder to inspect the result of your bare-earth classification. Press ‚g‘ on the keyboard to show only the ground points. After all points have loaded press ‚t‘ on the keyboard to triangulate only the points shown. Press ‚-‚ to make the points disappear and press ‚h‘ multiple times to iterate over the possible ways to shade the triangulation. Now press ‚=‘ to make the points reappear, press ‚u‘ to show the unclassified (non-ground) points, and press ‚c‘ a few times to iterate over the possible ways to color the points.

There are „dirt“ points below the ground and „air“ points far above the ground such as the one pointed at by the mouse cursor. We remove them with ‚lasheight.exe‘. This tool computes for each point the vertical distance to the triangulation of the ground points and stores it in the ‚user_data‘ field. We also need these heights for the following step of finding buildings and vegetation. Create a new folder called ‚.\LAStools\bin\tiles_height‘ for storing the cleaned LiDAR tiles with height above ground information. Run ‚lasheight.exe‘ in GUI mode and load the 4 ground-classified tiles by double-clicking ‚lasheight.exe‘ and using the ‚browse …‘ menu or by entering the command below:

C:\LAStools\bin>lasheight -i tiles_ground\fitch*.laz -gui

Configure the GUI with the settings shown below. By default ‚lasheight.exe‘ uses the points classified as ground to construct a TIN and then calculates the height of all other points in respect to this ground surface. With ‚drop_above 30‘ and  ‚drop_below -2‘ all points that are 30 meters above or 2 meters below the ground are removed from the compressed LAZ tiles output to the ‚tiles_height‘ folder. Run the process on multiple cores if possible.

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>lasheight -i tiles_ground\fitch*.laz ^
                          -drop_below -2 -drop_above 30 ^
                          -odir tiles_height -olaz ^
                          -cores 4

Next, create a new folder called ‚.\LAStools\bin\tiles_classified‘ for storing classified LiDAR tiles with building and vegetation points. Run ‚lasclassify.exe‘ in GUI mode and load the 4 height-cleaned tiles by double-clicking or by entering the command below:

C:\LAStools\bin>lasclassify -i tiles_height\fitch*.laz -gui

Set the GUI settings as shown below to identify buildings and trees. Mileage may vary on this step because automatic LiDAR classification is a hard problem. The default settings require a density of at least 2 pulses per square. Setting the step size for computing planarity (or ruggedness) to 3 meters has fewer false positives but misses some smaller buildings. For more details on tuning these parameters see the corresponding README.txt file.

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 command-line:

C:\LAStools\bin>lasclassify -i tiles_height\fitch*.laz ^
                            -step 3 ^
                            -odir tiles_classified -olaz ^
                            -cores 4

Run ‚lasview.exe‘ on ‚fitch_273920_4714320.laz‘ in the ‚tiles_classified‘ folder and be amazed about the result of your building and vegetation classification. Although not 100 percent correct, ‚lasclassify.exe‘ identifies all man-made structures as class 6 and marks all vegetation above two meters as class 5.

Create the last folder called ‚.\LAStools\bin\tiles_final‘ for storing the final LiDAR tiles without the 50 meter buffers that are ready to be shipped to the customer. Run ‚lastile.exe‘ in GUI mode and load the 4 classified tiles with:

C:\LAStools\bin>lastile -i tiles_classified\fitch*.laz -gui

Set the GUI settings as shown below to remove the buffers that were added in the very beginning and carried through all processing steps to avoid artifacts along the boundaries.

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 command-line that has one additional parameter shown in red. This sets the user data field of each point back to zero that was set by ‚lasheight.exe‘ to the height above ground (in decimeters). This number is meaningless to the customer and setting it to zero improves compression.

C:\LAStools\bin>lastile -i tiles_classified\fitch*.laz ^
                        -set_user_data 0 ^
                        -remove_buffer ^
                        -odir tiles_final -olaz

Finally let us run ‚lasinfo‘ on one of the generated tiles and also compute the density:

C:\LAStools\bin>lasinfo -i tiles_final\fitch_271920_4714320.laz -cd
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            0
  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000
  version major.minor:        1.0
  system identifier:          'LAStools (c) by Martin Isenburg'
  generating software:        'lastile (131011) commercial'
  file creation day/year:     0/0
  header size:                227
  offset to point data:       5689
  number var. length records: 4
  point data format:          0
  point data record length:   20
  number of point records:    2571486
  number of points by return: 2196212 0 375274 0 0
  scale factor x y z:         0.01 0.01 0.01
  offset x y z:               0 4000000 0
  min x y z:                  272044.80 4714466.57 82.34
  max x y z:                  272919.99 4715243.91 169.21
variable length header record 1 of 4:
  reserved             43707
  user ID              'LeicaGeo'
  record ID            1001
  length after header  5120
  description          ''
variable length header record 2 of 4:
  reserved             43707
  user ID              'LeicaGeo'
  record ID            1002
  length after header  22
  description          'MissionInfo'
variable length header record 3 of 4:
  reserved             43707
  user ID              'LeicaGeo'
  record ID            1003
  length after header  54
  description          'UserInputs'
variable length header record 4 of 4:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34735
  length after header  48
  description          'by LAStools of Martin Isenburg'
    GeoKeyDirectoryTag version 1.1.0 number of keys 5
      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
      key 3072 tiff_tag_location 0 count 1 value_offset 32619 - ProjectedCSTypeGeoKey: PCS_WGS84_UTM_zone_19N
      key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
      key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter
      key 4096 tiff_tag_location 0 count 1 value_offset 5103 - VerticalCSTypeGeoKey: VertCS_North_American_Vertical_Datum_1988
the header is followed by 2 user-defined bytes
LASzip compression (version 2.2r0 c2 50000): POINT10 2
LAStiling (idx 8, lvl 2, sub 0, bbox 271920 4.71232e+006 275920 4.71632e+006)
reporting minimum and maximum for all LAS point record entries ...
  X   27204480   27291999
  Y   71446657   71524391
  Z       8234      16921
  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     6
  scan_angle_rank   -13    21
  user_data           0     0
  point_source_ID     1     7
number of last returns: 2196301
covered area in square meters/kilometers: 408700/0.41
point density: all returns 6.29 last only 5.37 (per square meter)
      spacing: all returns 0.40 last only 0.43 (in meters)
overview over number of returns of given pulse: 1821027 750459 0 0 0 0 0
histogram of classification of points:
          286793  Unclassified (1)
          589019  Ground (2)
         1586840  High Vegetation (5)
          108834  Building (6)