Category Archives: hardware

Preparing Drone LiDAR from Snoopy by LidarUSA carrying a Velodyne HDL-32E

Posted on by

In March 2019 I was welcoming Nelson Mattie from LiDAR Latinoamerica to Samara who brought along his versatile Snoopy A-Series scanning system by LidarUSA that is based on the Velodyne HDL-32E scanner. We mounted it to his truck for a mobile scan of the core downtown block, Nelson carried it on his shoulder through “Samara Jungle” for a pedestrian scan and we strapped it onto a DJI Matrice 600 to scan this cute beach and surf town from above.

  • preparing hardware and software
  • scanning Samara Jungle on foot
  • scanning my compost pile
  • getting ready for take-off in a side-street of Samara

The scanning part was easy. Getting sensible data out of the ScanLookPC software proved to be quite an Odyssee as I neither had access to nor training with the ScanLookPC software. Hard surfaces such as the roof of the “New China” supermarket looked this wobbly when looking at the points from individual beams of this 32 beam scanner and turned into a complete fuzz when using all points from all beams.

  • returns from beam with ID 15
  • returns from beam with ID 31
  • returns from all beams

I could only scrutinize the LAS files and I found several unrelated errors, such as duplicate returns and non-unique GPS time stamps, but I was unable to fix the wobbles. Frustrated with the vendor support I was ready to give up when out-of-nowhere I suddenly got an email from Luis Hernandez Perez who also worked with these system. He sent me a link to properly exported flight strips and suggested “it was a problem of poor GNSS signal and the solution for me was use the base station IND1 (of IGS) that was 3.2 kilometers away”. After months of struggles I finally had LiDAR data from downtown Samara and look how crisp the roof of the supermarket suddenly was.

As always my standard quality checks are running lasview, lasinfo, lasgrid but most importantly lasoverlap which would reveal the typical flight line misalignments that we can find in most airborne surveys.

lasoverlap ^
-i Samara\Drone\00_raw\*.laz -faf ^
-step 1 -min_diff 0.1 -max_diff 0.2 -no_over ^
-o Samara\Drone\01_quality\overlap_10_20_before.png

As usually, I contacted Andre Jalobeanu from Bayesmap and asked for help with the alignment. A few days later he returned a new and improved set of flight lines to me that he had run through his stripalign software. So I performed the same quality check with lasoverlap once again.

lasoverlap ^
-i Samara\Drone\00_raw_aligned\*.laz -faf ^
-step 1 -min_diff 0.1 -max_diff 0.2 -no_over ^
-o Samara\Drone\01_quality\overlap_10_20_after.png

  • before stripalign
  • after stripalign

Vertical differences of more than 20 centimeters are mapped to saturated blue and red and the improvement in alignment though stripalign is impressive. Note that the road is not white because it was perfectly aligned but because there was no data. A few days before the scan, Samara got a new tar road from the municipality. As we were flying at 70 meters above ground – a little too high for this LiDAR scanner – we did not capture surfaces with low reflectivity and the fresh black tar did not reflect enough photons.

The Velodyne HDL-32E scanner rotates shooting laser beams from 32 different heads and the information which return comes from which head was stored into the “user data” field and into the “point_source ID” field of each point by the exporting software. The LAS format does not have a dedicated field for this information as the supported maximum number of different laser beams is 4 in the new point types 6 through 10 of the latest LAS 1.4 specification (where this field is called the “scanner channel”). But when we use lastile to turn the flight lines into square tiles we will override the “point_source ID” with the flight line number. Also the “user data” field is a fragile place to store important information as lasheight, for example, will store temporary data there. The “extra bytes” concept of the LAS format is perfect to store such an additional attribute and the ASPRS LAS Working Group is currently discussing to have standardized “extra bytes” for exactly such laser beam IDs.

We at rapidlasso have already implemented this a while ago into our LAStools software. So before processing the data any further we copy the beam index that ranges from 0 to 31 from the “user data” field into an unsigned 8 bit integer “extra byte” with these two las2las commands.

las2las ^
-i Samara\Drone\00_raw_aligned\*.laz ^
-add_attribute 1 "laser beam ID" "which beam ranged this return" ^
-odir Samara\Drone\00_raw_temp -olaz

las2las ^
-i Samara\Drone\00_raw_temp\*.laz ^
-copy_user_data_into_attribute 0 ^
-set_user_data 0 ^
-set_point_source 0 ^
-odir Samara\Drone\00_raw_ready -olaz

In a future article we will process these aligned and prepared flight lines into a number of products. We thank Nelson Mattie from LiDAR Latinoamerica, Luis Hernandez Perez and Andre Jalobeanu from Bayesmap to help me acquire and fix this data. Several others helped with experiments using their own software and data or contributed otherwise to the discussions in the LAStools user forum. Thanks, guys. This data will soon be available as open data but a sample lasinfo report is already below.

lasinfo (210128) report for 'Samara\Drone\00_raw_ready\flightline_01.laz'
reporting all LAS header entries:
file signature: 'LASF'
file source ID: 1
global_encoding: 0
project ID GUID data 1-4: 00000000-0000-0000-0000-000000000000
version major.minor: 1.2
system identifier: 'LAStools (c) by rapidlasso GmbH'
generating software: 'las2las (version 210315)'
file creation day/year: 224/2019
header size: 227
offset to point data: 527
number var. length records: 2
point data format: 1
point data record length: 29
number of point records: 6576555
number of points by return: 6576555 0 0 0 0
scale factor x y z: 0.001 0.001 0.001
offset x y z: 661826 1092664 14
min x y z: 660986.622 1092595.013 11.816
max x y z: 661858.813 1092770.575 95.403
variable length header record 1 of 2:
reserved 0
user ID 'ScanLook'
record ID 25
length after header 0
description 'ScanLook Point Cloud'
variable length header record 2 of 2:
reserved 0
user ID 'LASF_Spec'
record ID 4
length after header 192
description 'by LAStools of rapidlasso GmbH'
Extra Byte Descriptions
data type: 1 (unsigned char), name "laser beam ID", description: "which beam ranged this return", scale: 1 (not set), offset: 0 (not set)
LASzip compression (version 3.4r3 c2 50000): POINT10 2 GPSTIME11 2 BYTE 2
reporting minimum and maximum for all LAS point record entries …
X -839378 32813
Y -68987 106575
Z -2184 81403
intensity 4 255
return_number 1 1
number_of_returns 1 1
edge_of_flight_line 0 0
scan_direction_flag 0 0
classification 1 1
scan_angle_rank -12 90
user_data 0 0
point_source_ID 0 0
gps_time 537764.192416 537854.547507
attribute0 0 31 ('laser beam ID')
number of first returns: 6576555
number of intermediate returns: 0
number of last returns: 6576555
number of single returns: 6576555
overview over number of returns of given pulse: 6576555 0 0 0 0 0 0
histogram of classification of points:
6576555 unclassified (1)

No Sugarcoating: Sweet LiDAR from RiCOPTER carrying VUX-1UAV over Sugarcane

Posted on by

Recently we saw an interesting LiDAR data set talked about on social media by Chad Netto from Chustz Surveying in New Roads, Louisiana and asked for a copy. It is a LiDAR scan of a sugarcane plantation in Assumption Parish, Louisiana carried out with the VUX-1UAV by RIEGL mounted onto a RiCOPTER and guided by an Applanix IMU and a Trimble base station. That is probably one of the sweetest (but also one of the most expensive) UAV LiDAR system you can buy today. I received this LiDAR file and this trajectory file. In the following we talk a detailed look at this data set.

First we run lasinfo to get an idea of the contents of the data set. We create various histograms as those can often help understand an unfamiliar data set:

lasinfo -i sugarcane\181026_163424.laz ^
        -cd ^
        -histo gps_time 5 ^
        -histo intensity 64 ^
        -histo point_source 1 ^
        -histo z 5 ^
        -odix _info -otxt

You can download the resulting report here. For the 84,751,955 points we notice that

  1. both horizontal and vertical coordinates are stored in US survey feet
  2. with scale factors of 0.00025 this means a resolution of 76 micrometer
  3. there is no explicit flight line information (all point source IDs are zero)
  4. gaps in the GPS time stamp histogram are suggesting multiple lines

