LASmoons: Martin Buchauer

Martin Buchauer (recipient of three LASmoons)
Cartography & Geomedia Technology
University of Applied Science Munich, GERMANY

Background:
Salt marsh areas provide numerous services such as natural flood defenses, carbon sequestration, agricultural services, and are a valuable coastal habitat for flora, fauna and humans. However, they are not only threatened by the constant rise of sea levels caused by global warming but also by human settlement in coastal areas. A sensible local coastal development is important as it may help to support the development and progression of stressed salt marshes.

Looking South you can see the salt marsh area next to a famous golf course with St Andrews in the background.

Goal:

This research aims to visualize and extract vegetation metrics as well as the temporal analysis of four salt marsh data sets which are derived from terrestrial laser scanning. Located at the South and North shore of the Eden Estuary near St Andrews, Scotland, the scans were acquired in the summer and winter of 2016. Ground based laser scanning is an ideal method of fully reconstructing vegetation structures as well as having the ability to retrieve meaningful metrics such as height, area, and vegetation density. Although this technology has frequently been applied in the area of forestry, its application to salt marsh areas has not yet fully explored.

Data:
+
 TLS data acquired with a Leica HDS6100 (average density of 38000 points/m²)
+ ground control points (field data)

LAStools processing:
1) check the quality of the LiDAR data [lasinfo, lasoverlap, lasgrid]
2) merge and retile the original data with buffers [lastile]
3) classify point clouds into ground and non-ground [lasthin, lasground]
4) create digital terrain (DTM) and digital surface models (DSM) [lasthin, las2dem, blast2dem]

First Look with LAStools at LiDAR from Hovermap Drone by CSIRO

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

 

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.

Removing Low Noise from RIEGL’s VUX-1 UAV LiDAR flown in the Philippines

In this tutorial we are removing some “tricky” low noise from LiDAR point clouds in order to produce a high-resolution Digital Terrain Model (DTM). The data was flown above a tropical beach and mangrove area in the Philippines using a VUX-1 UAV based system from RIEGL mounted on a helicopter. The survey was done as a test flight by AB Surveying who have the capacity to fly such missions in the Philippines and who have allowed us to share this data with you for educational purposes. You can download the data (1 GB) here. It covers a popular twin beach knows as “Nacpan” near El Nido in Palawan (that we happen to have visited in 2014).

A typical beach fringed by coconut palms in Palawan, Philippines.

We start our usual quality check with a run of lasinfo. We add the ‘-cd’ switch to compute an average point density and the ‘-histo gps_time 1’ switch to produce a 1 second histogram for the GPS time stamps.

lasinfo -i lalutaya.laz ^
        -cd ^
        -histo gps_time 1 ^
        -odix _info -otxt

You can download the resulting lasinfo report here. It tells us that there are 118,740,310 points of type 3 (with RGB colors) with an average density of 57 last returns per square meter. The point coordinates are in the “PRS92 / Philippines 1” projection with EPSG code 3121 that is based on the “Clarke 1866” ellipsoid.

Datum Transform

We prefer to work in an UTM projection based on the “WGS 1984” ellipsoid, so we will first perform a datum transform based on the seven parameter Helmert transformation – a capacity that was recently added to LAStools. For this we first need a transform to get to geocentric or Earth-Centered Earth-Fixed (ECEF) coordinates on the current “Clarke 1866” ellipsoid, then we apply the Helmert transformation that operates on geocentric coordinates and whose parameters are listed in the TOWGS84 string of EPSG code 3121 to get to geocentric or ECEF coordinates on the “WGS 1984” ellipsoid. Finally we can convert the coordinates to the respective UTM zone. These three calls to las2las accomplish this.

las2las -i lalutaya.laz ^
        -remove_all_vlrs ^
        -epsg 3121 ^
        -target_ecef ^
        -odix _ecef_clark1866 -olaz

las2las -i lalutaya_ecef_clark1866.laz ^
        -transform_helmert -127.62,-67.24,-47.04,-3.068,4.903,1.578,-1.06 ^
        -wgs84 -ecef ^
        -ocut 10 -odix _wgs84 -olaz
 
las2las -i lalutaya_ecef_wgs84.laz ^
        -target_utm auto ^
        -ocut 11 -odix _utm -olaz

