Smooth DTM from Drone LiDAR off Velodyne HDL 32A mounted on DJI M600 UAV

Recently we attempted to do a small LiDAR survey by drone for a pet project of our CEO in our “code and surf camp” here in Samara, Costa Rica. But surveying is difficult when you are a novice and we ran into a trajectory issue. The dramatic “wobbles” were entirely our fault, but fortunately our mistakes also led to something useful: We found some LAS export bugs. Our laser scanner was a Velodyne HDL-32E integrated with a NovAtel INS into the Snoopy Series A HD made by LiDARUSA. The system was carried by a DJI Matrice 600 (M600) drone. We processed the trajectory with NovAtel Inertial Explorer (here we made the “wobbles” error) and finally exported the LAS and LAZ files with ScanLook PC (version 1.0.182) from LiDARUSA.

While we were investigating our “wobbles” (which clearly were our mistake) we also found five different LAS export bugs in ScanLook PC that seem to have started sometime after version 1.0.171 and will likely end with version 1.0.193. Below an illustration of a correct export from version 1.0.129 and a buggy export from version 1.0.182. In both instances you see the returns from one revolution of the Velodyne HDL-32E scanner head ordered by their GPS time stamps and colored to distinguish the 32 separate beams. In the buggy version, groups of around seven non-adjacent returns are given the same time stamp. This bug will only affect you, if correct GPS time stamps are important for your subsequent LiDAR processing or if your client explicitly asked for ASPRS specification compliant LAS files. We plan to publish another blog post detailing how to find this GPS time stamping bug (and the other four bugs we found).

During the many interactions we had working through “wobbles” and export bugs, we obtained a nice set of six flight lines from Seth Gulich of Bowman Consulting – a US American company based in Stuart, Florida – who flew an identical “Snoopy Series A HD” system also on a DJI Matrice 600 drone at approximately 100 feet above ground level above a model airplane airport in Palm Beach, Florida. You can download the data set here. In the following we will check the flight line alignment of this data set and then process it into a smooth DTM. All command lines used are summarized in this text file.

First we generate a lasinfo report that includes a number of histograms for on-the-fly merged flight lines with lasinfo and then use the z coordinate histogram from the lasinfo report to set reasonable min/max values for the elevation color ramp of lasview:

lasinfo -i 0_strips_raw\Velodyne*.laz -merged ^
        -cd ^
        -histo z 1 ^
        -histo user_data 1 ^
        -histo point_source 1 ^
        -o 1_quality\Velodyne_merged_info.txt

lasview -i 0_strips_raw\Velodyne*.laz ^
        -points 10000000 ^
        -set_min_max 25 75

The lasinfo report shows no information about the coordinate reference system. We found out experimentally that the horizontal coordinates seem to be EPSG code 2236 and that the vertical units are most likely be US survey feet. The warnings you will see in the lasinfo report have to do with the fact that the double-precision bounding box stored in the LAS header was populated with numbers that have many more decimal digits than the coordinates in the file, which only have millifeet resolution as all three scale factors are 0.001 (meaning coordinates have three decimal digits). The information which of the 32 lasers was collecting which point is stored in both the ‘user data’ and the ‘point source ID’ field which is evident from the histograms in the lasinfo report. We need to be careful not to override both fields in later processing.

Next we use lasoverlap to check how well the LiDAR points from the flight out and the flight back align vertically. This tool computes the difference of the lowest points for each square foot covered by multiple flight lines. Differences of less than a quarter of a foot are both times mapped to white, differences of more than one foot (more than half a foot) are mapped to saturated red or blue depending on whether the difference is positive or negative in the first run (in the second run):

lasoverlap -i 0_strips_raw\Velodyne*.laz ^
           -faf ^
           -min_diff 0.25 -max_diff 1.00 -step 1 ^
           -odir 1_quality -o overlap_025_100.png

lasoverlap -i 0_strips_raw\Velodyne*.laz ^
           -faf ^
           -min_diff 0.25 -max_diff 0.50 -step 1 ^
           -odir 1_quality -o overlap_025_050.png

We use a new feature of the LAStools GUI (as of version 180429) to closer inspect large red or blue areas. With lasmerge we clip out regions that looks suspect for closer examination with lasview. First we spatially index the flight lines to make this process faster. With the ‘-gui’ switch we start the tool in GUI mode with flight lines already loaded. Using the new PNG overlay roll-out on the left we add the ‘overlap_025_050_diff.png’ image from the quality folder created in the last step and clip out three areas.

lasindex -i 0_strips_raw\Velodyne*.laz
         -tile_size 10 -maximum -100 ^
         -cores 3

lasmerge -i 0_strips_raw\Velodyne*.laz -gui

You can also clip out these three areas using the command lines below:

lasmerge -i 0_strips_raw\Velodyne*.laz ^
         -faf ^
         -inside_tile 939500 889860 100 ^
         -o 1_quality\939500_889860.laz

lasmerge -i 0_strips_raw\Velodyne*.laz ^
         -faf ^
         -inside_tile 940400 889620 100 ^
         -o 1_quality\940400_889620.laz