First we use las2las to lower the insanely high resolution from 0.00025 US survey feet to something more reasonable for an airborne UAV scan, namely to 0.01 or 1 hundredths of a US survey foot or centi-US-survey-feet:

las2las -i sugarcane\181026_163424.laz ^
        -rescale 0.01 0.01 0.01 ^
        -odix _cft -olaz

I have already done this for you. The file that is online is already in “centi-US-survey-feet” because it reduced the file size from the original 678 MB file that we got from Netto to the 518 MB file that is online, meaning that you had 160 MB less data to download.

Next we use lassplit to recover the original flight lines as follows:

lassplit -i sugarcane\181026_163424.laz ^
         -recover_flightlines ^
         -odir sugarcane\0_recovered_strips ^
         -o assumption.laz

This results in 5 strips. We then use lassort to bring the strips back into their original acquisition order by sorting first based on the GPS time stamp (which brings all returns of one pulse together) and second on the return number (which sorts them in ascending order). We do this on 3 cores in parallel with this command:

lassort -i sugarcane\0_recovered_strips\*.laz ^
        -gps_time ^
        -return_number ^
        -odir sugarcane\1_sorted_strips -olaz ^
        -cores 3

We also create a spatial index for each of these strips using lasindex so that any area-of-interest query that we do later will be significantly accelerated. See the README file for the meaning of the parameters:

lasindex -i sugarcane\1_sorted_strips\*.laz ^
         -tile_size 10 -maximum -100 ^
         -cores 3

Then we check for flight line alignment using lasoverlap by comparing – per 2 feet by 2 feet area – the lowest elevation value of points from different strips wherever there is overlap. Cells with an absolute vertical difference of less than a quarter of a foot are mapped to white. Cells with vertical differences of more (or less) than a quarter foot are mapped to an increasingly red (or blue) color that is saturated red (or blue) when one full foot is reached.

lasoverlap -i sugarcane\1_sorted_strips\*.laz ^
           -files_are_flightlines ^
           -step 2.0 ^
           -min_diff 0.25 -max_diff 1.0 ^
           -o sugarcane\2_quality\overlap.png

The resulting image looks dramatic at first glance. But we have to remember that this is sugarcane. The penetration of the laser can vary greatly depending on the direction from which the beam hits the densely standing stalks. Large differences between flight lines can be expected where sugarcane stands tall. We need to focus our visual quality checks on the few open areas, namely the access roads and harvested areas.

Color-mapped highest vertical difference in lowest point per 2 feet by 2 feet area between overlapping flight lines.

We use las2las via its native GUI to cut out several suspicious-looking open areas with overly red or overly blue shading. By loading the resulting image into the GUI these areas-of-interest are easy to target and cut out.

las2las -i sugarcane\1_sorted_strips\*.laz -gui

Overlaying the difference image in the GUI of las2las to cut out suspicious areas for closer inspection.

We cut out four square 100 by 100 meter tiles in open areas that show a suspiciously strong pattern of red or blue colors for closer inspection. The command lines for these four square areas are given below and you can download them here:

  1. assumption_3364350_534950_100.laz
  2. assumption_3365600_535750_100.laz
  3. assumption_3364900_535500_100.laz
  4. assumption_3365500_535600_100.laz
las2las -i sugarcane\1_sorted_strips\*.laz ^
        -merged -faf ^
        -inside_tile 3364350 534950 100 ^
        -o sugarcane\assumption_3364350_534950_100.laz

las2las -i sugarcane\1_sorted_strips\*.laz ^
        -merged -faf ^
        -inside_tile 3365600 535750 100 ^
        -o sugarcane\assumption_3365600_535750_100.laz

las2las -i sugarcane\1_sorted_strips\*.laz ^
        -merged -faf ^
        -inside_tile 3364900 535500 100 ^
        -o sugarcane\assumption_3364900_535500_100.laz

las2las -i sugarcane\1_sorted_strips\*.laz ^
        -merged -faf ^
        -inside_tile 3365500 535600 100 ^
        -o sugarcane\assumption_3365500_535600_100.laz

In the image sequence below we scrutinize these differences which will lead us to notice two things:

  1. There are vertical miss-alignments of around one foot. These big difference can especially be observed between flight lines that cover an area with a very high point density and those that cover the very same area with a very low point density.
  2. There are horizontal miss-alignments of around one foot. Again these differences seem somewhat correlated to the density that these flight lines cover a particular area with.

Picking a square area with suspicious discoloration to clip into a separate file.
View down the access road. Points colored by elevation. Notice the lighter blue points.
One flight line only.
Another flight line. This one contributes those lighter blue points.
Picking a profile line in lasview (press ‘x’ and then select a line while holding down the ‘SHIFT’ key). Then press ‘x’ again.
The points of the sparser flight line are one foot above the denser one.
Picking another square area with suspicious discoloration. This one in the middle of an empty field.
The plow lines are clearly visible.
Hillshading only the points of one flight line.
Hillshading only the points of another flight line. Notice the horizontal shift.
Picking a profile line in lasview (press ‘x’ and then select a line while holding down the ‘SHIFT’ key). Then press ‘x’ again.
Profile line across the plow lanes of an empty sugarcane plantation.
Profile line for one flight line. The horizontal shift is about two thirds of a foot.
Profile line for another flight line. The horizontal shift is about two thirds of a foot.
Picking another square area with suspicious discoloration. This one includes an access road again.
Already the issue is visible in the miss-classification of the data. There is odd grass growing in the middle of the road.
Also the elevation coloring shows this pattern. Press ‘c’ to step through different coloring options.
The flight line coloring confirms that it is a miss alignment issue again. Press ‘c’ to step through different coloring options.
Picking a profile line in lasview (press ‘x’ and then select a line while holding down the ‘SHIFT’ key). Then press ‘x’ again.
The miss-alignment between the sparse points of one flight line and the dense points of the other is obvious.
The vertical difference is close to a foot.
Picking another suspicious area. We are mostly interested in the discoloration on the empty field..
Picking a profile line in lasview (press ‘x’ and then select a line while holding down the ‘SHIFT’ key). Then press ‘x’ again.
Profile line across the plow lanes of an empty sugarcane plantation.
Profile line with exaggerated z scale across the plow lanes of an empty sugarcane plantation. The horizontal shift is about one foot.
Picking a profile rectangle in lasview (press ‘x’ and then select a rectangle). Then press ‘x’ again.
Hillshading only the points of one flight line.
Hillshading only the points of another flight line. Notice the horizontal shift.
Using all flight lines together would give messy and incorrect results in this area.

For the most part the miss-aligned points come from a flight line that has only sparse coverage in that area. In a flat terrain the return density per area goes down the farther we are from the drone as those areas are only reached with higher and higher scan angles. Hence an immediate idea that comes to mind is to limit the scan angle to those closer to nadir and lower the range from -81 to 84 degrees reported in the lasinfo report to something like -75 to 75, -70 to 70, or -65 to 65 degrees. We can check how this will improve the alignment with these lasoverlap command lines:

lasoverlap -i sugarcane\1_sorted_strips\*.laz ^
           -files_are_flightlines ^
           -keep_scan_angle -75 75 ^
           -step 2.0 ^
           -min_diff 0.25 -max_diff 1.0 ^
           -o sugarcane\2_quality\overlap75.png

lasoverlap -i sugarcane\1_sorted_strips\*.laz ^
           -files_are_flightlines ^
           -keep_scan_angle -70 70 ^
           -step 2.0 ^
           -min_diff 0.25 -max_diff 1.0 ^
           -o sugarcane\2_quality\overlap70.png

lasoverlap -i sugarcane\1_sorted_strips\*.laz ^
           -files_are_flightlines ^
           -keep_scan_angle -65 65 ^
           -step 2.0 ^
           -min_diff 0.25 -max_diff 1.0 ^
           -o sugarcane\2_quality\overlap65.png