In these steps we implicitly reduced the resolution that each coordinate was stored with from quarter-millimeters (i.e. scale factors of 0.00025) to the default of centimeters (i.e. scale factors of 0.01) which should be sufficient for subsequent vegetation analysis. The LAZ files also compress better when coordinates a lower resolution so that the ‘lalutaya_utm.laz’ file is over 200 MB smaller than the original ‘lalutaya.laz’ file. The quantization error this introduces is probably still below the overall scanning error of this helicopter survey.

Flightline Recovery

Playing back the file visually with lasview suggests that it contains more than one flightline. Unfortunately the point source ID field of the file is not properly populated with flightline information. However, when scrutinizing the GPS time stamp histogram in the lasinfo report we can see an occasional gap. We highlight two of these gaps in red between GPS second 540226 and 540272 and GPS second 540497 and 540556 in this excerpt from the lasinfo report:

gps_time histogram with bin size 1
[...]
 bin 540224 has 125768
 bin 540225 has 116372
 bin 540226 has 2707
 bin 540272 has 159429
 bin 540273 has 272049
 bin 540274 has 280237
[...]
 bin 540495 has 187103
 bin 540496 has 180421
 bin 540497 has 126835
 bin 540556 has 228503
 bin 540557 has 275025
 bin 540558 has 273861
[...]

We can use lasplit to recover the original flightlines based on gaps in the continuity of GPS time stamps that are bigger than 10 seconds:

lassplit -i lalutaya_utm.laz ^
         -recover_flightlines_interval 10 ^
         -odir strips_raw -o lalutaya.laz

This operation splits the points into 11 separate flightlines. The points within each flightline are stored in the order that the vendor software – which was RiPROCESS 1.7.2 from RIEGL according to the lasinfo report – had written them to file. We can use lassort to bring them back into the order they were acquired in by sorting first on the GPS time stamp and then on the return number field:

lassort -i strips_raw\*.laz ^
        -gps_time -return_number ^
        -odir strips_sorted -olaz ^
        -cores 4

Now we turn the sorted flightlines into tiles (with buffers !!!) for further processing. We also erase the current classification of the data into ground (2) and medium vegetation (4) as a quick visual inspection with lasview immediately shows that those are not correct:

lastile -i strips_sorted\*.laz ^
        -files_are_flightlines ^
        -set_classification 0 ^
        -tile_size 250 -buffer 30 -flag_as_withheld ^
        -odir tiles_raw -o lalu.laz

Quality Checking

Next comes the standard check of flightline overlap and alignment check with lasoverlap. Once more it become clear why it is so important to have flightline information. Without we may have missed what we are about to notice. We create false color images load into Google Earth to visually assess the situation. We map all absolute differences between flightlines below 5 cm to white and all absolute differences above 30 cm to saturated red (positive) or blue (negative) with a gradual shading from white to red or blue for any differences in between. We also create an overview KML file that lets us quickly see in which tile we can find the points for a particular area of interest with lasboundary.

lasoverlap -i tiles_raw\*.laz ^
           -step 1 -min_diff 0.05 -max_diff 0.30 ^
           -odir quality -opng ^
           -cores 4

lasboundary -i tiles_raw\*.laz ^
            -use_tile_bb -overview -labels ^
            -o quality\overview.kml

The resulting visualizations show (a) that our datum transform to the WGS84 ellipsoid worked because the imagery aligns nicely with Google Earth and (b) that there are several issues in the data that require further scrutiny.

In general the data seems well aligned (most open areas are white) but there are blue and red lines crossing the survey area. With lasview have a closer look at the visible blue lines running along the beach in tile ‘lalu_765000_1252750.laz’ by repeatedly pressing ‘x’ to select a different subset and ‘x’ again to view this subset up close while pressing ‘c’ to color it differently:

lasview -i tiles_raw\lalu_765000_1252750.laz

These lines of erroneous points do not only happen along the beach but also in the middle of and below the vegetation as can be seen below:

Our initial hope was to use the higher than usual intensity of these erroneous points to reclassify them to some classification code that we would them exclude from further processing. Visually we found that a reasonable cut-off value for this tile would be an intensity above 35000:

lasview -i tiles_raw\lalu_765000_1252750.laz ^
        -keep_intensity_above 35000 ^
        -filtered_transform ^
        -set_classification 23