lasmerge -i 0_strips_raw\Velodyne*.laz ^
         -faf ^
         -inside_tile 940500 890180 100 ^
         -o 1_quality\940500_890180.laz

The reader may inspect the areas 939500_889860.laz, 940400_889620.laz, and 940500_890180.laz with lasview using profile views via hot keys ‘x’ and switching back and forth between the points from different flight lines via hot keys ‘0’, ‘1’, ‘2’, ‘3’, … for individual and ‘a’ for all flight lines as we have done it in previous tutorials [1,2,3]. Using drop-lines or rise-lines via the pop-up menu gives you a sense of scale. Removing points with lastrack that are horizontally too far from the trajectory could be one strategy to use fewer outliers. But as our surfaces are expected to be “fluffy” (because we have a Velodyne LiDAR system), we accept these flight line differences and continue processing.

Here the complete LAStools processing pipeline for creating an average ground model from the set of six flight lines that results in the hillshaded DTM shown below. The workflow is similar to those we have developed in earlier blog posts for Velodyne Puck based systems like the Hovermap and the Yellowscan and in the other Snoopy tutorial. All command lines used are summarized in this text file.

Hillshaded DTM with half foot resolution generated via average ground computation with LAStools.

In the first step we lastile the six flight lines into 250 by 250 feet tiles with 25 feet buffer while preserving flight line information. The flight line information will be stored in the “point source ID” field of each point and therefore override the beam ID that is currently stored there. But the beam ID is also stored in the “user data” field as the  lasinfo report had told us. We set all classifications to zero and add information about the horizontal coordinate reference system EPSG code 2236 and the vertical units (US Survey Feet).

lastile -i 0_strips_raw\*.laz ^
        -faf ^
        -set_classification 0 ^
        -epsg 2236 -elevation_survey_feet ^
        -tile_size 250 -buffer 25 -flag_as_withheld ^
        -odir 2_tiles_raw -o pb.laz

On three cores in parallel we then lassort the points in the tiles into a space-filling curve order which will accelerate later operations.

lassort -i 2_tiles_raw\*.laz ^
        -odir 2_tiles_sorted -olaz ^
        -cores 3

Next we use lasthin to classify the point whose elevation is closest to the 5th elevation percentile among all points falling into its cell with classification code 8. We run lasthin multiple times and each time increase the cell size from 1, 2, 4, 8 to 16 foot. We do this because we have requested the 5th elevation percentile to only be computed when there are at least 20 points in the cell. Percentiles are statistical measures and need a reasonable sample size to be stable. Because drone flights are very dense in the center and more sparse at the edges this increase in cell size assures that we have a good selection of points classified with classification code 8 across the entire survey area.

lasthin -i 2_tiles_sorted\*.laz ^
        -step 1 -percentile 5 20 -classify_as 8 ^
        -odir 3_tiles_thinned_p05_step01 -olaz ^
        -cores 3

lasthin -i 3_tiles_thinned_p05_step01\*.laz ^
        -step 2 -percentile 5 20 -classify_as 8 ^
        -odir 3_tiles_thinned_p05_step02 -olaz ^
        -cores 3

lasthin -i 3_tiles_thinned_p05_step02\*.laz ^
        -step 4 -percentile 5 20 -classify_as 8 ^
        -odir 3_tiles_thinned_p05_step04 -olaz ^
        -cores 3

lasthin -i 3_tiles_thinned_p05_step04\*.laz ^
        -step 8 -percentile 5 20 -classify_as 8 ^
        -odir 3_tiles_thinned_p05_step08 -olaz ^
        -cores 3

lasthin -i 3_tiles_thinned_p05_step08\*.laz ^
        -step 16 -percentile 5 20 -classify_as 8 ^
        -odir 3_tiles_thinned_p05_step16 -olaz ^
        -cores 3

Then we let lasground_new run on only the points classified with classification code 8 (i.e. by ignoring the points still classified with code 0) which classifies them into ground (code 2) and non-ground (code 1).

lasground_new -i 3_tiles_thinned_p05_step16\*.laz ^
              -ignore_class 0 ^
              -town ^
              -odir 4_tiles_ground_low -olaz ^
              -cores 3

The ground points we have computed form somewhat of a lower envelope of the “fluffy” points of a Velodyne scanner. With lasheight we now draw all the points near the ground – namely those from 0.1 foot below to 0.4 foot above the ground – into a new classification code 6 that we term “thick ground”. The ‘-do_not_store_in_user_data’ switch prevent the default behavior of lasheight from happening, which would override the beam ID information that it stored in the ‘user data’ field with approximate height value.

lasheight -i 4_tiles_ground_low\*.laz ^
          -classify_between -0.1 0.4 6 ^
          -do_not_store_in_user_data ^
          -odir 4_tiles_ground_thick -olaz ^
          -cores 3

A few close-up shots of the resulting “thick ground” are shown in the picture gallery below.

We then use lasgrid to average the (orange) thick ground points onto a regular grid with a cell spacing of half a foot. We do not grid the tile buffers by adding the ‘-use_tile_bb’ switch.