lasoverlap -i sugarcane\1_sorted_strips\*.laz ^
           -files_are_flightlines ^
           -keep_scan_angle -60 60 ^
           -step 2.0 ^
           -min_diff 0.25 -max_diff 1.0 ^
           -o sugarcane\2_quality\overlap60.png

This simple technique significantly improves the difference image. Based on these images would suggest to only use returns with a scan angle between -70 and 70 degrees for any subsequent analysis. This seems to remove most of the discoloring in open areas without loosing too many points. Note that only using returns with a scan angle between -60 and 60 degrees means that some flight lines are no longer overlapping each other.

Note also that by limiting the scan angle we get suddenly get white areas even in incredible dense vegetation. The more horizontal a laser shoot is the more likely it will only hit higher up sugarcane plants and the less likely it will penetrate all the way to the ground. The white areas coincide with where laser pulses are close to nadir which is in the flight line overlap areas that directly below the drone’s trajectory.

Can we improve alignment further? Not with LAStools, so I turned to Andre Jalobeanu, a specialist on that particular issue, who I have known for many years. Andre has developed BayesStripAlign – a software by his company BayesMap that is quite complementary to LAStools and does exactly what the name suggests: it align strips. I gave Andre the five flight lines and he aligned them for me. Below the new difference images:

Before strip alignment (all returns).
After strip alignment (all returns).
Before strip alignment (returns with scan angle between -70 and 70 degrees)
After strip alignment (returns with scan angle between -70 and 70 degrees)

We cut out the very same four square areas from the realigned strips for closer inspection and you may investigate them on your own. You can download them here.

  1. assumption_3364350_534950_100_realigned.laz
  2. assumption_3365600_535750_100_realigned.laz
  3. assumption_3364900_535500_100_realigned.laz
  4. assumption_3365500_535600_100_realigned.laz

In the image sequence below we are just looking at the last of the four square areas and you can see that most of the miss-alignment we saw earlier between the flight lines was removed.

Profile line across the plow lanes of an empty sugarcane plantation.
Profile line with exaggerated z scale across the plow lanes of an empty sugarcane plantation.
Hillshaded plow lanes of the returns of one flight line.
Hillshaded plow lanes of the returns of another flight line.

We would like to thank Chad Netto from Chustz Surveying to make this data set available to us and Andre Jalobeanu from BayesMap to align the flight lines for us.

First Look with LAStools at LiDAR from Hovermap Drone by CSIRO

Posted on by

Last December we had a chance to visit the team of Dr. Stefan Hrabar at CSIRO in Pullenvale near Brisbane who work on a drone LiDAR system called Hovermap. This SLAM-based system is mainly developed for the purpose of autonomous flight and exploration of GPS-denied environments such as buildings, mines and tunnels. But as the SLAM algorithm continuously self-registers the scan lines it produces a LiDAR point cloud that in itself is a nice product. We started our visit with a short test flight around the on-site tower. You can download the LiDAR data and the drone trajectory of this little survey here:

The Hovermap system is based on the Velodyne Puck Lite (VLP-16) that is much cheaper and more light-weight than many other LiDAR systems. One interesting tidbit in the Hovermap setup is that the scanner is installed such that the entire Puck is constantly rotating as you can see in this video. But  the Velodyne Puck is also known to produce somewhat “fluffy” surfaces with a thickness of a few centimeters. In a previous blog post with data from the YellowScan Surveyor system (that is also based on the Puck) we used a “median ground” surface to deal with the “fluff”. In the following we will have a look at the LiDAR data produced by Hovermap and how to process it with LAStools.

LiDAR data of CSIRO tower acquired during test flight of Hovermap system.

As always we start with a lasinfo report that computes the average density ‘-cd’ and histograms for the intensity and the GPS time:

lasinfo -i CSIRO_Tower\results.laz ^
        -cd ^
        -histo intensity 16 -histo gps_time 2 ^
        -odir CSIRO_Tower\quality -odix _info -otxt

A few excerpts of the resulting lasinfo report that you can download here are below:

lasinfo (180409) report for 'CSIRO_Tower\results.laz'
[...]
 number of point records: 16668904
 number of points by return: 0 0 0 0 0
 scale factor x y z: 0.0001 0.0001 0.0001
 offset x y z: -5.919576153930379 22.785394470724583 9.535698734939086
 min x y z: -138.6437 -125.2552 -34.1510
 max x y z: 126.8046 170.8260 53.2224
WARNING: full resolution of min_x not compatible with x_offset and x_scale_factor: -138.64370561381907
WARNING: full resolution of min_y not compatible with y_offset and y_scale_factor: -125.25518631070418
WARNING: full resolution of min_z not compatible with z_offset and z_scale_factor: -34.150966206894068
WARNING: full resolution of max_x not compatible with x_offset and x_scale_factor: 126.80455330595831
WARNING: full resolution of max_y not compatible with y_offset and y_scale_factor: 170.82597525215334
WARNING: full resolution of max_z not compatible with z_offset and z_scale_factor: -34.150966206894068
[...]
 gps_time 121.288045 302.983110
WARNING: 2 points outside of header bounding box
[...]
covered area in square units/kilounits: 51576/0.05
point density: all returns 323.19 last only 318.40 (per square units)
 spacing: all returns 0.06 last only 0.06 (in units)
WARNING: for return 1 real number of points by return is 16424496 but header entry was not set.
WARNING: for return 2 real number of points by return is 244408 but header entry was not set.
[...]
real max z larger than header max z by 0.000035
real min z smaller than header min z by 0.000035
[...]

Most of these warnings have to do with poorly chosen offset values in the LAS header that have many decimal digits instead of being nice round numbers. The points are stored with sub-millimeter resolution (scale factors of 0.0001) which is unnecessarily precise for a UAV flying a Velodyne Puck where the overall system error can be expected to be on the order of a few centimeters. Also the histogram of return numbers in the LAS header was not populated. We can fix these issues with one call to las2las:

las2las -i CSIRO_Tower\results.laz ^
        -rescale 0.01 0.01 0.01 ^
        -auto_reoffset ^
        -odix _fixed -olaz

If you create another lasinfo report on the fixed file you will see that all the warnings have gone. The file size is now also only 102 MB instead of 142 MB because centimeter coordinate compress much better than sub-millimeter coordinates.

The average density of 318 last return per square meter reported by lasinfo is not that useful for a UAV survey because it does account for the highly varying distribution of LiDAR returns in the area surveyed. With lasgrid we can get a much more clear picture of that.

lasgrid -i CSIRO_Tower\results_fixed.laz ^
        -last_only ^
        -step 0.5 -use_bb -density ^
        -false -set_min_max 0 1500 ^
        -o CSIRO_Tower\quality\density_0_1500.png

LiDAR density: blue is close to zero and red is 1500 or more last returns / sqr mtr

The red dot in the point density indicated an area with over 1500 last returns per square meter. No surprise that this is the take-off and touch-down location of the copter drone. Naturally this spot is completely over-scanned compared to the rest of the area. We can remove these points with the help of the timestamps by cutting off the start and the end of the recording.

The total recording time including take-off, flight around the tower, and touch-down was around 180 seconds or 3 minutes as the lasinfo report tells us. Note that the recorded time stamps are neither “GPS Week Time” nor “Adjusted Standard GPS Time” but an internal system time. By visualizing the trajectory of the UAV with lasview while binning the timestamps into the intensity field we can easily determine what interval of timestamps describes the actual survey flight. First we convert the drone trajectory from the textual ASCII format to the LAZ format with txt2las:

txt2las -i CSIRO_Tower\results_traj.txt ^
        -skip 1 ^
        -parse txyz ^
        -set_classification 12 ^
        -olaz

lasview -i CSIRO_Tower\results_traj.laz ^
        -bin_gps_time_into_intensity 1

Binning timestamps into intensity allows visually determining start and end of survey.

Using lasview and pressing <i> while hovering over those points of the trajectory that appear to be the survey start and end we determine visually that the timestamps between 164 to 264 correspond to the actual survey flight over the area of interest with the take-off and touch-down maneuvers excluded. We use las2las to cut out the relevant part and re-run lasgrid:

