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.

LASmoons: Stéphane Henriod

Stéphane Henriod (recipient of three LASmoons)
National Statistical Committee of the Kyrgyz Republic
Bishkek, Kyrgyzstan

This pilot study is part of the International Climate Initiative project called “Ecosystem based Adaptation to Climate change in the high mountainous regions of Central Asia” that is funded by the Federal Ministry for the Environment, Nature Conservation, Building and Nuclear Safety (BMU) of Germany and implemented by the Deutsche Gesellschaft für Internationale Zusammenarbeit (GIZ) GmbH in Kyrgyzstan, Tajikistan and Kazakhstan.


The ecosystems in high mountainous regions of Central Asia are characterized by a unique diversity of flora and fauna. In addition, they are the foundation of the livelihoods of the local population. Specific benefits include clean water, pasture, forest products, protection against floods and landslides, maintenance of soil fertility, and ecotourism. However, the consequences of climate change such as melting glaciers, changing river runoff regimes, and weather anomalies including sharp temperature fluctuations and non-typical precipitation result in negative impacts on these ecosystems. Coupled with unwise land use, these events damage fragile mountain ecosystems and reduce their regeneration ability undermining the local population’s livelihoods. Therefore, people living in rural areas and directly depending on natural resources must adapt to adverse impacts of climate change. This can be done through a set of measures, known in the world practice as ecosystem-based adaptation (EbA) approach. It promotes the sustainable use of natural resources to sustain and enhance the livelihood of the population depending on those resources.

lasmoons_Stephane_Henriod_2 Goal:
In two selected pilot regions of Kyrgyzstan and Tajikistan, planned measures will concentrate on climate-informed management of ecosystems in order to maintain their services for the rural population. EbA always starts with identifying the vulnerability of the local population. Besides analyzing the socio-economic situation of the local population, this includes (1) assessing the ecological conditions of the ecosystems in the watershed and the related ecosystem services people benefit from, (2) identifying potential disaster risks, and (3) analyzing glacier dynamics with its response to water runoff. As a baseline to achieve this and to get spatially explicit, remote sensing based techniques and mapping activities need to be utilized.

A first UAV (unmanned aerial vehicle) mission has taken place in the Darjomj watershed of the Bartang valley in July 2016. RGB-NIR images as well as a high-resolution Digital Surface Model have been produced that now need to be segmented and analysed in order to produce comprehensive information. The main processing that will take advantage of LAStools is the generation of a DTM from the DSM that will then be used for identifying risk areas (flood zones, landslides and avalanches, etc.). The results of this approach will ultimately be compared with lower-cost satellite images (RapidEye, Planet, Sentinel).

+ High-resolution RGB and NIR image (10 cm) from a SenseFly Ebee
+ High-resolution DSM (10 cm) from a SenseFly Ebee

LAStools processing:
classify DSM points obtained via dense-matching photogrammetry into a SenseFly Ebee imagery into ground and non-ground points via processing pipelines as described here and here [lastile, lassort, lasnoise, lasground]
2) create a DTM [las2dem, lasgrid, blast2dem]
3) produce 3D visualisations to facilitate the communication around adaptation to climate change [lasview]

Removing low noise from Semi-Global Matching (SGM) Points

At PhoWo and INTERGEO 2015 rapidlasso was spending quality time with VisionMap who make the A3 Edge camera that provides fine resolution images from high altitudes and can quickly cover large areas. Under the hood of their LightSpeed software is the SURE dense-matching algorithm from nframes that turns images into photogrammetric point clouds. We were asked whether LAStools is able to create bare-earth DTM rasters from such points. If you have read our first, second, or third blog post on the topic you know that our asnwer was a resounding “YES!”. But we ran into an issue due to the large amount of low noise. Maybe the narrow angle between images at a high flying altitude affects the semi-global matching (SGM) algorithm. Either way, in the following we show how we use lascanopy and lasheight to mark low points as noise in a preprocessing step.