lasgrid -i 4_tiles_ground_thick\*.laz ^
        -keep_class 6 ^
        -step 0.5 -average ^
        -use_tile_bb ^
        -odir 5_tiles_gridded_mean_ground -olaz &
        -cores 3

Finally we use blast2dem to merge all the averaged ground point grids into one file, interpolate across open areas without ground points, and compute the hillshaded DTM shown above. All command lines used are summarized in this text file.

blast2dem -i 5_tiles_gridded_mean_ground\*.laz ^
          -merged ^
          -step 0.5 ^
          -hillshade ^
          -o dtm.png

We thank Seth Gulich of Bowman Consulting for sharing this LiDAR data set with us. It was flown with a DJI Matrice 600 drone carrying a “Snoopy A series HD” LiDAR system from LidarUSA.

New Step-by-Step Tutorial for Velodyne Drone LiDAR from Snoopy by LidarUSA

The folks from Harris Aerial gave us LiDAR data from a test-flight of one of their drones, the Carrier H4 Hybrid HE (with a 5kg maximum payload and a retail price of US$ 28,000), carrying a Snoopy A series LiDAR system from LidarUSA in the countryside near Huntsville, Alabama. The laser scanner used by the Snoopy A series is a Velodyne HDL 32E that has 32 different laser/detector pairs that fire in succession to scan up to 700,000 points per second within a range of 1 to 70 meters. You can download the raw LiDAR file from the 80 second test flight here. As always, the first thing we do is to visualize the file with lasview and to generate a textual report of its contents with lasinfo.

lasview -i Velodyne001.laz -set_min_max 680 750

It becomes obvious that the drone must have scanned parts of itself (probably the landing gear) during the flight and we exploit this fact in the later processing. The information which of the 32 lasers was collecting which point is stored into the ‘point source ID’ field which is usually used for the flightline information. This results in a psychedelic look in lasview as those 32 different numbers get mapped to the 8 different colors that lasview uses for distinguishing flightlines.

The lasinfo report we generate computes the average point density with ‘-cd’ and includes histograms for a number of point attributes, namely for ‘user data’, ‘intensity’, ‘point source ID’, ‘GPS time’, and ‘scan angle rank’.

lasinfo -i Velodyne001.laz ^
        -cd ^
        -histo user_data 1 ^
        -histo point_source 1 ^
        -histo intensity 16 ^
        -histo gps_time 1 ^
        -histo scan_angle_rank 5 ^
        -odir quality -odix _info -otxt

You can download the resulting report here and it will tell you that the information which of the 32 lasers was collecting which point was stored both into the ‘user data’ field and into the ‘point source ID’ field. The warnings you see below have to do with the fact that the double-precision bounding box stored in the LAS header was populated with numbers that have many more decimal digits than the coordinates in the file, which only have millimeter (or millifeet) resolution as all three scale factors are 0.001 (meaning three decimal digits).

WARNING: stored resolution of min_x not compatible with x_offset and x_scale_factor: 2171988.6475160527
WARNING: stored resolution of min_y not compatible with y_offset and y_scale_factor: 1622812.606925504
WARNING: stored resolution of min_z not compatible with z_offset and z_scale_factor: 666.63504345017589
WARNING: stored resolution of max_x not compatible with x_offset and x_scale_factor: 2172938.973065129
WARNING: stored resolution of max_y not compatible with y_offset and y_scale_factor: 1623607.5209975131
WARNING: stored resolution of max_z not compatible with z_offset and z_scale_factor: 1053.092674726669

Both the “return number” and the “number of returns” attribute of every points in the file is 2. This makes it appear as if the file would only contain the last returns of those laser shots that produced two returns. However, as the Velodyne HDL 32E only produces one return per shot we can safely conclude that those numbers should all be 1 instead of 2 and that this is just a small bug in the export software. We can easily fix this with las2las.

reporting minimum and maximum for all LAS point record entries ...
[...]
 return_number 2 2
 number_of_returns 2 2
[...]

The lasinfo report lacks information about the coordinate reference system as there is no VLR that stores projection information. Harris Aerial could not help us other than telling us that the scan was near Huntsville, Alamaba. Measuring certain distances in the scene like the height of the house or the tree suggests that both horizontal and vertical units are in feet, or rather in US survey feet. After some experimenting we find that using EPSG 26930 for NAD83 Alabama West but forcing the default horizontal units from meters to US survey feet gives a result that aligns well with high-resolution Google Earth imagery as you can see below:

lasgrid -i flightline1.laz ^
        -i flightline2.laz ^
        -merged ^
        -epsg 26930 -survey_feet ^
        -step 1 -highest ^
        -false -set_min_max 680 750 ^
        -o testing26930usft.png

Using EPSG code 26930 but with US survey feet instead of meters results in nice alignment with GE imagery.