las2las -i CSIRO_Tower\results_fixed.laz ^
        -keep_gps_time 164 264 ^
        -o CSIRO_Tower\results_survey.laz

lasgrid -i CSIRO_Tower\results_survey.laz ^
        -last_only ^
        -step 0.5 -use_bb -density ^
        -false -set_min_max 0 1500 ^
        -o CSIRO_Tower\quality\density_0_1500_survey.png

LiDAR density after removing take-off and touch-down maneuvers.

The other set of point we are less interested in are those occasional hits far from the scanner that sample the area too sparsely to be useful for anything. We use lastrack to reclassify points as noise (7) that exceed a x/y distance of 50 meters from the trajectory and then use lasgrid to create another density image without the points classified as noise..

lastrack -i CSIRO_Tower\results_survey.laz ^
         -track CSIRO_Tower\results_traj.laz ^
         -classify_xy_range_between 50 1000 7 ^
         -o CSIRO_Tower\results_xy50.laz

lasgrid -i CSIRO_Tower\results_xy50.laz ^
        -last_only -keep_class 0 ^
        -step 0.5 -use_bb -density ^
        -false -set_min_max 0 1500 ^
        -o CSIRO_Tower\quality\density_0_1500_xy50.png

LiDAR density after removing returns farther than 50 m from trajectory.

We process the remaining points using a typical tile-based processing pipeline. First we run lastile to create tiling of 200 meter by 200 meter tiles with 20 buffers while dropping the noise points::

lastile -i CSIRO_Tower\results_xy50.laz ^
        -drop_class 7 ^
        -tile_size 200 -buffer 20 -flag_as_withheld ^
        -odir CSIRO_Tower\tiles_raw -o eta.laz

Because of the high sampling we expect there to be quite a few duplicate point where all three coordinate x, y, and z are identical. We remove them with a call to lasduplicate:

lasduplicate -i CSIRO_Tower\tiles_raw\*.laz ^
             -unique_xyz ^
             -odir CSIRO_Tower\tiles_unique -olaz ^
             -cores 4

This removes between 12 to 25 thousand point from each tile. Then we use lasnoise to classify isolated points as noise:

lasnoise -i CSIRO_Tower\tiles_unique\*.laz ^
         -step_xy 0.5 -step_z 0.1 -isolated 5 ^
         -odir CSIRO_Tower\tiles_denoised_temp -olaz ^
         -cores 4

Aggressive parameters assure most noise point below ground are found.

This classifies between 13 to 23 thousand point from each tile into the noise classification code 7. We use rather aggressive settings to make sure we get most of the noise points that are below the terrain. Getting a correct ground classification in the next few steps is the main concern now even if this means that many points above the terrain on wires, towers, or vegetation will also get miss-classified as noise (at least temporarily). Next we use lasthin to classify a subset of points with classification code 8 on which we will then run the ground classification. We classify each point that is closest to the 5th percentile in elevation per 25 cm by 25 cm grid cell given there are at least 20 non-noise points in a cell. We then repeat this while increasing the cell size to 50 cm by 50 cm and 100 cm by 100 cm.

lasthin -i CSIRO_Tower\tiles_denoised_temp\*.laz ^
        -ignore_class 7 ^
        -step 0.25 -percentile 5 20 -classify_as 8 ^
        -odir CSIRO_Tower\tiles_thinned_025 -olaz ^
        -cores 4

lasthin -i CSIRO_Tower\tiles_thinned_025\*.laz ^
        -ignore_class 7 ^
        -step 0.50 -percentile 5 20 -classify_as 8 ^
        -odir CSIRO_Tower\tiles_thinned_050 -olaz ^
        -cores 4

lasthin -i CSIRO_Tower\tiles_thinned_025\*.laz ^
        -ignore_class 7 ^
        -step 1.00 -percentile 5 20 -classify_as 8 ^
        -odir CSIRO_Tower\tiles_thinned_100 -olaz ^
        -cores 4

 

All points of tile ‘eta_0_0.laz’.
TIN of points classified in three calls to lasthin.
TIN of ground points classified by lasground.

Then we ground classify the points that were classified into the temporary classification code 8 in the previous step using lasground.

lasground -i CSIRO_Tower\tiles_thinned_100\*.laz ^
          -ignore_class 7 0 ^
          -town -ultra_fine ^
          -odir CSIRO_Tower\tiles_ground -olaz ^
          -cores 4

The resulting ground points are a lower envelope of the “fluffy” sampled surfaces produced by the Velodyne Puck scanner. We use lasheight to thicken the ground by moving all points between 1 cm below and 6 cm above the TIN of these “low ground” points to a temporary classification code 6 representing a “thick ground”. We also undo the overly aggressive noise classifications above the ground by setting all higher points back to classification code 1 (unclassified).

lasheight -i CSIRO_Tower\tiles_ground\*.laz ^
          -classify_between -0.01 0.06 6 ^
          -classify_above 0.06 1 ^
          -odir CSIRO_Tower\tiles_ground_thick -olaz ^
          -cores 4

Profile view for 25 centimeter wide strip of open terrain. Top: Green points are low ground. Orange points are thickened ground with 5 cm drop lines. Middle: Brown points are median ground computed from thick ground. Bottom: Comparing low ground points (in green) with median ground points (in brown).

From the “thick ground” we then compute a “median ground” using lasthin in a similar fashion as we used it before. A profile view for a 25 centimeter wide strip of open terrain illustrates the workflow of going from “low ground” via “thick ground” to “median ground” and shows the slight difference in elevation between the two.

lasthin -i CSIRO_Tower\tiles_ground_thick\*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.25 -percentile 50 10 -classify_as 2 ^
        -odir CSIRO_Tower\tiles_ground_median_025 -olaz ^
        -cores 4

lasthin -i CSIRO_Tower\tiles_ground_median_025\*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.50 -percentile 50 10 -classify_as 2 ^
        -odir CSIRO_Tower\tiles_ground_median_050 -olaz ^
        -cores 4

lasthin -i CSIRO_Tower\tiles_ground_median_050\*.laz ^
        -ignore_class 0 1 7 ^
        -step 1.00 -percentile 50 10 -classify_as 2 ^
        -odir CSIRO_Tower\tiles_ground_median_100 -olaz ^
        -cores 4

Then we use lasnoise once more with more conservative settings to remove the noise points that are sprinkled around the scene.

lasnoise -i CSIRO_Tower\tiles_ground_median_100\*.laz ^
         -step_xy 1.0 -step_z 1.0 -isolated 5 ^
         -odir CSIRO_Tower\tiles_denoised -olaz ^
         -cores 4

While we classify the scene into building roofs, vegetation, and everything else with lasclassify we also move all (unused) classifications to classification code 1 (unclassified). You may play with the parameters of lasclassify (see README) to achieve better a building classification. However, those buildings the laser can peek into (either via a window or because they are gazebo-like structures) will not be classified correctly. unless you remove the points that are under the roof somehow.

lasclassify -i CSIRO_Tower\tiles_denoised\*.laz ^
            -ignore_class 7 ^
            -change_classification_from_to 0 1 ^
            -change_classification_from_to 6 1 ^
            -step 1 ^
            -odir CSIRO_Tower\tiles_classified -olaz ^
            -cores 4

A glimpse at the final classification result is below. A hillshaded DTM and a strip of classified points. Of course the tower was miss-classified as vegetation given that it looks just like a tree to the logic used in lasclassify.

The hillshaded DTM with a strip of classified points.

Finally we remove the tile buffers (that were really important for tile-based processing) with lastile:

lastile -i CSIRO_Tower\tiles_classified\*.laz ^
        -remove_buffer ^
        -odir CSIRO_Tower\tiles_final -olaz ^
        -cores 4

And publish the LiDAR point cloud as version 1.6 of Potree using laspublish:

laspublish -i CSIRO_Tower\tiles_final\*.laz ^
           -i CSIRO_Tower\results_traj.laz ^
           -only_3D -elevation -overwrite -potree16 ^
           -title "CSIRO Tower" ^
           -description "HoverMap test flight, 18 Dec 2017" ^
           -odir CSIRO_Tower\tiles_portal -o portal.html -olaz