However, while this method seems successful on the tile shown above it fails miserably on others such as ‘lalu_764250_1251500.laz’ where large parts of the beach are very reflective and result in high intensity returns to to the dry white sand:

lasview -i tiles_raw\lalu_764250_1251500.laz ^
        -keep_intensity_above 35000 ^
        -filtered_transform ^
        -set_classification 23

Low Noise Removal

In the following we describe a workflow that can remove the erroneous points below the ground so that we can at least construct a high-quality DTM from the data. This will not, however, remove the erroneous points above the ground so a subsequent vegetation analysis would still be affected. Our approach is based on two obervations (a) the erroneous points affect only a relatively small area and (b) different flightlines have their erroneous points in different areas. The idea is to compute a set of coarser ground points separately for each flightline and – when combining them in the end – to pick higher ground points over lower ones. The combined points should then define a surface that is above the erroneous below-ground points so that we can mark them with lasheight as not to be used for the actual ground classification done thereafter.

The new huge_las_file_extract_ground_points_only.bat example batch script that you can download here does all the work needed to compute a set of coarser ground points for each flightline. Simply edit the file such that the LAStools variable points to your LAStools\bin folder and rename it to end with the *.bat extension. Then run:

huge_las_file_extract_ground_points_only strips_sorted\lalutaya_0000001.laz strips_ground_only\lalutaya_0000001.laz
huge_las_file_extract_ground_points_only strips_sorted\lalutaya_0000002.laz strips_ground_only\lalutaya_0000001.laz
huge_las_file_extract_ground_points_only strips_sorted\lalutaya_0000003.laz strips_ground_only\lalutaya_0000001.laz
...
huge_las_file_extract_ground_points_only strips_sorted\lalutaya_0000009.laz strips_ground_only\lalutaya_0000009.laz
huge_las_file_extract_ground_points_only strips_sorted\lalutaya_0000010.laz strips_ground_only\lalutaya_0000010.laz
huge_las_file_extract_ground_points_only strips_sorted\lalutaya_0000011.laz strips_ground_only\lalutaya_0000011.laz

The details on how this batch script works – a pretty standard tile-based multi-core processing workflow – are given as comments in this batch script. Now we have a set of individual ground points computed separately for each flightline and some will include erroneous points below the ground that the lasground algorithm by its very nature is likely to latch on to as you can see here:

The trick is now to utilize the redundancy of multiple scans per area and – when combining flightlines – to pick higher rather than lower ground points in overlap areas by using the ground point closest to the 75th elevation percentile per 2 meter by 2 meter area with at least 3 or more points with lasthin:

lasthin -i strips_ground_only\*.laz -merged ^
        -step 2 -percentile 75 3 ^
        -o lalutaya_ground_only_2m_75_3.laz

There are still some non-ground points in the result as ground-classifying of flightlines individually often results in vegetation returns being included in sparse areas along the edges of the flight lines but we can easily get rid of those:

lasground_new -i lalutaya__ground_only_2m_75_3.laz ^
              -town -hyper_fine ^
              -odix _g -olaz

We sort the remaining ground points into a space-filling curve order with lassort and spatially index them with lasindex so they can be efficiently accessed by lasheight in the next step.

lassort -i lalutaya__ground_only_2m_75_3_g.laz ^
        -keep_class 2 ^
        -o lalutaya_ground.laz

lasindex -i lalutaya_ground.laz

Finally we have the means to robustly remove the erroneous points below the ground from all tiles. We use lasheight with the ground points we’ve just so painstakingly computed to classify all points 20 cm or more below the ground surface they define into classification code 23. Later we simply can ignore this classification code during processing:

lasheight -i tiles_raw\*.laz ^
          -ground_points lalutaya_ground.laz ^
          -do_not_store_in_user_data ^
          -classify_below -0.2 23 ^
          -odir tiles_cleaned -olaz ^
          -cores 4

Rather than trying to ground classify all remaining points we run lasground on a thinned subset of all points. For this we mark the lowest point in every 20 cm by 20 cm grid cell with some temporary classification code such as 6.

lasthin -i tiles_cleaned\*.laz ^
        -ignore_class 23 ^
        -step 0.20 -lowest -classify_as 6 ^
        -odir tiles_thinned -olaz ^
        -cores 4