We use the fact that the drone has scanned itself to extract an (approximate) trajectory by isolating those LiDAR returns that have hit the drone. Via a visual check with lasview (by hovering with the cursor over the lowest drone hits and pressing hotkey ‘i’) we determine that the lowest drone hits are all above 719 feet. We use two calls to las2las to split the point cloud vertically. In the same call we also change the resolution from three to two decimal digits, fix the return number issue, and add the missing geo-referencing information:

las2las -i Velodyne001.laz ^
        -rescale 0.01 0.01 0.01 ^
        -epsg 26930 -survey_feet -elevation_survey_feet ^
        -set_return_number 1 ^
        -set_number_of_returns 1 ^
        -keep_z_above 719 ^
        -odix _above719 -olaz

las2las -i Velodyne001.laz ^
        -rescale 0.01 0.01 0.01 ^
        -epsg 26930 -survey_feet -elevation_survey_feet ^
        -set_return_number 1 ^
        -set_number_of_returns 1 ^
        -keep_z_below 719 ^
        -odix _below719 -olaz

We then use the manual editing capabilities of lasview to change the classifications of the LiDAR points that correspond to drone hits from 1 to 12, which is illustrated by the series of screen shots below.

lasview -i Velodyne001_above719.laz

The workflow illustrated above results in a tiny LAY file that is part of the LASlayers functionality of LAStools. It only encodes the few changes in classifications that we’ve made to the LAZ file without re-writing those parts that have not changed. Those interested may use laslayers to inspect the structure of the LAY file:

laslayers -i Velodyne001_above719.laz

We can apply the LAY file on-the-fly with the ‘-ilay’ option, for example, when running lasview:

lasview -i Velodyne001_above719.laz -ilay

To separate the drone-hit trajectory from the remaining points we run lassplit with the ‘-ilay’ option and request to split by classification with this command line:

lassplit -i Velodyne001_above719.laz -ilay ^
         -by_classification -digits 3 ^
         -olaz

This gives us two new files with the three-digit appendices ‘_001’ and ‘_012’. The latter one contains those points we marked as being part of the trajectory. We now want to use lasview to – visually – find a good moment in time where to split the trajectory into multiple flightlines. The lasinfo report tells us that the GPS time stamps are in the range from 418,519 to 418,602. In order to use the same trick as we did in our recent article on processing LiDAR data from the Hovermap Drone, where we mapped the GPS time to the intensity for querying it via lasview, we first need to subtract a large number from the GPS time stamps to bring them into a suitable range that fits the intensity field as done here.

lasview -i Velodyne001_above719_012.laz ^
        -translate_gps_time -418000 ^
        -bin_gps_time_into_intensity 1
        -steps 5000

The ‘-steps 5000’ argument makes for a slower playback (press ‘p’ to repeat) to better follow the trajectory.

Hovering with the mouse over a point that – visually – seems to be one of the turning points of the drone and pressing ‘i’ on the keyboard shows an intensity value of 548 which corresponds to the GPS time stamp 418548, which we then use to split the LiDAR point cloud (without the trajectory) into two flightlines:

las2las -i Velodyne001_below719.laz ^
        -i Velodyne001_above719_001.laz ^
        -merged ^
        -keep_gps_time_below 418548 ^
        -o flightline1.laz

las2las -i Velodyne001_below719.laz ^
        -i Velodyne001_above719_001.laz ^
        -merged ^
        -keep_gps_time_above 418548 ^
        -o flightline2.laz

Next we use lasoverlap to check how well the LiDAR points from the flight out and the flight back align vertically. This tool computes the difference of the lowest points for each square foot covered by both flightlines. Differences of less than a quarter of a foot are mapped to white, differences of more than half a foot are mapped to saturated red or blue depending on whether the difference is positive or negative:

lasoverlap -i flightline1.laz ^
           -i flightline2.laz ^
           -faf ^
           -min_diff 0.25 -max_diff 0.50 -step 1 ^
           -odir quality -o overlap_025_050.png

We then use a new feature of the LAStools GUI (as of version 180429) to closer inspect larger red or blue areas. We want to use lasmerge and clip out any region that looks suspect for closer examination with lasview. We start the tool in the GUI mode with the ‘-gui’ command and the two flightlines pre-loaded. Using the new PNG overlay roll-out on the left we add the ‘overlap_025_050_diff.png’ image from the quality folder created in the last step and clip out three areas.

lasmerge -i flightline1.laz -i flightline2.laz -gui

You can also clip out these three areas using the command lines below:

lasmerge -i flightline1.laz -i flightline2.laz ^
         -faf ^
         -inside 2172500 1623160 2172600 1623165 ^
         -o clip2500_3160_100x005.laz

lasmerge -i flightline1.laz -i flightline2.laz ^
         -faf ^
         -inside 2172450 1623450 2172550 1623455 ^
         -o clip2450_3450_100x005.laz

lasmerge -i flightline1.laz -i flightline2.laz ^
         -faf ^
         -inside 2172430 1623290 2172530 1623310 ^
         -o clip2430_3290_100x020.laz