Note that we also added the trajectory of the drone because it looks nice and gives a nice illustration of how the UAV was scanning the scene.

Via Potree we can publish and explore the final point cloud using any modern Web browser.

We would like to thank the entire team around Dr. Stefan Hrabar for taking time out of their busy schedules just a few days before Christmas.

Processing Drone LiDAR from YellowScan’s Surveyor, a Velodyne Puck based System

Posted on by

Points clouds from UAVs have become a common sight. Cheap consumer drones equipped with cameras produce points from images with increasing quality as photogrammetry software is improving. But vegetation is always a show stopper for point clouds generated from imagery data. Only an active sensing technique such as laser scanning can penetrate through the vegetation and generate points on the ground under the canopy of a forested area. Advances in UAV technology and the miniaturization of LiDAR systems have allowed lasers-scanning solutions for drones to enter the market.

Last summer we attended the LiDAR for Drone 2017 Conference by YellowScan and processed some data sets flown with their Surveyor system that is built around the Velodyne VLP-16 Puck LiDAR scanner and the Applanix APX15 single board GNSS-Inertial solution. One common challenge observed in LiDAR data generated by the Velodyne Puck is that surfaces are not as “crisp” as those generated by other laser scanners. Flat and open terrain surfaces are described by a layer of points with a “thickness” of a few centimeter as you can see in the images below. This visualization uses a 10 meter by 5 meter cut-out of from this data set with the coordinate range [774280,774290) in x and [6279463,6279468) in y. Standard ground classification routines will “latch onto” the lowermost envelope of these thick point layers and therefore produce a sub-optimal Digital Terrain Model (DTM).

Triangulation of all last returns of mostly flat 10 meter by 5 meter area. Points are colored by flightline ID.
Zooming shows the vertical spread of the points. Points are colored by flightline ID.
Zooming further in and pink like with length of 8 cm for scale. Points are colored by flightline ID.
Triangulation of all last returns of mostly flat 10 meter by 5 meter area
Triangulation of highest returns on 10 cm by 10 cm grid and lowest returns on 10 cm by 10 cm grid with 5 cm “rise lines”
Triangulation of highest returns on 10 cm by 10 cm grid and lowest returns on 10 cm by 10 cm grid with 10 cm “rise lines”
Triangulation of lowest returns on 10 cm by 10 cm grid and highest returns on 10 cm by 10 cm grid with 5 cm “drop lines”
Triangulation of lowest returns on 10 cm by 10 cm grid and highest returns on 10 cm by 10 cm grid with 10 cm “drop lines”

In part this “thickness” can be reduced by using fewer flightlines as the “thickness” of each flightline by itself is lower but it is compounded when merging all flightlines together. However, deciding which (subset of) flightlines to use for which part of the scene to generate the best possible ground model is not an obvious tasks either and even per flightline there will be a remaining “thickness” to deal with as can be seen in the following set of images.

All returns from flightline 1 triangulated.
All returns from flightline 2 triangulated.
All returns from flightline 4 triangulated.
All returns from flightline 5 triangulated.
All returns from flightline 6 triangulated.
All returns from all flightlines triangulated.

In the following we show how to deal with “thickness” in a layer of points describing a ground surface. We first produce a “lowest ground” which we then widen into a “thick ground” from which we then derive “median ground” points that create a plausible terrain representation when interpolated by a Delaunay triangulation and rasterized onto a DTM. Step by step we process this example data set captured in a “live demo” during the LiDAR for Drone 2017 Conference – the beautiful Château de Flaugergues in Montpellier, France where the event took place. You can download this data via this link if you would like to repeat these processing steps:

Once you decompress the RAR file (e.g. with the UnRar.exe freeware) you will find six raw flight strips in LAS format and the trajectory of the UAV in ASCII text format as it was provided by YellowScan.

E:\LAStools\bin>dir Flaugergues
06/27/2017 08:03 PM 146,503,985 Flaugergues_test_demo_ppk_L1.las
06/27/2017 08:02 PM  91,503,103 Flaugergues_test_demo_ppk_L2.las
06/27/2017 08:03 PM 131,917,917 Flaugergues_test_demo_ppk_L3.las
06/27/2017 08:03 PM 219,736,585 Flaugergues_test_demo_ppk_L4.las
06/27/2017 08:02 PM 107,705,667 Flaugergues_test_demo_ppk_L5.las
06/27/2017 08:02 PM  74,373,053 Flaugergues_test_demo_ppk_L6.las
06/27/2017 08:03 PM   7,263,670 Flaugergues_test_demo_ppk_traj.txt

As usually we start with quality checking by visual inspection with lasview and by creating a textual report with lasinfo.

E:\LAStools\bin>lasview Flaugergues_test_demo_ppk_L1.las

The raw LAS file “Flaugergues_test_demo_ppk_L1.las” colored by elevation.

E:\LAStools\bin>lasinfo Flaugergues_test_demo_ppk_L1.las
lasinfo (171011) report for Flaugergues_test_demo_ppk_L1.las
reporting all LAS header entries:
 file signature: 'LASF'
 file source ID: 1
 global_encoding: 1
 project ID GUID data 1-4: 00000000-0000-0000-0000-000000000000
 version major.minor: 1.2
 system identifier: 'YellowScan Surveyor'
 generating software: 'YellowReader by YellowScan'
 file creation day/year: 178/2017
 header size: 227
 offset to point data: 297
 number var. length records: 1
 point data format: 3
 point data record length: 34
 number of point records: 4308932
 number of points by return: 4142444 166488 0 0 0
 scale factor x y z: 0.001 0.001 0.001
 offset x y z: 774282 6279505 92
 min x y z: 774152.637 6279377.623 82.673
 max x y z: 774408.344 6279541.646 116.656
variable length header record 1 of 1:
 reserved 0
 user ID 'LASF_Projection'
 record ID 34735
 length after header 16
 description ''
 GeoKeyDirectoryTag version 1.1.0 number of keys 1
 key 3072 tiff_tag_location 0 count 1 value_offset 2154 - ProjectedCSTypeGeoKey: RGF93 / Lambert-93
reporting minimum and maximum for all LAS point record entries ...
 X -129363 126344
 Y -127377 36646
 Z -9327 24656
 intensity 0 65278
 return_number 1 2
 number_of_returns 1 2
 edge_of_flight_line 0 0
 scan_direction_flag 0 0
 classification 0 0
 scan_angle_rank -120 120
 user_data 75 105
 point_source_ID 1 1
 gps_time 219873.160527 219908.550379
 Color R 0 0
 G 0 0
 B 0 0
number of first returns: 4142444
number of intermediate returns: 0
number of last returns: 4142444
number of single returns: 3975956
overview over number of returns of given pulse: 3975956 332976 0 0 0 0 0
histogram of classification of points:
 4308932 never classified (0)

Nicely visible are the circular scanning patterns of the Velodyne VLP-16 Puck. We also notice that the trajectory of the UAV can be seen in the lasview visualization because the Puck was scanning the drone’s own landing gear. The lasinfo report tells us that point coordinates are stored with too much resolution (mm) and that points do not need to be stored using point type 3 (with RGB colors) because all RGB values are zero. We fix this with an initial run of las2las and also compress the raw strips to the LAZ format on 4 CPUs in parallel.

las2las -i Flaugergues\*.las ^
        -rescale 0.01 0.01 0.01 ^
        -auto_reoffset ^
        -set_point_type 1 ^
        -odir Flaugergues\strips_raw -olaz ^
        -cores 4

Next we do the usual check for flightline alignment with lasoverlap (README) which we consider to be by far the most important quality check. We compare the lowest elevation from different flightline per 25 cm by 25cm cell in all overlap areas. We consider a vertical difference of up to 5 cm as acceptable (color coded as white) and mark differences of over 30 cm (color coded as saturated red or blue).

lasoverlap -i Flaugergues\strips_raw\*.laz -faf ^
           -step 0.25 ^
           -min_diff 0.05 -max_diff 0.3 ^
           -odir Flaugergues\quality -o overlap.png