Finally we can run lasground to compute the ground classification considering all points with classification code 6 by ignoring all points with classification codes 23 and 0.

lasground_new -i tiles_thinned\*.laz ^
              -ignore_class 23 0 ^
              -city -hyper_fine ^
              -odir tiles_ground_new -olaz ^
              -cores 4

And finally we can create a DTM with a resolution of 25 cm using las2dem and the result is truly beautiful:

las2dem -i tiles_ground_new\*.laz ^
        -keep_class 2 ^
        -step 0.25 -use_tile_bb ^
        -odir tiles_dtm_25cm -obil ^
        -cores 4

We have to admit that a few bumps are left (see mouse cursor below) but adjusting the parameters presented here is left as an exercise to the reader.

We would again like to acknowledge AB Surveying whose generosity has made this blog article possible. They have the capacity to fly such missions in the Philippines and who have allowed us to share this data with you for educational purposes.

LASmoons: Sebastian Kasanmascheff

Sebastian Kasanmascheff (recipient of three LASmoons)
Forest Inventory and Remote Sensing
Georg-August-Universität Göttingen, GERMANY

Background:
Forest inventories are the backbone of forest management in Germany. In most federal forestry administrations in Germany, they are performed every ten years in order to assure that logging activities are sustainable. The process involves trained foresters who visit each stand (i.e. an area where the forest is similar in terms of age structure and tree species) and perform angle count sampling as developed by Walter Bitterlich in 1984. In a second step the annual growth is calculated using yield tables and finally a harvest volume is derived. There are three particular reasons to investigate how remote sensing can be integrated in the current inventory system:

  1. The current process does not involve random sampling of the sampling points and thus does not offer any measure of the accuracy of the data.
  2. Forest engineers hardly ever rely on the inventory data as a stand-alone basis for logging planning. Most often they rely on intuition alone and on the total volume count that they have to deliver for a wider area every year.
  3. In the last ten years, the collection of high-resolution LiDAR data has become more cost-effective and most federal agencies in Germany have access to it.

In order to be able to integrate the available remote-sensing data for forest inventories in Germany, it is important to tell apart different tree species as well as estimate their volumes.

Hesse is one of the most forested federal states in Germany.

Goal:
The goal of this project is to perform an object-based classification of conifer trees in Northern Hesse based on high-resolution LiDAR and multi-spectral orthophotos. The first step is to delineate the tree crowns. The second step is to perform a semi-automated classification using the spectral signature of the different conifer species.

Data:
+
 DSM (1m), DTM (1m), DSM (0.2 m) of the study area
+ Stereo images with 0.2 m resolution
+ high-resolution LiDAR data (average 10 points/m²)
+ forest inventory data
+ vector files of the individual forest stands
+ ground control points (field data)
All of this data is provided by the Hessian Forest Agency (HessenForst).

LAStools processing:
1) merge and clip the LAZ files [las2las]
2) classify ground and non-ground points [lasground]
3) remove low and high outliers [lasheight, lasnoise]
4) identify buildings within the study area [lasclassify]
5) create a normalized point cloud [lasheight]
6) create a highest-return canopy height model (CHM) [lasthin, las2dem]
7) create a pit-free (CHM) with the spike-free algorithm [las2dem]

LASmoons: Manuel Jurado

Manuel Jurado (recipient of three LASmoons)
Departamento de Ingeniería Topográfica y Cartografía
Universidad Politécnica de Madrid, SPAIN

Background:
The availability of LiDAR data is creating a lot of innovative possibilities in different fields of science, education, and other field of interests. One of the areas that has been deeply impacted by LiDAR is cartography and in particular one highly specialized sub-field of cartography in the domain of recreational and professional orienteering running: the production of high-quality maps for orienteering races (Ditz et al., 2014). These are thematic maps with a lot of fine detail which demands many hours of field work for the map maker. In order to reduce the fieldwork, LiDAR data obtained from Airborne Research Australia (ARA) is going to be used in order to obtain DEM and to extract features that must be included in these maps. The data will be filtered and processed with the help of LAStools.

Final map with symbolism typical for use in orienteering running