A closer inspection of the three cut out slices explains the red and blue areas in the difference image created by lasoverlap. We find a small systematic error in two of the slices. In slice ‘clip2500_3160_100x005.laz‘ the green points from flightline 1 are on average slightly higher than the red points from flightline 2. Vice-versa in slice ‘clip2450_3450_100x005.laz‘ the green points from flightline 1 are on average slightly lower than the red points from flightline 2. However, the error is less than half a foot and only happens near the edges of the flightlines. Given that our surfaces are expected to be “fluffy” anyways (as is typical for Velodyne LiDAR systems), we accept these differences and continue processing.

The strong red and blue colors in the center of the difference image created by lasoverlap is easily explained by looking at slice ‘clip2430_3290_100x020.laz‘. The scanner was “looking” under a gazebo-like open roof structure from two different directions and therefore always seeing parts of the floor in one flightline that were obscured by the roof in the other.

While working with this data we’ve also implemented a new feature for lastrack that computes the 3D distance between LiDAR points and the trajectory and allows storing the result as an additional per point attribute with extra bytes. Those can then be visualized with lasgrid. Here an example:

lastrack -i flightline1.laz ^
         -i flightline2.laz ^
         -track Velodyne001_above719_012.laz ^
         -store_xyz_range_as_extra_bytes ^
         -odix _xyz_range -olaz ^
         =cores 2

lasgrid -i flightline*_xyz_range.laz -merged ^
        -drop_attribute_below 0 1 ^
        -attribute0 -lowest ^
        -false -set_min_max 20 200 ^
        -o quality/closest_xyz_range_020ft_200ft.png

lasgrid -i flightline*_xyz_range.laz -merged ^
        -drop_attribute_below 0 1 ^
        -attribute0 -highest ^
        -false -set_min_max 30 300 ^
        -o quality/farthest_xyz_range_030ft_300ft.png

Below the complete processing pipeline for creating a median ground model from the “fluffy” Velodyne LiDAR data that results in the hillshaded DTM shown here. The workflow is similar to those we have developed in earlier blog posts for Velodyne Puck based systems like the Hovermap and the Yellowscan.

Hillshaded DTM with a resolution of 1 foot generated via a median ground computation by the LAStools processing pipeline detailed below.

lastile -i flightline1.laz ^
        -i flightline2.laz ^
        -faf ^
        -tile_size 250 -buffer 25 -flag_as_withheld ^
        -odir tiles_raw -o somer.laz

lasnoise -i tiles_raw\*.laz ^
         -step_xy 2 -step 1 -isolated 9 ^
         -odir tiles_denoised -olaz ^
          -cores 4

lasthin -i tiles_denoised\*.laz ^
        -ignore_class 7 ^
        -step 1 -percentile 0.05 10 -classify_as 8 ^
        -odir tiles_thinned_1_foot -olaz ^
        -cores 4

lasthin -i tiles_thinned_1_foot\*.laz ^
        -ignore_class 7 ^
        -step 2 -percentile 0.05 10 -classify_as 8 ^
        -odir tiles_thinned_2_foot -olaz ^
        -cores 4

lasthin -i tiles_thinned_2_foot\*.laz ^
        -ignore_class 7 ^
        -step 4 -percentile 0.05 10 -classify_as 8 ^
        -odir tiles_thinned_4_foot -olaz ^
        -cores 4

lasthin -i tiles_thinned_4_foot\*.laz ^
        -ignore_class 7 ^
        -step 8 -percentile 0.05 10 -classify_as 8 ^
        -odir tiles_thinned_8_foot -olaz ^
        -cores 4

lasground -i tiles_thinned_8_foot\*.laz ^
          -ignore_class 1 7 ^
          -town -extra_fine ^
          -odir tiles_ground_lowest -olaz ^
          -cores 4

lasheight -i tiles_ground_lowest\*.laz ^
          -classify_between -0.05 0.5 6 ^
          -odir tiles_ground_thick -olaz ^
          -cores 4

lasthin -i tiles_ground_thick\*.laz ^
        -ignore_class 1 7 ^
        -step 1 -percentile 0.5 -classify_as 2 ^
        -odir tiles_ground_median -olaz ^
        -cores 4

las2dem -i tiles_ground_median\*.laz ^
        -keep_class 2 ^
        -step 1 -use_tile_bb ^
        -odir tiles_dtm -obil ^
        -cores 4

blast2dem -i tiles_dtm\*.bil -merged ^
          -step 1 -hillshade ^
          -o dtm_hillshaded.png

We thank Harris Aerial for sharing this LiDAR data set with us flown by their Carrier H4 Hybrid HE drone carrying a Snoopy A series LiDAR system from LidarUSA.

RIEGL Becomes LASzip Sponsor for LAS 1.4 Extension

PRESS RELEASE (for immediate release)
August 31, 2015
rapidlasso GmbH, Gilching, Germany