We obtained a USB stick containing a 2.42 GB file called “valparaiso_DSM_SURE_100.las” containing about 100 million points spaced 10 cm apart generated by SURE and stored with an (unnecessary high) resolution of millimeters (aka “resolution fluff”) as the third digit of all coordinates was always zero:

las2txt -i F:\valparaiso_DSM_SURE_100.las -stdout | more
255991.440 6339659.230 89.270
255991.540 6339659.240 89.270
255991.640 6339659.240 88.660
255991.740 6339659.230 88.730

We first compressed the bulky 2.42 GB LAS file into a compact 0.23 GB LAZ to our local hard drive – a file that is 10 times smaller and that will be 10 times faster to copy:

laszip -i F:\valparaiso_DSM_SURE_100.las ^
       -rescale 0.01 0.01 0.01 ^
       -o valparaiso_DSM_SURE_100.laz ^

Then we tiled the 100 million points into 250 meter by 250 meter tiles with 25 meter buffer using lastile. We use the new option ‘-flag_as_withheld’ to mark all buffer points with the withheld flag so they can be easily removed on-the-fly via the ‘-drop_withheld’ command-line filter (also see the README file).

lastile -i valparaiso_DSM_SURE_100.laz ^
        -tile_size 250 -buffer 25 ^
        -flag_as_withheld ^
        -odir valparaiso_tiles_raw -o val.laz
250 meter by 250 meter tiling with 25 meter buffer

250 meter by 250 meter tiles with 25 meter buffer

Before processing millions to billions of points we experiment with different options to find what works best on a smaller area, namely the tile “val_256750_6338500.laz” pointed to above. Using the workflow from this blog posts did not give perfect results due to the large amount of low noise. Although many low points were marked as noise (violett) by lasnoise, too many ended up classified as ground (brown) by lasground as seen here:
excessive low noise affects ground classification

excessive low noise affects ground classification

We use lascanopy – a tool very popular with forestry folks – to compute four BIL rasters where each 5m by 5m grid cell contains the 5th, 10th, 15th, and 20th percentile of the elevation values from all points falling into a cell (also see the README file):
lascanopy -i val_256750_6338500.laz ^
          -height_cutoff -1000 -step 5 ^
          -p 5 10 15 20 ^
The four resulting rasters can be visually inspected and compared with lasview:
lasview -i val_256750_6338500_*.bil -files_are_flightlines
comparing 5th and 10th elevation percentiles

comparing the 5th and the 10th elevation percentiles

By pressing the hot keys <0>, <1>, <2> and <3> to switch between the different percentiles and <t> to triangulate them into a surface, we can see that for this example the 10th percentile works well while the 5th percentile is still affected by the low noise. Hence we use the 10th percentile elevation surface and classify all points below it as noise with lasheight (also see the README file).
lasheight -i val_256750_6338500.laz ^
          -ground_points val_256750_6338500_p10.bil ^
          -classify_below -0.5 7 ^
          -odix _denoised -olaz
We visually confirm that all low points where classified as noise (violett). Note the obvious “edge artifact” along the front boundary of the tile. This is why we always recommend to use a buffer in tile-based processing.
points below 10th percentile surface marked as noise

points below 10th percentile surface marked as noise

At the end of the blog post we give the entire command sequence that first computes a 10th percentile raster with 5m resolution for the entire file with lascanopy and then marks all points of each tile below the10th percentile surface as noise with lasheight. When we classify all points into ground and non-ground points with lasground we ignore all points classified as noise. Here are the results:

DTM extracted from SGM points despite low noise

DTM extracted from dense-matching points despite low noise

corresponding DSM with all buildings and vegetaion included

corresponding DSM with all buildings and vegetaion included

Above you see the generated DTM and the corresponding DSM. So yes, LAStools can create DTMs from points that are result of dense-matching photogrammetry … even when there is a lot of low noise. There are many other ways to mix and match the modules of LAStools for more refined workflows. Sometimes declaring all points below the 10th percentile surface as noise may be too agressive. In a future blog post we will look how to combine lascanopy and lasnoise for a more adaptive approach.