Color codes number of overlapping flightlines.
Color codes vertical difference between flightlines.

The vertical difference in open areas between the flightlines is slightly above 5 cm which we consider acceptable in this example. Depending on the application we recommend to investigate further where these differences come from and what consequences they may have for post processing. We also create a color-coded visualization of the last return density per 25 cm by 25 cm cell using lasgrid (README) with blue meaning less than 100 returns per square meter and red meaning more than 4000 returns per square meter.

lasgrid -i Flaugergues\strips_raw\*.laz -merged ^
        -keep_last ^
        -step 0.25 ^
        -point_density ^
        -false -set_min_max 100 4000 ^
        -odir Flaugergues\quality -o density_100_4000.png

Color coded density of last returns per square meter for each 25 cm by 25 cm cell. Blue means 100 or less last returns per square meter. Red means 4000 or more last returns per square meter

As usual we start the LiDAR processing by reorganizing the flightlines into square tiles. Because of the variability in the density that is evident in the visualization above we use lastile (README) to create an adaptive tiling that starts with 200 m by 200 m tiles and then iterate to refine those tiles with over 10 million points down to smaller 25 m by 25 m tiles.

lastile -i Flaugergues\strips_raw\*.laz ^
        -apply_file_source_ID ^
        -tile_size 200 -buffer 8 -flag_as_withheld ^
        -refine_tiling 10000000 ^
        -odir Flaugergues\tiles_raw -o flauge.laz

lastile -i Flaugergues\tiles_raw\flauge*_200.laz ^
        -refine_tiles 10000000 ^
        -olaz ^
        -cores 4

lastile -i Flaugergues\tiles_raw\flauge*_100.laz ^
        -refine_tiles 10000000 ^
        -olaz ^
        -cores 4

lastile -i Flaugergues\tiles_raw\flauge*_50.laz ^
        -refine_tiles 10000000 ^
        -olaz ^
        -cores 4

Subsequent processing is faster when the points have a spatially coherent order. Therefore we rearrange the points into standard space-filling z-order using a call to lassort (README). We run this in parallel on as many cores as it makes sense (i.e. not using more cores than there are physical CPUs).

lassort -i Flaugergues\tiles_raw\flauge*.laz ^
        -odir Flaugergues\tiles_sorted -olaz ^
        -cores 4

Next we classify those points as noise that are isolated on a 3D grid of 1 meter cell size using lasnoise. See the README file of lasnoise for a description on the exact manner in which the isolated points are classified. We do this to eliminate low noise points that would otherwise cause trouble in the subsequent processing.

lasnoise -i Flaugergues\tiles_sorted\flauge*.laz ^
         -step 1 -isolated 5 ^
         -odir Flaugergues\tiles_denoised -olaz ^
         -cores 4

Next we mark the subset of lowest points on a 2D grid of 10 cm cell size with classification code 8 using lasthin (README) while ignoring the noise points with classification code 7 that were marked as noise in the previous step.

lasthin -i Flaugergues\tiles_denoised\flauge*.laz ^
        -ignore_class 7 ^
        -step 0.1 -lowest ^
        -classify_as 8 ^
        -odir Flaugergues\tiles_lowest -olaz ^
        -cores 4

Considering only the resulting points marked with classification 8 we then create a temporary ground classification that we refer to as the “lowest ground”. For this we run lasground (README) with a set of suitable parameters that were found by experimentation on two of the most complex tiles from the center of the survey.

lasground -i Flaugergues\tiles_lowest\flauge*.laz ^
          -ignore_class 0 7 ^
          -step 5 -hyper_fine -bulge 1.5 -spike 0.5 ^
          -odir Flaugergues\tiles_lowest_ground -olaz ^
          -cores 4

We then “thicken” this “lowest ground” by classifying all points that are between 2 cm below and 15 cm above the lowest ground to a temporary classification code 6 using the lasheight (README) tool. Depending on the spread of points in your data set you may want to tighten this range accordingly, for example when processing the flightlines acquired by the Velodyne Puck individually. We picked our range based on the visual experiments with “drop lines” and “rise lines” in the lasview viewer that are shown in images above.

lasheight -i Flaugergues\tiles_lowest_ground\flauge*.laz ^
          -do_not_store_in_user_data ^
          -classify_between -0.02 0.15 6 ^
          -odir Flaugergues\tiles_thick_ground -olaz ^
          -cores 4

Lowest ground points triangulated. Other points colored by flightline ID.
Thick ground points in orange. Other points in black.
Zoom on lowest ground point triangulation. Other points colored by flightline ID.
Thick ground points in orange. Other points in black.

The final ground classification is obtained by creating the “median ground” from the “thick ground”. This uses a brand-new option in the lasthin (README) tool of LAStools. The new ‘-percentile 50 10’ option selects the point that is closest to the specified percentile of 50 of all point elevations within a grid cell of a specified size given there are at least 10 points in that cell. The selected point either survives the thinning operation or gets marked with a specified classification code or flag.

lasthin -i Flaugergues\tiles_thick_ground\flauge*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.1 -percentile 50 10 ^
        -classify_as 8 ^
        -odir Flaugergues\tiles_median_ground_10_10cm -olaz ^
        -cores 4

lasthin -i Flaugergues\tiles_median_ground_10_10cm\%NAME%*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.2 -percentile 50 10 ^
        -classify_as 8 ^
        -odir Flaugergues\tiles_median_ground_10_20cm -olaz ^
        -cores 4

lasthin -i Flaugergues\tiles_median_ground_10_20cm\%NAME%*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.4 -percentile 50 10 ^
        -classify_as 8 ^
        -odir Flaugergues\tiles_median_ground_10_40cm -olaz ^
        -cores 4

lasthin -i Flaugergues\tiles_median_ground_10_40cm\flauge*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.8 -percentile 50 10 ^
        -classify_as 8 ^
        -odir Flaugergues\tiles_median_ground_10_80cm -olaz ^
         -cores 4

We now compare a triangulation of the median ground points with a triangulation of the highest and the lowest points per 10 cm by 10 cm cell to demonstrate that – at least in open areas – we really have computed a median ground surface.

Triangulating median ground returns into a TIN.
Triangulating lowest returns on 10 cm by 10 cm grid into a TIN.
Triangulating highest returns on 10 cm by 10 cm grid into a TIN.
Comparing triangulation of median ground points to lowest returns on a 10 cm by 10 cm grid with 5 cm rise lines.
Comparing triangulation of median ground points to highest returns on a 10 cm by 10 cm grid with 5 cm drop lines.

Finally we raster the tiles with the las2dem (README) tool onto binary elevation grids in BIL format. Here we make the resolution dependent on the tile size, giving the 25 meter and 50 meter tiles the highest resolution of 10 cm and rasterize the 100 meter and 200 meter tiles at 20 cm and 40 cm respectively.

las2dem -i Flaugergues\tiles_median_ground_10_80cm\*_25.laz ^
        -i Flaugergues\tiles_median_ground_10_80cm\*_50.laz ^
        -keep_class 8 ^
        -step 0.1 -use_tile_bb ^
        -odir Flaugergues\tiles_dtm -obil ^
        -cores 4

las2dem -i Flaugergues\tiles_median_ground_10_80cm\*_100.laz ^
        -keep_class 8 ^
        -step 0.2 -use_tile_bb ^
        -odir Flaugergues\tiles_dtm -obil ^
        -cores 4

las2dem -i Flaugergues\tiles_median_ground_10_80cm\*_200.laz ^
        -keep_class 8 ^
        -step 0.4 -use_tile_bb ^
        -odir Flaugergues\tiles_dtm -obil ^
        -cores 4

Because all LAStools can read BIL files via on the fly conversion from rasters to points we can visually inspect the resulting elevation rasters with the lasview (README) tool. By adding the ‘-faf’ or ‘files_are_flightlines’ argument we treat the BIL files as if they were different flightlines which allows us to assign different color to points from different files to better inspect the transitions between tiles. The ‘-points 10000000’ argument instructs lasview to load up to 10 million points into memory instead of the default 5 million.