We are happy to announce that RIEGL Laser Measurement Systems, Austria has become a sponsor of the award-winning LASzip compressor. Their contribution at the Silver level will kick-off the actual development phase of the “native LAS 1.4 extension” that had been discussed with the LiDAR community over the past two years. This “native extension” for LAS 1.4 complements the existing “compatibility mode” for LAS 1.4 that was supported by Gold sponsor NOAA and Bronze sponsors Quantum Spatial and Trimble Geospatial. The original sponsor who initiated and financed the open sourcing of the LASzip compressor was USACE – the US Army Corps of Engineers (see http://laszip.org).

The existing “LAS 1.4 compatibility mode” in LASzip was created to provide immediate support for compressing the new LAS 1.4 point types by rewriting them as old point types and storing their new information as “Extra Bytes”. As an added side-benefit this has allowed legacy software without LAS 1.4 support to readily read these newer LAS files as most of the important fields of the new point types 6 to 10 can be mapped to fields of the older point types 1, 3, or 5.

In contrast, the new “native LAS 1.4 extension” of LASzip that is now sponsored in part by RIEGL will utilize the “natural break” in the format due to the new point types of LAS 1.4 to introduce entirely new features such as “selective decompression”, “rewritable classifications and flags”, “integrated spatial indexing”, … and other functionality that has been brain-stormed with the community since rapidlasso GmbH had issued the open “call for input” on native LASzip compression for LAS 1.4 in January 2014. We invite you to follow the progress or contribute to the development via the discussions in the “LAS room“.

silverLASzip_m60_512_275

About rapidlasso GmbH:
Technology powerhouse rapidlasso GmbH specializes in efficient LiDAR processing tools that are widely known for their high productivity. They combine robust algorithms with efficient I/O and clever memory management to achieve high throughput for data sets containing billions of points. The company’s flagship product – the LAStools software suite – has deep market penetration and is heavily used in industry, government agencies, research labs, and educational institutions. Visit http://rapidlasso.com for more information.

About RIEGL:
Austrian based RIEGL Laser Measurement Systems is a performance leader in research, development and production of terrestrial, industrial, mobile, bathymetric, airborne and UAS-based laser scanning systems. RIEGL’s innovative hard- and software provides powerful solutions for nearly all imaginable fields of application. Worldwide sales, training, support and services are delivered from RIEGL‘s Austrian headquarters and its offices in Vienna, Salzburg, and Styria, main offices in the USA, Japan, and in China, and by a worldwide network of representatives covering Europe, North and South America, Asia, Australia and Africa. Visit http://riegl.com for more information.

Discriminating Vegetation from Buildings

I came across an interesting blog article by Jarlath O’Neil-Dunne from the University of Vermont on how LiDAR return information can be used as a simple way to discriminate vegetated areas from buildings. He first computes a normalized first-return DSM and a normalized last-return DSM that he subtracts from another to highlight the vegetation. He writes “This is because the height difference of the first and last returns for buildings is often identical, whereas for trees it is typically much greater.”

Side note: I am not entirely happy with the terminology of a “Normalized Digital Terrain Model (nDTM)”. Jarlath writes: “A similar approach is used to create a Normalized Digital Terrain Model (nDTM).  A DTM is generated from the last returns. The DEM is then subtracted from the DTM to create the nDTM.” I like to reserve the term “Digital Terrain Model (DTM)” for bare-earth terrain computed from returns classified as ground.

Below I radically simplify Jarlath workflow by eliminating the two normalization steps. This not only saves the creation of 3 temporary rasters but also removes the requirement to have ground-classified LiDAR:

  1. Create a first-return frDSM
  2. Create a last-return lrDSM
  3. Subtract the lrDSM from the frDSM to get a return-difference rdDEM

This rdDEM has non-zero heights in all areas where the LiDAR produced more than one return. This happens most often and most pronounced in vegetated areas. Here is how to implement this with las2dem of LAStools:

las2dem -i ..\data\fusa.laz -first_only -o frDSM.bil
las2dem -i ..\data\fusa.laz -last_only -o lrDSM.bil
lasdiff -i frDSM.bil -i lrDSM.bil -o rdDEM.laz
lasview -i rdDEM.laz
rdDEM

The return-difference rdDEM shows the height difference between first and last returns.

Does this work well for you? The results on the “fusa.laz” data set are not entirely convincing … maybe because the vegetation was too dense (leaf-on?) so that the LiDAR penetration is not as pronounced. You can switch back and forth between the first-return and the last-return DSM by loading both *.bil files into lasview with the ‘-files_are_flightlines’ option and then press hotkeys ‘0’ and ‘1’ to toggle between the points and ‘t’ to triangulate the selected DSM.

lasview -i frDSM.bil lrDSM.bil -files_are_flightlines
first-return DSM

first-return DSM

last-return DSM

last-return DSM

We should point out that for Jarlath the return difference raster rdDEM is just one part of the pipeline that is followed by an object-based approach in which they integrate the spectral information from aerial imagery and then use iterative expert systems to further improve the tree canopy classification.

Nevertheless, we believe that our way of classifying vegetation and buildings via a pipeline of lasground, lasheight, and lasclassify gives a better and more robust initial guess than multi-return height differences towards what is vegetation and what are buildings. Below you see this is implemented using the new LASlayers concept:

lasground -i ..\data\fusa.laz -city -extra_fine -olay
lasheight  -i ..\data\fusa.laz -ilay -olay
lasclassify -i ..\data\fusa.laz -ilay -olay 
lasview -i ..\data\fusa.laz -ilay
Automated building and vegetation classification with lasclassify.

Automated building and vegetation classification with lasclassify.

Using lasgrid there are many ways that can easily turn the classified point cloud into a raster so that it can be used for subsequent exploitation together with other image data using a raster processing software. An example is shown below.

lasgrid -i ..\data\fusa.laz -ilay -keep_class 5 ^
        -step 0.5 -subcircle 0.1 -occupancy -fill 1 -false ^
        -use_bb -o vegetation.tif
lasgrid -i ..\data\fusa.laz -ilay -keep_class 6 ^
        -step 0.5 -subcircle 0.1 -occupancy -fill 1 -gray ^
        -use_bb -o buildings.tif
gdalwarp vegetation.tif buildings.tif classified.tif

classified

Alternatively we can use lasboundary to create a shapefile describing either the vegetation or the buildings.

lasboundary -i ..\data\fusa.laz -ilay -keep_class 5 ^
            -disjoint -concavity 1.5 -o vegetation.shp
lasboundary -i ..\data\fusa.laz -ilay -keep_class 6 ^
            -disjoint -concavity 1.5 -o buildings.shp
SHP file generated with lasboundary with polygons describing the vegetation.

SHP file generated with lasboundary with polygons describing the vegetation.

SHP file generated with lasboundary with polygons describing the buildings.

SHP file generated with lasboundary with polygons describing the buildings.

Tutorial: editing LAS or LAZ files “by hand” with lasview

This tutorial describes how to manually edit LiDAR using the new inspection and editing functionality available in ‘lasview.exe’ with the latest release of LAStools (version 140301). We will work with the familiar ‘fusa.laz’ sample LiDAR data set from the LAStools distribution that was recently reported to have shown strange symptoms assumed to be side-effects of the LAZ cloning experiments in the ESRI labs … (-;

Inspecting LiDAR files with cross sections

Copy ‘fusa.laz’ from the folder ‘lastools\data’ to the folder ‘lastools\bin’. Run ‘lasview.exe’ so that it loads ‘fusa.laz’. Either do this via the GUI by double-clicking ‘lasview.exe’, loading ‘fusa.laz’ via the ‘browse …’ menu, and then clicking the ‘VIEW’ button or by entering the command below:

C:\lastools\bin>lasview -i fusa.laz

Press the <x> key to toggle to “select cross” where you can pick a rectangular “cross” section. The default cross section is a profile extending across the bottom of the bounding box.

tutorial4_lasview_01_select_cross

By pressing the <x> key again you toggle back to actually view the cross section. Holding down <ALT> you can rotate the view to look at the cross section from the side. Holding down <CTRL> you can zoom in and out. Holding down <SHIFT> you can translate up and down or left and right. Increase or decrease the size of the points pressing <=> or <->. Hover with the mouse over a point and press <i> to inspect its coordinates and attributes.

tutorial4_lasview_02_cross_inspect_point

Traverse the LiDAR file visually by moving the cross section with the arrow keys <UP> <DOWN> <LEFT> and <RIGHT>. You can move either in the “select cross” view and see the picked rectangle move or in the “cross” view and “walk” through the LiDAR. Hold down the <SHIFT> key simultaneously to take bigger steps or the <ALT> key to take smaller steps. Inspect other points by hovering over them with the cursor and pressing <i>. The point information disappears when pressing <i> with the cursor over the background.

tutorial4_lasview_03_traverse_with_arrow_keysToggle back to the “select cross” view with <x> and pick approximately the same rectangle as shown below:

tutorial4_lasview_04_pick_missclassified_roof

Changing Classifications and Deleting Points

Continuing the steps above, toggle back to the “cross” view by pressing <x>. Note that part of the roof of the house has been miss-classified as vegetation while others are left unclassified. Press <e> to turn on the “EDIT” mode and right-click to select “reclassify points as building (6)” via the pop-up menu.

tutorial4_lasview_05_reclassify_menu

Now use the cross-hair cursor to draw a polygonal fence around all points that should be reclassified. Press <ESC> to remove the last vertex of the polygon if you miss-placed it by a mistake.

tutorial4_lasview_06_reclassify_fenceOnce you are happy with your polygon press <r> to register the edit. A note appears informing you how many points had their classification changed. In the top right corner an “undo” counter appears informing you how many changes you can undo by pressing <CTRL-u>. Try it. Immediataly the changes disappear and a “redo” counter appears instead. Press <CTRL-o> to redo the change you have just undone.

tutorial4_lasview_07_edit_result

Press <CTRL-s> to save this edit as a tiny LAY file using the recently introduced LASlayers concept. In case there was already an existing LAY file (that was not applied with ‘-ilay’ when starting lasview) you will be warned and have to press <CTRL-f> to force overwriting it as shown below.

tutorial4_lasview_07a_first_save

Press the key sequence <SHIFT-b>, <t>, and <a> to get the same visuals above.

tutorial4_lasview_07a_second_save

Press <SHIFT-t> to remove the triangulation again. After saving an edit it can no longer be undone via <CTRL-u>. Instead you will have to strip off this particular layer with the layer management available through “laslayers.exe” as described here. Now press <x> to toggle to the “cross select” view.

tutorial4_lasview_08_move_select_cross

Use the <DOWN> arrow to move the selected cross section to the area shown above that has a few unclassified points in the middle of the roof. Press <x> to go back to the “cross” view and try to understand why these points are not part of the roof. Looks like they are from the top of a chimney, and antenna, or a satellite dish as they do not fit the otherwise planar roof.

Assume we need to remove them for some reason. Pan, translate, and zoom the view such that these points can be easily surrounded by a polygon. Now press <d> to enter the “DELETE” mode, fence in these points, and press <r> to register the deletion.

tutorial4_lasview_09_delete_points

It can be tricky to place a clearly seperating polygon and you may be worried about deleting a few orange building points as well. Press <u> to only display the unclassified points before pressing <r> to register the deletion.

tutorial4_lasview_10_delete_points_unclassified_only

Press <a> to see all points again, then delete the other two points by finding a good view point, pressing <d>, and drawing a polygon.

tutorial4_lasview_11_delete_other_points

After registering this deletion of two points your “undo” counter should be at two. Press <CTRL-u> twice to undo this and the last deletion, then press <CTRL-o> twice to redo them both.

tutorial4_lasview_12_delete_other_points_undo_count_2

Now press <CTRL-s> to save this deletion as another layer. It will be appended to the LAY file that already contains one layer with the roof re-classification edit we did first. Press <t> to triangulate the points in the “cross” view. See how nicly flat the triangulated roofs are now that we deleted these 6 chimney points.

tutorial4_lasview_13_second_save

Look at the size of the tiny LAY file called ‘fusa.lay’ that is in the same folder as the ‘fusa.laz’ file. It contains all the edits we have done so far and mine is only 681 bytes in size. The original LAZ file has not changed. Maybe this is all you want for now. You could send only this tiny LAY file to a colleague elsewhere and he or she could apply those changes locally when needed using the ‘-ilay’ switches. For more on this see the LASlayers page.

C:\lastools\bin>lasview -i fusa.laz
C:\lastools\bin>lasview -i fusa.laz -ilay 1
C:\lastools\bin>lasview -i fusa.laz -ilay 2

However, you may want to eventually apply the changes and produce a new LAZ file. This will be a lot slower as it requires rewriting the entire file. It will also make changes permanent. Press <CTRL-a> and a new file is produced called ‘fusa_1.laz’ that has 6 points less than ‘fusa.laz’ and 69 points with a different classification as “building”. One more thing, press <CTRL-x> if you want to toggle between the “cross” section view and the default view.

tutorial4_lasview_14_applied_laslayers You need to have a license to LAStools to save edits for file that contain 1 million points or more.

Warming up for ILMF 2014, rapidlasso puts lean yet plush “LASlayers” on LiDAR

PRESS RELEASE
February 14, 2014
rapidlasso GmbH, Gilching, Germany

As a sweet foretaste to ILMF 2014, the creators of LAStoolsLASzip, and PulseWaves are announcing “LASlayers” already on Valentine’s Day. The new functionality nicely complements their popular and widely-used LiDAR compressor making the compressed LAZ files editable for most practical purposes. LASlayers significantly reduce I/O load for writing modification to LAS or LAZ files, especially when batch-processing LiDAR tiles on many cores in parallel or when sending changes to LiDAR files across bandwidth-limited networks.

Conceptually LASlayers add

LASlayers store modifications and additional  attributes to raw LAS or LAZ files in small LAY files avoiding to replicate data that has not changed.

Most point attributes (e.g. coordinates, intensities, scan angles, GPS times, etc.) are not modified when processing LiDAR. LASlayers avoids re-writing the unchanged portions of LAS or LAZ by storing only the actual differences layer by layer to a new “LAY” file. Changing the point classifications or deleting a few points, for example, can be done with LAY files that are just a tiny fraction of the size of a new LAS or LAZ file. Adding new attributes such as RGB colors or the height-above-ground with LASlayers means only this new information needs to be written.

This also provides simultaneous access to different versions of the data: a LiDAR server or a Web portal can store only a single copy of the raw LiDAR and apply LASlayers as needed on-the-fly, for example, to replace ellipsoidal with orthometric elevations or to add RGB colors.

Even users of other LiDAR processing software can readily take advantage of LASlayers with the new “laslayers.exe” tool that computes the difference between a raw and a modified LAS or LAZ file and expresses it as a LAY file (assuming the point order has not changed). A typical use case is the exchange of modifications to LiDAR files between remote locations such as a vendor in Australia or Canada and a data processing center in China or India. Instead of up- and downloading the entire modified LAS or LAZ files, only the much smaller LAY files need to be send back and forth.

A fully featured prototype of LASlayers is available (10 MB including data) together with three simple exercises that illustrate the concept and allow interested parties to test it already today on their own data.