Goal:
The goal of this project is to extract either point (boulders, mounds), linear (contours, erosion gullies, cliffs) and area features (vegetation density) that should be drawn in a orienteering map derived from high-resolution LiDAR. Different LiDAR derived raster images are being created: 0.5m DTM, vegetation density (J. Ryyppo, 2013), slope, Sky-View factor (Ž. Kokalj et al., 2011), and shaded relief. The area used is in Renmark, South Australia and the produced map is going to be used for the Australian Orienteering Championships 2018.

Sky-View factor of DTM for same area as shown above.

Data:
+
4 square kilometers of airborne LiDAR data produced by Airborne Research Australia at 18 pulses per square meter using the full waveform scanning LiDAR Q680i-S laser scanner from RIEGL
+ 60 hours of check and validation work in the field

LAStools processing:
1) tile into 500 by 500 meter tiles with 20 meter buffer [lastile]
2) classify isolated points as noise [lasnoise]
3) classify point clouds into ground and non-ground [lasground]
4) create a Digital Terrain Model (DTM) [las2dem]
5) normalize height of points above the ground [lasheight]
6) compute vegetation density metrics [lascanopy]
7) create hillshades of the raster DTMs [blast2dem or GDAL]

References:
Ditz, Robert, Franz Glaner, and Georg Gartner. (2014). “Laser Scanning and Orienteering Maps.” Scientific Journal of Orienteering 19.1.
JRyyppo, Jarkko. (2013). “Karttapullautin vegetation mapping guide”.
Kokalj, Žiga, Zaksek, Klemen, and Oštir, Krištof. (2011). Application of sky-view factor for the visualization of historic landscape features in lidar-derived relief models. Antiquity. 85. 263-273.

LASmoons: Chris J. Chandler

Chris J. Chandler (recipient of three LASmoons)
School of Geography
University of Nottingham, UNITED KINGDOM

Background:
Wetlands provide a range of important ecosystem services: they store carbon, regulate greenhouse gas emissions, provide flood protection as well as water storage and purification. Preserving these services is critical to achieve sustainable environmental management. Currently, mangrove forests are protected in Mexico, however, fresh water wetland forests, which also have high capacity for storing carbon both in the trees and in the soil, are not protected under present legislation. As a result, coastal wetlands in Mexico are threatened by conversion to grazing areas, drainage for urban development and pollution. Given these threats, there is an urgent need to understand the current state and distribution of wetlands to inform policy and protect the ecosystem services provided by these wetlands.
In this project we will combine field data collection, satellite data (i.e. optical remote sensing, radar and LiDAR remote sensing) and modelling to provide an integrated technology for assessing the value of a range of ecosystem services, tested to proof of concept stage based on carbon storage. The outcome of the project will be a tool for mapping the value of a range of ecosystem services. These maps will be made directly available to local stakeholders including policy makers and land users to inform policy regarding forest protection/legislation and aid development of financial incentives for local communities to protect these services.

Wetland classification in the Chiapas region of Mexico

Goal:
At this stage of the project we have characterized wetlands for three priority areas in Mexico (Pantanos de Centla, La Encrucijada and La Mancha). Next stage is the up scaling of the field data at the three study sites using LiDAR data for producing high quality Canopy Height Model (CHM), which has been of great importance for biomass estimation (Ferraz et al., 2016). A high quality CHM will be achieved using LAStools software.

Data:
+
LiDAR provided by the Mexican National Institute of Statistics and Geography (INEGI)
+ average height: 5500 m, mirror angle: +/- 30 degrees, speed: 190 knots
+ collected with Cessna 441, Conquest II system at 1 pts/m².

LAStools processing:
1)
create 1000 meter tiles with 35 meter buffer to avoid edge artifacts [lastile]
2) classify point clouds into ground and non-ground [lasground]
3) normalize height of points above the ground [lasheight]
4) create a Digital Terrain and Surface Model (DTM and DSM) [las2dem]
5) generate a spike-free Canopy Height Model (CHM) as described here and here [las2dem]
6) compute various metrics for each plot and the normalized tiles [lascanopy]

References:
Ferraz, A., Saatchi, S., Mallet, C., Jacquemoud S., Gonçalves G., Silva C.A., Soares P., Tomé, M. and Pereira, L. (2016). Airborne Lidar Estimation of Aboveground Forest Biomass in the Absence of Field Inventory. Remote Sensing, 8(8), 653.

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

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).

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.

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

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

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.

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.

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.