:: compute 10th percentile for entire area
lascanopy -i valparaiso_DSM_SURE_100.laz ^
          -height_cutoff -1000 -step 5 ^
          -p 10 ^

:: tile input into 250 meter tiles with buffer
lastile -i valparaiso_DSM_SURE_100.laz ^
        -tile_size 250 -buffer 25 ^
        -flag_as_withheld ^
        -odir valparaiso_tiles_raw -o val.laz

:: mark points below as noise
lasheight -i valparaiso_tiles_raw/*.laz ^
          -ground_points valparaiso_DSM_SURE_100_p10.bil ^
          -classify_below -0.5 7 ^
          -odir valparaiso_tiles_denoised -olaz ^
          -cores 4

:: ground classify while ignoring noise points
 lasground -i valparaiso_tiles_denoised\*.laz ^
          -ignore_class 7 ^
          -town -bulge 0.5 ^
          -odir valparaiso_tiles_ground -olaz ^
          -cores 4 

:: create 50 cm DTM rasters in BIL format
las2dem -i valparaiso_tiles_ground\*.laz ^
        -keep_class 2 ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir valparaiso_tiles_dtm -obil ^
        -cores 4 

:: average 50 cm DTM values into single 1m DTM 
lasgrid -i valparaiso_tiles_dtm\*.bil -merged ^
        -step 1.0 -average ^
        -o valparaiso_dtm.bil

:: create hillshade adding in UTM 19 southern
blast2dem -i valparaiso_dtm.bil ^
          -hillshade -utm 19M ^
          -o valparaiso_dtm_hill.png

:: create DSM hillshade with same three steps
las2dem -i valparaiso_tiles_raw\*.laz ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir valparaiso_tiles_dsm -obil ^
        -cores 4
lasgrid -i valparaiso_tiles_dsm\*.bil -merged ^
        -step 1.0 -average ^
        -o valparaiso_dsm.bil
blast2dem -i valparaiso_dsm.bil ^
          -hillshade -utm 19M ^
          -o valparaiso_dsm_hill.png

Creating DTMs from dense-matched points of UAV imagery from SenseFly’s eBee

Tim Sutton and his team at Kartoza work on flood modelling and risk assessment using Inasafe. They have been trying to generate a DTM from point cloud data derived via dense-matching from UAV imagery collected by an eBee of SenseFly in the “unplanned developments” or “slums” North West of Dar es Salaam, the capital city of Tanzania. Tim’s team was stuck after “other software” produced this result:

results for ground points classification with other software

poor ground classification of “Tandale” with “other software”

Tim reached out to us at rapidlasso asking whether LAStools could handle this better. After all, we had published two blog articles – namely this one and that one – showing how to generate DTMs from the point clouds generated by the dense-matching photogrammetry software of Pix4D. Below the workflow we devised and the results we produced for Tim and his team.

We obtained 3 different data sets of areas called “Tandale”, “Borahatward”, and “Bugurunni”. We added one new option to our lasground software called ‘-bulge 1.0’ (see README) to improve the removal of smaller buildings and got this result for “Tandale”.

ground classification with LAStools

DTM of “Tandale” from ground points classified with LAStools

Before you point out the “facetted” look of this DTM keep in mind that “Tandale” is a densely populated poor area. A first hand account of the rough life in this area can be found here. Most dense-matching points are on corrugated roofs that become voids that need to be interpolated across in the DTM. Take a look at the corresponding DSM where all objects are still present.

original data

DSM of “Tandale” from all dense-matching points

Below we give a detailed description at the example of the “Bugurunni” data set of the workflow that was used to generate DTMs for the three data sets. At the end of this article you will see some more results.

We first use lassort to quantize, sort, and compress on 4 cores the seven spatially incoherent LAS files of the “Bugurunni” data set (totalling 4.5 GB with excessive resolution of millimeters) into LASzip-compressed files with a more reasonable resolution of centimeters and points ordered along a space-filling curve. We also add the missing projection information with ‘-utm 37M’. The resulting 7 LAZ files occupy only 0.7 GB meaning we get a compression of 9 : 1. The option ‘-odir’ specifies the output directory.

lassort -i bugurunni_densified_*.las ^
        -rescale 0.01 0.01 0.01 ^
        -utm 37M ^
        -odir bugurunni_strips -olaz ^
        -cores 4

Next we tile the sorted strips into 500 meter by 500 meter tiles with 50 meter buffer using lastile. We use the new option ‘-flag_as_withheld’ to mark all buffer points with the withheld flag so they can easily be removed on-the-fly with the ‘-drop_withheld’ command-line filter (see the README file for more on this).

lastile -i bugurunni_strips\*.laz ^
        -files_are_flightlines ^
        -tile_size 500 -buffer 50 ^
        -flag_as_withheld ^
        -o bugurunni_raw\bugu.laz
Using lasnoise on 4 cores we classify isolated points that might hinder ground-classification as noise (class 7). The parameters ‘-isolated 15’ means that all points surrounded by less than 15 other points in their 3 by 3 by 3 = 27 cells neighborhood in a 3D grid are considered isolated. The size of each grid cell is specified with ‘-step_xy 2 -step_z 1’  as 2 meter by 2 meter by 1 meter. These parameters were found experimentally (see the README file for more on this).
lasnoise -i bugurunni_raw\*.laz ^
         -step_xy 2 -step_z 1 ^
         -isolated 15 ^
         -odir bugurunni_noise -olaz ^
         -cores 4
Then we run lasground on 4 cores to classify the ground points with options ‘-metro’ and ‘-bulge 1.0’. The option ‘-metro’ is a convenient short-hand for ‘-step 50’ that will remove all objects on the terrain (e.g. large buildings) that have an extend of 50 meters or less. The option ‘-bulge 1.0’ instructs lasground to be conservative and only add points that are 1 meter or less above a smoothed version of the initial ground estimate (see the README file for more on this)..
lasground -i bugurunni_noise\*.laz ^
          -ignore_class 7 ^
          -metro -bulge 1.0 ^
          -odir bugurunni_ground -olaz ^
          -cores 4
Now we use las2dem to raster a DTM from only those points that were classified as ground. The option ‘-step 0.5’ sets the output grid resolution to 0.5 meters, ‘-kill 200’ interpolates across voids of up to 200 meters, and ‘-use_tile_bb’ rasters only the original 500 meter by 500 meter tile interior but not the 50 meter buffer. This assures that the resulting raster tiling aligns without artifacts across tile boundaries. The option ‘-obil’ chooses BIL as the output raster format.
las2dem -i bugurunni_ground\*.laz ^
        -keep_class 2 ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir bugurunni_dtm -obil ^
        -cores 4
As a simply form of anti-aliasing we average each four pixels of 0.5 meter resolution into one pixel of 1.0 meter resolution with lasgrid as all LAStools can read BIL files via on-the-fly conversion to points.
lasgrid -i bugurunni_dtm\*.bil -merged ^
        -step 1.0 -average ^
        -o bugurunni_dtm.bil

Finally we create a hillshade of the DTM adding back the projection that was “lost” in the BIL file generation so that blast2dem – the extremely scalable BLAST version of las2dem – can automatically produce a KML file for display in Google Earth.

blast2dem -i bugurunni_dtm.bil ^
          -hillshade -utm 37M ^
          -o bugurunni_dtm_hill.png

For comparison we also create a DSM with the same three steps but using all points.

las2dem -i bugurunni_raw\*.laz ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir bugurunni_dsm -obil ^
        -cores 4
lasgrid -i bugurunni_dsm\*.bil -merged ^
        -step 1.0 -average ^
        -o bugurunni_dsm.bil
blast2dem -i bugurunni_dsm.bil ^
          -hillshade -utm 37M ^
          -o bugurunni_dsm_hill.png
DTM of "Bugurunni" from ground points classified with LAStools

DTM of “Bugurunni” from ground points classified with LAStools

Above you see the generated DTM and below the corresponding DSM. So yes, LAStools can create DTMs from points that are result of dense-matching photogrammetry … under one assumption: there is not too much vegetation.

DSM of "Bugurunni" from all dense-matching points

DSM of “Bugurunni” from all dense-matching points

Below also the results for the “Borahatward” data. In a future blog post we will see how to deal with the excessive low noise sometimes present in dense-matching points.

DTM of "Bo"

DTM of “Borahatward” from ground points classified with LAStools

DSM of "Borahatward" from all dense-matching points

DSM of “Borahatward” from all dense-matching points

Generating a DTM from Dense-Matching Photogrammetry

Point clouds from dense-matching photogrammetry are popular. Stored in the LAS format – or its compressed LAZ twin – these “near-LiDAR” points are easily ingested into LiDAR processing software to, for example, generate a Digital Surface Model (DSM). “Can LAStools also create a Digital Terrain Model (DTM) from such data?” is a question we often get asked by private email or via the LAStools user forum. At INTERGEO Euroasia 2014 we gave a demo of this at the UltraCam booth. We repeated this demo with different data at the Pix4D booth during the Geospatial World Forum 2014. Now the folks from Visionmap have also become curious …


Above the “topography_cadastre_switzerland_densified.las” point cloud we got at the Pix4D booth. You may have already seen this colorized point cloud at some other Pix4D event. The original points were in a highly incoherent spatial order that negativly affected all subsequent processing. Therefore we first use lassort to reorder the points into a space-filling Hilbert curve. At the same time we rescale the coordinate resolution to centimeters and compress the file into the LAZ format.

lassort -i topography_cadastre_switzerland_densified.las ^
        -rescale 0.01 0.01 0.01 -olaz
Next we use lasground to classify the ground points with options ‘-city’ and ‘-ultra_fine’. The option ‘-city’ is a convenient short-hand for ‘-step 25’ that will remove all objects on the terrain that have an extend of 25 meters or less. The option ‘-ultra_fine’ instructs lasground to spend more time on finding a good initial ground estimate. The option ‘-odix _g’ add the appendix ‘_g’ to the output file name.
lasground -i topography_cadastre_switzerland_densified.laz ^
          -city -ultra_fine ^
          -odix _g -olaz
Above you see only the points that were classified as ground (classification 2) and below you see how triangulating these points interpolates the voids where buildings and vegetation used to be.
Now we use las2dem to raster a DTM from only those points that were classified as ground. The option ‘-step 0.5’ sets the output grid resolution to 0.5 meters, ‘-ocut 2’ removes the ‘_g’ from the file name, ‘-odix _dtm’ adds the appendix ‘_dtm’ to the file name, and ‘-obil’ chooses BIL as the output raster format.
las2dem -i topography_cadastre_switzerland_densified_g.laz ^
        -keep_class 2 -thin_with_grid 0.125 -extra_pass ^
        -step 0.5 ^
        -ocut 2 -odix _dtm -obil
For comparison we also use las2dem to create the corresponding DSM. The option ‘-thin_with_grid 0.125’ makes the density of the points that are triangulated into a TIN more “compatible” with the 0.5 by 0.5 grid that the TIN is rastered with by keeping only the first point that falls within each 0.125 by 0.125 area or a maximum of 16 points per 0.5 by 0.5 raster cell. The option ‘-extra_pass’ counts all surviving points in a first pass to minimize the memory footprint needed for Delaunay TIN construction.
las2dem -i topography_cadastre_switzerland_densified_g.laz ^
        -thin_with_grid 0.125 -extra_pass ^
        -step 0.5 ^
        -ocut 2 -odix _dsm -obil


Did you know that lasview can visualize BIL files via on-the-fly conversion from grids to points? Above you see the generated DTM and below the corresponding DSM. So yes, LAStools can create DTMs from points that are result of dense-matching photogrammetry … under one assumption: there is not too much vegetation.


clone wars and drone fights

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.