lasview -i Flaugergues\tiles_dtm\*.bil -faf ^
        -points 10000000

Final raster tiles in BIL format of three different sizes form seamless DTM.

For visual comparison we also produce a DSM and create hillshades. Note that the workflow for DSM creation shown below produces a “highest DSM” that will always be a few centimeter above the “median DTM”. This will be noticeable only in open areas of the terrain where the DSM and the DTM should coincide and their elevation should be identical.

highest DSM
median DTM 

lasthin -i Flaugergues\tiles_denoised\flauge*.laz ^
        -keep_z_above 110 ^
        -filtered_transform ^
        -set_classification 18 ^
        -ignore_class 7 18 ^
        -step 0.1 -highest ^
        -classify_as 5 ^
        -odir Flaugergues\tiles_highest -olaz ^
        -cores 4

las2dem -i Flaugergues\tiles_highest\*_25.laz ^
        -i Flaugergues\tiles_highest\*_50.laz ^
        -keep_class 5 ^
        -step 0.1 -use_tile_bb ^
        -odir Flaugergues\tiles_dsm -obil ^
        -cores 4

las2dem -i Flaugergues\tiles_highest\*_100.laz ^
        -keep_class 5 ^
        -step 0.2 -use_tile_bb ^
        -odir Flaugergues\tiles_dsm -obil ^
        -cores 4

las2dem -i Flaugergues\tiles_highest\*_200.laz ^
        -keep_class 5 ^
        -step 0.4 -use_tile_bb ^
        -odir Flaugergues\tiles_dsm -obil ^
        -cores 4

We thank YellowScan for challenging us to process their drone LiDAR with LAStools in order to present results at their LiDAR for Drone 2017 Conference and for sharing several example data sets with us, including the one used here.

Density and Spacing of LiDAR

Posted on by

Recently I worked with LiDAR from an Optech Gemini scanner in the Philippines and with LiDAR from a RIEGL Q680i scanner in Thailand. The two devices scan the terrain below with a very different pattern: the Optech uses an oscillating mirror producing zig-zag scan lines, whereas the RIEGL uses a rotating polygon producing parallel scan lines. In the following we investigate the point distribution in a small piece cut from each flightline.

First we compute an estimate for the average point density and the average point spacing with lasinfo:

D:\lastools\bin>lasinfo -i optech.laz ^
                        -nh -nv -nmm -cd
number of last returns: 2330671
covered area in square units/kilounits: 1280956/1.28
point density: all returns 2.25 last only 1.82 (per unit^2)
      spacing: all returns 0.67 last only 0.74 (in units)

The option ‘-nh’ asks lasinfo not to print the header (aka ‘no header’), the option ‘-nv’ means not to print the variable length records (aka ‘no vlr’), the option ‘-nmm’ supresses the output of minimum and maximum point values (aka ‘no min max’), and the option ‘-cd’ requests that an average point density is computed (aka ‘compute density’) from which the average point spacing is derived.

D:\lastools\bin>lasinfo -i riegl.laz ^
                        -nh -nv -nmm -cd
number of last returns: 3707217
covered area in square units/kilounits: 889816/0.89
point density: all returns 4.58 last only 4.17 (per unit^2)
      spacing: all returns 0.47 last only 0.49 (in units)

What we really want to know are pulse density and pulse spacing which we get by only counting one return from every pulse. Commonly one uses the last return, which is reported as ‘last only’ values by lasinfo. According to lasinfo the Optech and the RIEGL scan have an average pulse density of 1.82 and 4.17 [pulses per square meter] and an average pulse spacing of 0.74 and 0.49 [meters] respectively.

Single number averages do not capture anything about the actual distribution of the last returns in the swath of the scan. Let us compute density rasters with lasgrid. Because the density of the two scans is different we use 3 by 3 meter cells for the Optech, which should average 3 * 3 * 1.82 = 16.38 points per cell, and 2 by 2 meter cells for the RIEGL, which should average 2 * 2 * 4.17 = 16.68 points per cell.

lasgrid -i optech.laz -last_only ^
        -density -step 3 ^
        -odix _density3x3 -obil
lasview -i optech_density3x3.bil

lasgrid -i riegl.laz -last_only ^
        -density -step 2 ^
        -odix _density2x2 -obil
lasview -i riegl_density2x2.bil

A visual comparison of the two resulting density rasters in BIL format with lasview shows first differences in pulse distribution for the two scanners.

Using 'lasview' to inspect a BIL density raster: Elevation and colors show number of points per 3 by 3 meter cell for Optech (top) and per 2 by 2 meter cell for RIEGL (bottom).

Using ‘lasview’ to inspect a BIL density raster: Elevations and colors illustrate the number of last returns per 3 by 3 meter cell for Optech (top) and per 2 by 2 meter cell for RIEGL (bottom). The images use differently scaled color ramps.

It is noticable that the number of last returns per cell increases at the edges of the swath for the Optech whereas it decreases for the RIEGL. For the Optech it is higher due to the slowing down of the oscillating mirror and the reversing of scan direction when reaching either side of the scanline. This decreases the distances between pulses at the edges of the zig-zag scan and increases the number of last returns falling into cells there. For the rotating polygon scan of the RIEGL these pulse distances are more uniform. Here the numbers are lower because many cells along the edges of the swath are only partly covered by the scan and therefore receive fewer last returns.

lasgrid -i optech.laz -last_only ^
        -density -step 3 ^
        -false -set_min_max 0 20 ^
        -odix _density3x3 -opng

lasgrid -i riegl.laz -last_only ^
        -density -step 2 ^
        -false -set_min_max 0 20 ^
        -odix _density2x2 -opng

Here another visualization of the pulse distribution by letting lasgrid false-color density rasters to a fixed range of 0 to 20. This range is based on the near identical expected values of around 16.38 and 16.68 last returns per cell based on the average densities that lasinfo reported and the cell sized used here.

A more quantitative check of how well distributed the pulse are is to generate histograms. You can use lasinfo with the ‘-histo z 1’ option to print a histogram of the z values of the BIL rasters and import these statistics into your favorite software to create a chart.

lasinfo -i optech_density3x3.bil ^
        -nh -nv -nmm -histo z 1
lasinfo -i riegl_density2x2.bil ^
        -nh -nv -nmm -histo z 1

Both distributions have their peaks at the expected 16 to 17 last returns per cell although the Optech distribution has a significantly wider spead. The odd two-peak distribution of the Optech scan seems puzzling at first but looking at the density image shows that this is an aliasing artifact of putting a fixed 3 by 3 meter grid over zig-zagging scanlines.

The next experiments uses the new ‘-edge_shortest’ and ‘-edge_longest’ options available in las2dem that compute for each point the length of its shortest or longest edge in a Delaunay triangulation and then rasters these lengths. Below you see the command lines to do this and generate histograms of the rastered edge lengths.

las2dem -i optech.laz -last_only ^
        -step 1.5 ^
        -edge_longest ^
        -odix _edge_longest -obil
lasinfo -i optech_edge_longest.bil ^
        -nh -nv -nmm -histo z 0.1 ^
        -o optech_edge_shortest_0_10.txt
las2dem -i optech.laz -last_only ^
        -step 1.5 ^
        -edge_shortest ^
        -odix _edge_shortest -obil
lasinfo -i optech_edge_shortest.bil ^
        -nh -nv -nmm -histo z 0.05 ^
        -o optech_edge_shortest_0_05.txt
las2dem -i riegl.laz -last_only ^
        -edge_longest ^
        -odix _edge_longest -obil
lasinfo -i riegl_edge_longest.bil ^
        -nh -nv -nmm -histo z 0.05 ^
        -o riegl_edge_longest_0_05.txt
las2dem -i riegl.laz -last_only ^
        -edge_shortest ^
        -odix _edge_shortest -obil
lasinfo -i riegl_edge_shortest.bil ^
        -nh -nv -nmm -histo z 0.05 ^
        -o riegl_edge_shortest_0_05.txt

These resulting histograms do look quite different.

The more important histograms are those of longest edge lengths. They illustrate the observed spacing between pulses showing us the maximal distance of each pulse from its surrounding pulses. A perfect pulse distribution that samples the ground evenly in all directions would have one narrow peak. The Optech scan is far from this ideal with a wide and flat peak that tells us that pulses are spaced apart anywhere from 1.2 to 2.2 meters. Visualizing the scan pattern shows that the 2.2 meter spacings happen at the tips of zig-zag and the 1.2 meter spacings at nadir. The RIEGL scan has much narrower peak with most pulses spaced apart at most 60 to 80 centimeters.

The histograms of shortest edge lengths illustrate how close pulses are spaced. Again, it helps to visualize the zig-zag scan pattern to explain the odd distribution of shortest edge lengths in the Optech scan: as we get closer to the edges of the flightline the pulse spacing along the scanline becomes increasingly dense. This is represented on the left side in the histogram that is slowly leveling off to zero. Again, the RIEGL scan has a narrower peak with most pulses spaced no closer than 35 to 55 centimeters.

To confirm our findings we illustrate where the scanners produce longest or shortest edges using ranges that we find in the histograms above and create false-colored rasters:

las2dem -i optech.laz -last_only ^
        -step 1.5 -edge_longest ^
        -false -set_min_max 1 2.5 ^
        -odix _edge_longest -opng
las2dem -i optech.laz -last_only ^
        -step 1.5 -edge_shortest ^
        -false -set_min_max 0 0.7 ^
        -odix _edge_shortest -opng
las2dem -i riegl.laz -last_only ^
        -edge_longest ^
        -false -set_min_max 0.5 1.0 ^
        -odix _edge_longest -opng
las2dem -i riegl.laz -last_only ^
        -edge_shortest ^
        -false -set_min_max 0.1 0.7 ^
        -odix _edge_shortest -opng

The color-codings clearly illustrate that the Optech has a much wider pulse spacing on both edges of the flightline but also a much narrower pulse spacing. That reads contradictory but one are the “within-zig-zag” spacings and the other are the “between-zig-zag” spacings. The RIEGL has an overall much more even distribution in pulse spacings.

As a final confirmation to what causes the pulse spacings to be both widest and narrowest at the edges of the flightline for the Optech scan we scrutinizing the distribution of last returns and their triangulation visually.

Now it should be clear. At the edge of the flightline the Optech scanner has increasingly close-spaced pulses within each zig-zagging scan line but also increasingly distant-spaced pulses between subsequent pairs of zig-zagging scanlines. At the very edge of the flightline the scanner is almost sampling a “linear area” twice but leaves unsampled gaps that are twice as wide as in the center of the flightline. I believe this is one of the reasons why people “cut off” the edges of the flightlines based on the scan angle when using oscillating mirror systems. Now let’s took at the RIEGL.

There is no obvious difference in pulse distribution at the edges and the center of the flightline. The entire width of the swath seems more or less evenly sampled.

The above results show that “average point density” and “average point spacing” do not capture the whole picture. A quality check that properly verifies that a scan has the desired point / pulse density and spacing should measure the actual spacings – especially when operating a scanner with oscillating mirrors. We have done it here by computing for each last return its longest surrounding edge in the Delaunay triangulation, rasterizing these edge length, and then generating a histogram of the raster values.

 

Good follow-up reads are Dr. Ullrich’s “Impact of point distribution on information content of point clouds of airborne LiDAR” presented at ELMF 2013 and “Assessing the Information Content of LiDAR Point Clouds” from PhoWo 2013.

clone wars and drone fights

Posted on by

NEWSFLASH: update on Feb 5th and 6th and 17th (see end of article)

The year 2014 shapes to be a fun one in which the LiDAR community will see some major laser battles. (-: First, ESRI starts a “lazer clone war” with the open source community saddening your very own LAS clown here at rapidlasso with a proprietary LAZ clone. Then AHAB and Optech start to squirmish over the true penetration depth of their bathymetric systems in various forums. And now RIEGL – just having launched their Q1560 “crossfire” for battling Leica and Optech in the higher skies – is heading into an all-out “laser drone war” with FARO and Velodyne over arming UAVs and gyrocopters with lasers by rolling out their new (leaked?) RIEGL VUX-1 scanner … (-;

Some specs below, but no details yet on power-consumption.

  • about 3.85 kilograms
  • 225mm x 180mm x 125mm
  • survey-grade accuracy of 25mm
  • echo signal digitisation
  • online waveform processing
  • measurement rate up to 500 kHz
  • field of view 300 degrees
  • internal 360 GB SSD storage or real-time data via LAN-TCP/IP
  • collects data at altitudes or ranges of more than 1,000 ft
  • will be on display during ILMF 2014 (Feb 17-19, Denver, USA)

The new RIEGL VUX-1 laser scanner for UAVs ...

The new RIEGL VUX-1 laser scanner for UAVs

Below a proof-of-concept concept video by Phoenix Aerial on how such a system (here based on a Velodyne HDL-32E scanner) may look and operate.

A similar system introduced by Airbotix is shown below.

Aibot X6 with laser scanner from Velodyne HDL-32E

Aibot X6 with laser scanner from Velodyne HDL-32E

But before trusting these kind of drones with your expensive hardware you should make plenty of test flights with heavy payloads. Here a suggestion to make this more fun:

UPDATE (February 5th): It is reported by SPAR Point Group that the RIEGL VUX-1 LiDAR scanner for drones is to be shown at International LiDAR Mapping Forum in Denver later this month. Given the ceremony that the Q1560 “crossfire” was introduced with, it surely seems as if Geospatial World Magazine spilled the news early …

UPDATE (February 6th): RIEGL officially press-released the VUX-1 LiDAR. Its weight is below 4kg and its range up to 1000 ft …

UPDATE (February 17th): Today the new VUX-1 finally made it’s much anticipated world-premiere during ILMF 2014 in Denver. For more images of the official unveiling see RIEGL’s blog.

The VUX-1 LiDAR scanner for UAVs unveiled at ILMF 2014.

The VUX-1 LiDAR scanner for UAVs unveiled at ILMF 2014.

new scanner cross examines terrain

Posted on by

It seems Leica and Optech may have come a bit under “crossfire” (-: by the other big news (besides adding PulseWaves support to RiPROCESS) at the RIEGL LiDAR 2013 user conference in Vienna, namely the unveiling of the new LMS-Q1560 airborne laser scanner. Below you see the management team doing the ceremonial act.Image

The specification looks great and will especially make those happy looking for dense surveys in complex terrain or for more LiDAR returns from building facades in urban environments. By employing two units that “crossfire” at a particular angle while one unit looks forward and one backward, the combined system can provide a better point spacing on the ground and a better point coverage on vertical surfaces.

The new system combines two LMS-Q780 units each firing at rates of up 400 kHz via a shared rotating mirror in the center. The scan lines of the two laser-beam-emitting units are not perpendicular to the flight direction but angled at +12 and -12 degrees (actual number may be different) in the x-y plane to assure that they independently and consistently sample the covered terrain – unaffected from flying height or aircraft speed. Furthermore, the two fanning sheets of laser pulses are angled at +7 and -7 degrees (actual number may be different) in flight direction, meaning that one unit is forward-looking and the other backward-looking. This increases the likelihood that the sides of buildings are scanned as well. Colloquially put, those back-facing facades perpendicular to the flight direction that are not seen by the first round of laser pulses from the forward-looking unit will be hit by the second round of laser pulses from the  backward-looking unit and vice-versa. The above details may not be 100% correct but are what I (wrongly?) remember from the technical presentation that Dr. Ullrich gave (see below).

Image

The overall ability to sample so rapidly is realized by combining the “crossfire” that I have just described with the “rapidfire” of the earlier Q780 units. Each unit is able to operate with up to 10 pulses simultaneously in the air by subsequently resolving ambiguities off-line with the multiple-time-around (MTA) processing technique introduced for earlier models (see the product descriptions for the Q680i or the Q780 for more on that topic).

Hence, this new scanner with “crossfire technology” promises to give your terrain a dense cross examination … and yes, it is a full waveform scanner. (-: