Generating DTM for “fluffy” Livox MID-40 LiDAR via “median ground” points

In the last article about the Livox MID-40 LiDAR scan of Samara we quality checked the data, aligned the flight lines and cleaned the remaining spurious scan lines. In this article we will process this data into the standard products. A focus will be on generating a smooth “median ground” surface from the “fluffy” scanner data. You can get the flight lines here and follow along with the processing after downloading LAStools. Below is one result of this work.

The first processing step will be to tile the strips into tiles that contain fewer points for faster and also parallel processing. One quick “flat terrain” trick first. Often there are spurious points that are far above or below the terrain. For a relatively flat area these can be easily be identified by computing a histogram elevation values with lasinfo and then eliminated with simple drop-filters on the z coordinate.

lasinfo ^
-i Samara\Livox\03_strips_cleaned\*.laz ^
-merged ^
-histo z 1 ^
-odir Samara\Livox\01_quality ^
-o strips_cleaned_merged_info.txt

The relevant excerpts of the output of the lasinfo report are shown below:

z coordinate histogram with bin size 1.000000
bin -104 has 1
bin 5 has 1

bin 11 has 273762
bin 12 has 1387999
bin 13 has 5598767
bin 14 has 36100225
bin 15 has 53521371
bin 59 has 60308
bin 60 has 26313
bin 61 has 284
bin 65 has 10
bin 66 has 31
bin 67 has 12
bin 68 has 1
bin 83 has 3
bin 84 has 4
bin 93 has 31
bin 94 has 93
bin 95 has 17


The few points below 11 meters and above 61 meters in elevation can be considered spurious. In the initial tiling step with lastile we add simple elevation filters to drop those points on-the-fly from the buffered tiles. The importance of buffers when processing LiDAR in tiles is discussed in this article. With lastile we create tiles of size 125 meters with a buffer of 20 meters, while removing the points identified as spurious with the appropriate filters. Because the input strips have their “file source ID” in the LAS header correctly set, we use ‘-apply_file_source_ID’ to set the “point source ID” of every point to this value. This preserves the information of which point comes from which flight line.

lastile ^
-i Samara\Livox\03_strips_cleaned\*.laz ^
-drop_z_below 11 -drop_z_above 61 ^
-apply_file_source_ID ^
-tile_size 125 -buffer 20 ^
-odir Samara\Livox\05_tiles_buffered -o desman.laz

This produces 49 buffered tiles that will now be processed similarly to the workflow outlined for another lower-priced system that generates similarly “fluffy” point clouds on hard surfaces, the Velodyne HLD-32E, described here and here. What do we mean with “fluffy”? We cut out a 1 meter slice across the road with the new ‘-keep_profile’ filter and las2las and inspect it with lasview.

las2las ^
-i Samara\Livox\05_tiles_buffered\desman_331750_1093000.laz ^
-keep_profile 331790 1093071 331799 1093062 1 ^
-o slice.laz

lasview ^
-i slice.laz ^
-color_by_flightline ^
-kamera 0 274.922 43.7695 0.00195313 -0.0247396 1.94022 ^
-point_size 9

In the view below we pressed hot key twice ‘]’ to exaggerate the z scale. The “fuzziness” is that thickness of the point cloud in the middle of this flat tar road. It is around 20 to 25 centimeters and is equally evident in both flight lines. What is the correct ground surface through this 20 to 25 centimeter “thick” road? We will compute a “mean ground” that roughly falls into the middle of this “fluffy” surface,

slice of 1 meter width across the tar road in front of Omael store

The next three lasthin runs mark a sub set of low candidate points for our lasground filtering. In every 25 cm by 25 cm, every 33 cm by 33 cm and every 50 cm by 50 cm area we reclassify the point closest to the 10th percentile as class 8. In the first call to lasthin we put all other points into class 1.

lasthin ^
-i Samara\Livox\05_tiles_buffered\desman_*.laz ^
-set_classification 1 ^
-step 0.25 -percentile 10 20 -classify_as 8 ^
-odir Samara\Livox\06_tiles_thinned_01 -olaz ^
-cores 4

lasthin ^
-i Samara\Livox\06_tiles_thinned_01\desman_*.laz ^
-step 0.3333 -percentile 10 20 -classify_as 8 ^
-odir Samara\Livox\06_tiles_thinned_02 -olaz ^
-cores 4

lasthin ^
-i Samara\Livox\06_tiles_thinned_02\desman_*.laz ^
-step 0.5 -percentile 10 20 -classify_as 8 ^
-odir Samara\Livox\06_tiles_thinned_03 -olaz ^
-cores 4

Below you can see the resulting points of the 10th percentile classified as class 8 in red.

Operating only on the points classified as 8 (i.e. ignoring those classified as 1) we then run a ground classification with lasground using the following command line, which creates a “low ground” classification. .

lasground ^
-i Samara\Livox\06_tiles_thinned_03\desman_*.laz ^
-ignore_class 1 ^
-town -ultra_fine ^
-ground_class 2 ^
-odir Samara\Livox\07_tiles_ground_low -olaz ^
-cores 4

Since this is an open road this classifies most of the red points as ground points.

Using lasheight we then create a “thick ground” by pulling all those points into the ground surface that are between 5 centimeter below and 17 centimeter above the “low ground”. For visualization purposes we temporarily use class 6 to capture this thickened ground.

lasheight ^
-i Samara\Livox\07_tiles_ground_low\desman_*.laz ^
-classify_between -0.05 0.17 6 ^
-odir Samara\Livox\07_tiles_ground_thick -olaz ^
-cores 4

The “thick ground” is shown below in orange.

We go back to lasthin and reclassify in every 50 cm by 50 cm area the point closest to the 50th percentile as class 8. This is what we call the “median ground”.

lasthin ^
-i Samara\Livox\07_tiles_ground_thick\desman_*.laz ^
-ignore_class 1 ^
-step 0.5 -percentile 50 -classify_as 8 ^
-odir Samara\Livox\07_tiles_ground_median -olaz ^
-cores 4

The final “median ground” points are shown in red below. These are the points we will use to eventually compute the DTM.

We complete the fully automated classification available in LAStools by running lasclassify with the following options. See the README file for what these options mean. Note that we move the “thick ground” from the temporary class 6 to the proper class 2. The “median ground” continues to be in class 8.

lasclassify ^
-i Samara\Livox\07_tiles_ground_median\desman_*.laz ^
-change_classification_from_to 6 2
-rugged 0.3 -ground_offset 1.5 ^
-odir Samara\Livox\08_tiles_classified -olaz ^
-cores 4

Before the resulting tiles are published or shared with others we should remove the temporary buffers, which is done with lastile – the same tool that created the buffers.

lastile ^
-i Samara\Livox\08_tiles_classified\desman_*.laz ^
-remove_buffer ^
-odir Samara\Livox\09_tiles_final -olaz ^
-cores 4

And then we can publish the points via a Potree 3D Webportal using laspublish.

laspublish ^
-i Samara\Livox\09_tiles_final\desman_*.laz ^
-elevation ^
-title "Samara Mangroves" ^
-odir Samara\Livox\99_portal -o SamaraMangroves.html -olaz ^

Below a screenshot of the resulting Potree 3D Web portal rendered with Potree Desktop. Inspecting the classification will reveal a number of errors that could be tweaked manually with lasview. How the point colors were generated is not described here but I used Google satellite imagery and mapped it with lascolor to the points. The elevation colors are mapped from 14 meters to 25 meters. The intensity image may help us understand why the black tar road on the left hand side that runs from the “Las Palmeras Condos” to the beach in “Cangrejal” has no samples. It seems the intensity is lower on this side which indicates that the drone may have flown higher here – too high to for the road to reflect enough photons. The yellow view of return type indicates that despite it’s multi-return capability, the Livox MID-40 LiDAR is mostly collecting single returns.

The penetration capability of the Livox MID-40 LiDAR was less good than we had hoped for. Below thick vegetation we have too few points on the ground to give us a good digital terrain model. In the visualization below you can see that below the dense vegetation there are large black areas which are completely void of points.

Now we produce the standard product DTM and DSM at a resolution of 50 cm. Because the total area is not that big we generate temporary tiles in “raster LAZ” with las2dem and merge them into a single GeoTiff with blast2dem.

las2dem ^
-i Samara\Livox\07_tiles_ground_median\desman_*.laz ^
-keep_class 8 ^
-step 0.5 -use_tile_bb ^
-odir Samara\Livox\10_tiles_dtm_50cm -olaz ^
-cores 4

blast2dem ^
-i Samara\Livox\10_tiles_dtm_50cm*.laz -merged ^
-step 0.5 -hillshade ^
-o Samara\Livox\dtm_50cm.png

blast2dem ^
-i Samara\Livox\10_tiles_dtm_50cm*.laz -merged ^
-step 0.5 ^
-o Samara\Livox\dtm_50cm.tif

lasthin ^
-i Samara\Livox\07_tiles_ground_median\desman_*.laz ^
-step 0.5 -percentile 95 ^
-odir Samara\Livox\11_tiles_highest -olaz ^
-cores 4

las2dem ^
-i Samara\Livox\11_tiles_highest\desman_*.laz ^
-step 0.5 -use_tile_bb ^
-odir Samara\Livox\12_tiles_dsm_50cm -olaz ^
-cores 4

blast2dem ^
-i Samara\Livox\12_tiles_dsm_50cm*.laz -merged ^
-step 0.5 -hillshade ^
-o Samara\Livox\dsm_50cm.png

blast2dem ^
-i Samara\Livox\12_tiles_dsm_50cm*.laz -merged ^
-step 0.5 ^
-o Samara\Livox\dsm_50cm.tif

A big “Thank You!” to Nelson Mattie from LiDAR Latinoamerica for bringing his fancy drone to Samara and to Andre Jalobeanu from Bayesmap for his help in aligning the data. You can download the flight lines here and do the above processing on your own after downloading LAStools.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

lasinfo (210128) report for 'Samara\Drone\00_raw_ready\flightline_01.laz'
reporting all LAS header entries:
file signature: 'LASF'
file source ID: 1
global_encoding: 0
project ID GUID data 1-4: 00000000-0000-0000-0000-000000000000
version major.minor: 1.2
system identifier: 'LAStools (c) by rapidlasso GmbH'
generating software: 'las2las (version 210315)'
file creation day/year: 224/2019
header size: 227
offset to point data: 527
number var. length records: 2
point data format: 1
point data record length: 29
number of point records: 6576555
number of points by return: 6576555 0 0 0 0
scale factor x y z: 0.001 0.001 0.001
offset x y z: 661826 1092664 14
min x y z: 660986.622 1092595.013 11.816
max x y z: 661858.813 1092770.575 95.403
variable length header record 1 of 2:
reserved 0
user ID 'ScanLook'
record ID 25
length after header 0
description 'ScanLook Point Cloud'
variable length header record 2 of 2:
reserved 0
user ID 'LASF_Spec'
record ID 4
length after header 192
description 'by LAStools of rapidlasso GmbH'
Extra Byte Descriptions
data type: 1 (unsigned char), name "laser beam ID", description: "which beam ranged this return", scale: 1 (not set), offset: 0 (not set)

LASzip compression (version 3.4r3 c2 50000): POINT10 2 GPSTIME11 2 BYTE 2
reporting minimum and maximum for all LAS point record entries …
X -839378 32813
Y -68987 106575
Z -2184 81403
intensity 4 255
return_number 1 1
number_of_returns 1 1
edge_of_flight_line 0 0
scan_direction_flag 0 0
classification 1 1
scan_angle_rank -12 90
user_data 0 0
point_source_ID 0 0
gps_time 537764.192416 537854.547507
attribute0 0 31 ('laser beam ID')
number of first returns: 6576555
number of intermediate returns: 0
number of last returns: 6576555
number of single returns: 6576555
overview over number of returns of given pulse: 6576555 0 0 0 0 0 0
histogram of classification of points:
6576555 unclassified (1)

Strip Aligning of Drone LiDAR flown with Livox MID-40 over destroyed Mangrove

September 11th 2020 seemed like a fitting day to hunt down – with a powerful drone – those who destroy our common good. The latest DJI M300 RTK drone came to visit me in Samara, Guanacaste, Costa Rica and it was carrying the gAirHawk GS-MID40 UAV laser scanning system by Geosun featuring the light-weight Livox Mid 40 LiDAR. The drone is owned and operated by my friends at LiDAR Latinoamerica.

We flew a two-sortie mission covering a destroyed mangrove lagoon that was illegally poisoned, cut-down and filled in with the intention to construct a fancy resort in its place some 25 years ago. For future environmental work I wanted to get a high-resolution baseline scan with detailed topography of the meadow and what now-a-days remains of the mangroves that are part of the adjacent “Rio Lagarto” estuary. Recently the area was illegally treated with herbicides to eliminate the native herbs and the wild flowers and improve grazing conditions for cattle. Detailed topography can show how the heavy rains have washed these illegal substances into the ocean and further prove that the application of agro-chemicals in this meadow causes harm to marine life.

Here you can see a sequence of video about the LiDAR system, the preparations and the survey flights. Shortly after the flight I obtained the LiDAR from Nelson Mattie, the CEO of LiDAR Latinoamerica and ran the usual quality checks with LAStools.

lasinfo ^
-i Samara\Livox\00_raw_laz\*.laz ^
-histo intensity 16 ^
-histo gps_time 10 ^
-histo z 5 ^
-odir Samara\Livox\01_quality -odix _info -otxt ^
-cores 3

lasgrid ^
-i Samara\Livox\00_raw_laz\*.laz ^
-utm 16north ^
-merged ^
-keep_last ^
-step 0.5 ^
-density ^
-false -set_min_max 100 1000 ^
-odir Samara\Livox\01_quality ^
-o density_050cm_100_1000.png

For the density image, lasgrid counts how many last return from all flight lines fall into each 50 cm by 50 cm area, computes the desnity per square meter and maps this number to a color ramp that goes from blue via cyan, yellow and orange to red. The overall density of our survey is in the hundred of laser pulses per square meters with great variations where flight line overlap and at the survey boundary. The start and landing area as well as the place where the first flight ended and the second flight started are the two red spots of maximum density that can easily be picked out.

blue: 100 or fewer laser pulses per square meters, red: 1000 or more laser pulses per square meter

lasoverlap ^
-i Samara\Livox\00_raw_laz\*.laz ^
-utm 16north ^
-merged -faf ^
-step 0.5 ^
-min_diff 0.10 -max_diff 0.25 ^
-elevation -lowest ^
-odir Samara\Livox\01_quality ^
-o overlap_050cm_10cm_25cm.png

For the overlap image lasoverlap counts how many different flight lines overlap each 50 cm by 50 cm area and maps the counter to a color: 1 = blue, 2 = cyan, 3 = yellow, 4 = orange, and 5 of more = red. Here the result suggests that the 27 delivered LAS files do not actually correspond to the logical flight lines but that the files are chopped up in some other way. We will have Andre Jalobeanu from Bayesmap repair this for us later.

number of flight lines covering each area: blue = 1, cyan = 2, yellow – 3, orange = 4, red = 5 or more

For the difference image, lasoverlap finds the maximal vertical difference between the lowest points from all flight lines that overlap for each 50 cm by 50 cm area and maps it to a color. If this difference is less than 10 cm, the area is colored white. Differences of 25 cm or more are colored either red or blue. All open areas such as roads, meadows and rooftops should be white here we definitely have way to much red and blue in the open areas.

vertical differences below 10 cm are white but red or blue if above 25 cm

There is way too much red and blue in areas that are wide open or on roof tops. We inspect this in further detail by taking a closer look at some of these red and blue areas. For this we first spatially index the strips with lasindex so that area-of-interest queries are accelerated, then load the strips into the GUI of lasview and add the difference image into the background via the overlay option.

lasindex ^
-i Samara\Livox
\00_raw_laz\*.laz ^
-tile_size 10 -maximum -100 ^
-cores 3

lasview ^
-i Samara\Livox
\00_raw_laz\*.laz ^

using the difference image as an overlay to inspect troublesome areas

This way is easy to lasview or clip out (with las2las) those areas that look especially troublesome. We do this here for the large condominium “Las Palmeras” whose roofline and pool provide perfect features to illustrate the misalignment. As you can see in the image sequence below, there is a horizontal shift of up to 1 meter that can be nicely visualized with a cross section drawn perpendicular across the gable of the roof and – due to the inability to get returns from water – in the area without points where the pool is.

The misalignments between flight lines are too big for the data to be useful as is, so we do what we always do when we have this problem: We write an email to Andre Jalobeanu from Bayesmap and ask for help.

After receiving the LAZ files and the trajectory file Andre repaired the misalignment in two steps. The first call to his software stripalign in mode ‘-cut’ recovered a proper set of flight lines and removed most of the LiDAR points from the moments when the drone was turning. The second call to his software stripalign in mode ‘-align’ computed the amount of misalignment in this set of flight lines and produced a new set of flight lines with these errors corrected as much as possible. The results are fabulous.

lasmerge ^
-i Samara_MID40\*.laz ^
-o samaramid40.laz

stripalign ^
-uav -cut ^
-i samaramid40.laz ^
-po Samara_MID40\*.txt -po_parse ntxyzwpk ^
-G2 -cut_dist 50 ^
-O Samara_MID40\cut

stripalign ^
-uav -align ^
-i Samara_MID40\cut\*.laz ^
-po Samara_MID40\*.txt -po_parse ntxyzwpk ^
-A -G2 -full -smap -rmap -sub 2 ^
-O Samara_MID40\corr

As you can see above, the improvements are incredible. The data seems now sufficiently aligned to be useful for being processed into a number of products. One last thing to do is the removal of spurious scan lines that seem to stem from an unusual movement of the drone, like the beginning or the end of a turn.

We use lasview with option ‘-load_gps_time’ to determine the GPS time stamps of these spurious scan lines and remove them manually using las2las with option ‘-drop_gps_time_between t1 t2’ or similar. As the points are ordered in acquisition order, we can simply replay the flight by pressing ‘p’ and step forward and backward with ‘s’ and ‘S’.

Using lasview with hot keys ‘i’, ‘p’, ‘s’ and ‘S’ we find the GPS time of points from the last reasonable scan line.

Once we determined a suitable set of GPS times to remove from a flight lines we first verify our findings once more visually using lasview before actually creating the final cut with las2las.

lasview ^
-i Samara\Livox
\02_strips_aligned\samaramid40_c_13_i_13.laz ^
-drop_gps_time_below 283887060 ^
-drop_gps_time_above 283887123 ^
-filtered_transform ^
-set_classification 8 ^

visualizing which points we keep by mapping them on-the-fly to classification 8 with a filtered transform

las2las ^
-i Samara\Livox
\02_strips_aligned\samaramid40_c_13_i_13.laz ^
-drop_gps_time_below 283887060 ^
-drop_gps_time_above 283887123 ^
-odix _cut -olaz

After spending several hours of manually removing these spurious scan lines as well as deciding to remove a few short scan lines in areas of exzessive overlap we have a sufficiently aligned and cleaned data set to start the actual post-processing.

A big “Thank You!” to Andre Jalobeanu from Bayesmap for his help in aligning the data and to Nelson Mattie from LiDAR Latinoamerica for bringing his fancy drone to Samara. You can download the data here.

final density after removing turns, spurious scan lines and redundant scan lines

Swiss add liberal “Open LiDAR” and break with conservative stereotypes like Bank Secrecy, Yodeling and Punctuality

My new favorite Swiss Miss, Viola Amherd, said today “Geodaten gehören heute zur Infrastruktur wie die Straßen und die Eisenbahn”, which says that “today, geodata are part of the infrastructure like roads and railways”. The federal councilor of Switzerland announced that programmers and planners, whether private or professional, can now download the data free of charge and use it for their projects. There are almost no limits to innovation for information projects.

Hello Germany? Hello Bavaria? Hello Austria? Where is your Open-Data-Autobahn? Whatever happened to unlimited speeds and endless Fahrvergnügen … (-;

This means in Switzerland large amounts of LiDAR are now available for download from this portal here. There is amazing data there, including high-resolution ortho-photos and land cover data. However. we went straight for the LAS files and got ourselves a few tiles near the Sazmartinshorn in the example below.

The tiles can be selected in a number of ways via an interactive map and come as individually zipped LAS files. We downloaded the nine tiles indicated above (2.16 GB), unzipped the bulky LAS files (4.78 GB) and compressed them to the compact LAZ format (638 MB) with laszip. Using LAZ instead of zipped LAS would lower storage size and transmission bandwidth by a factor of 3.5. Something the stereotypically frugal people of Switzerland may want to consider … (-;

Then we process the data with a few typical command lines with the result shown below. The first uses blast2dem to create a hill shaded 1 meter DTM from the points classified as ground.

blast2dem ^
-i swisssurface3d_laz\*.laz -merged ^
-keep_class 2 -thin_with_grid 0.5 ^
-step 1.0 -hillshade ^
-o swiss_dtm_1m_hillshade.jpg

hillshade of 1 meter DTM computed with BLAST

We use lasgrid to visualize the varying last return density per 2 meter by 2 meter area across the surveyed area with a false coloring that maps 10 or fewer pulses per square meter to blue and 20 or more pulses per square meters to red.

lasgrid ^
-i swisssurface3d_laz\*.laz -merged ^
-keep_last ^
-step 2.0 ^
-density ^
-false -set_min_max 10 20 ^
-o swiss_dendity_2m_10_20.jpg

last return density per 2 meter area. blue = 10 or less, red = 20 or more

A lasinfo report reveals that the scanner used was a RIEGL and the returning pulse width was quantized in tenths of a nanosecond into the user data field. We use lasgrid to visualize the range of the pulse width between 4.0 and 6.0 nanoseconds with a false coloring. Make sure to drop the points for which no pulse width was recorded (i.e. user data is zero) to avoid artifacts in the visualization.

lasgrid ^
-i swisssurface3d_laz\*.laz -merged ^
-drop_user_data 0 -keep_last ^
-step 2.0 ^
-user_data -lowest ^
-false -set_min_max 40 60 ^
-o swiss_pulsewidth_40_60.jpg

shortest last return pulse width per 2 meter area. blue = 4.0 ns or less, red = 6.0 ns or more

Finally we created a portal with laspublish to visualize the point cloud data interactively with Potree. The four screenshots below highlight only a few of the abilities for visualizing and measuring the point cloud.

laspublish ^
-i swisssurface3d_laz\*.laz ^
-elevation ^
-odir swisssurface3d_portal ^
-title Sazmartinshorn ^
-o sazmartinshorn.html ^
-olaz -overwrite

Colored by elevation with a distance and two height measurements.
The intensity coloring reveals some scanner artifact drawn across the mountain flank.
No surprise that the return type here is predominantly yellow single returns.
Mountains scream for a coloring by elevation, here mapped from 1700 to 2600 meters.

The open data license can be found here and we are hereby naming the source.

LASmoons: Leonidas Alagialoglou

Leonidas Alagialoglou (recipient of three LASmoons)
Multimedia Understanding Group, Aristotle University of Thessaloniki
Thessaloniki, GREECE

Canopy height is a fundamental geometric tree parameter in supporting sustainable forest management. Apart from the standard height measurement method using LiDAR instruments, other airborne measurement techniques, such as very high-resolution passive airborne imaging, have also shown to provide accurate estimations. However, both methods suffer from high cost and cannot be regularly repeated.

Preliminary results of predicted CHE based on multi-temporal satellite images against ground-truth LiDAR measurements. The 3rd column depicts pixel-wise absolute error of prediction. Last column depicts pixel-wise uncertainty estimation of the prediction (in means of 3 standard deviations).

In our study, we attempt to substitute airborne measurements with widely available satellite imagery. In addition to spatial and spectral correlations of a single-shot image, we seek to exploit temporal correlations of sequential lower resolution imagery. For this we use a convolutional variant of a recurrent neural network based model for estimating canopy height, based on a temporal sequence of Sentinel-2 images. Our model’s performance using sequential space borne imagery is shown to outperform the compared state-of-the-art methods based on costly airborne single-shot images as well as satellite images.

Digital Terrain Model of a part of the study area

The experimental study area of approximately 940 squared km is includes two national parks, Bavarian Forest National Park and Šumava National Park, which are located at the border between Germany and Czech Republic. LiDAR measurements of the area from 2017 and 2019 will be used as ground truth height measurements that have been provided by the national park’s authorities. Temporal sequences of Sentinel-2 imagery will be acquired from the Copernicus hub for canopy height estimation.

LAStools processing:
Accurate conversion of LAS files into DEM and DSM in order to acquire ground truth canopy height model.
1) Remove noise [lasthin, lasnoise]
2) Classify points into ground and non-ground [lasground, lasground_new]
3) Create DTMs and DSMs [lasthin, las2dem]

LASmoons: Zak Kus

Zak Kus (recipient of three LASmoons)
Topology Enthusiast
San Francisco, USA

While LiDAR data enables a lot of research and innovation in a lot of fields, it can also be used to create unique and visceral art. Using the high resolution data available, a 3D printer, and a long tool chain, we can create a physical, 3D topological map of the San Francisco bay area that shows off both the city’s hilly geology, and its unique skyline.


Test print of San Francisco’s Golden Gate Park.


Test print of San Francisco’s Golden Gate Park.

The ultimate goal of this project is to create an accurate, unique physical map of San Francisco, and the surrounding areas, which will be given to a loved one as a birthday gift. Using the data from the 2010 ARRA-CA GoldenGate survey, we can filter and process the raw lidar data into a DEM format using LAStools, which can be converted using a python script into a “water tight” 3D printable STL file.

While the data works fairly well out of the box, it does require a lot of manual editing, to remove noise spikes, and to delineate the coast line from the water in low lng areas. Interestingly, while many sophisticated tools exist to edit STLs that could in theory be used to clean up and prepare the files at the STL stage, few are capable of even opening files with so much detailed data. Using LAStools to manually classify, and remove unwanted data is the only way to achieve the desired level of detail in the final piece.

LiDAR data provided through USGS OpenTopography, using the ARRA-CA GoldenGate 2010 survey
+ Average point density of 3.33 pts/m^2 (though denser around SF)
+ Covers 2638 km^2 in total (only a ~100 km^2 subset is used)

LAStools processing:
Remove noise [lasnoise]
2) Manually clean up shorelines and problematic structures [lasview, laslayers]
3) Combine multiple tiles (to fit 3d printer) [lasmerge]
4) Create DEMs (asc format) for external tool to process [las2dem]

LASmoons: Martin Romain

Martin Romain (recipient of three LASmoons)
Marshall Islands Conservation Society
Majuro, Republic of the MARSHALL ISLANDS

As a low-lying coastal nation, the Republic of the Marshall Islands (RMI) is at the forefront of exposure to climate change impacts. RMI has a strong dependence on natural resources and biodiversity not only for food and income but also for culture and livelihood. However, these resources are threatened by rising sea levels and associated coastal hazards (king tides, storm surges, wave run-up, saltwater intrusion, erosion). This project aims at addressing the lack of technical capacity and available data to implement effective risk reduction and adaptation measures, with a particular focus on inundation mapping and local evacuation planning in population centers.


Typical low-lying coastal area of the Republic of the Marshall

This project intends to use LAStools to generate a DEM of the inhabited sections of 3 remote atolls (Aur, Ebon, Likiep) and 1 island (Mejit). The resulting DEM will be used to produce an inundation exposure model (and map) under variable sea level rise projections for each site. The ultimate goal is to integrate the results into each site’s disaster risk reduction strategy (long-term outcome) and present it through community consultations in schools, community centers, and council houses.

Aerial imagery of 11.5 square kilometers of land (6.3% of total national landmass) using DJI Matrice 200 V2 & DJI Zenmuse X5S with a minimum overlap of 75/75 and maximum altitude of 120m.

LAStools processing:
1) tile large point cloud into tiles with buffer [lastile]
2) remove noise points [lasthin, lasnoise]
3) classify points into ground and non-ground [lasground]
4) create Digital Terrain Models and Digital Surface Models [lasthin, las2dem]

Potential LAStools pipelines:
Removing Excessive Low Noise from Dense-Matching Point Clouds
Digital Pothole Removal: Clean Road Surface from Noisy Pix4D Point Cloud
Creating DTMs from dense-matched points of UAV imagery from SenseFly’s eBee

Coco Loco for LiDAR: Open Data and 3D GIS Carnival Cruise into Sao Paulo

The municipality of São Paulo just dropped a MASSIVE amount of open LiDAR upon the geospatial community – a whopping 33 billion points. São Paulo city now has their own Entwine Point Tiles (EPT) service available for everyone interested in downloading, visualizing, or otherwise accessing this amazing data set. The EPT point cloud service can be accessed here and an interactive 3D portal that allows you to explore this massive data with nothing but a Web browser can be accessed here. It simple serves these EPT tiles via Potree.


What you see above is the beautiful Matarazzo Building (Edifício Matarazzo in Portuguese) where the City Hall of Sao Paulo is located. On May 27 2019 we had a one day workshop in the similarly beautiful Martinelli building (Edifício Martinelli in Portuguese) to brainstorm how to process the data and I seem to remember being rather VOCAL about repeating the success of Open LiDAR in Guadalajara and also providing easy access to the data for Sao Paolo. KUDOS to those who made it happen.

This data was originally collected for the Urban Agriculture project “São Paulo: Growing Farmers’ Income, Shrinking Urban Sprawl” that was supported after winning the 2016 MAYORS CHALLENGE by Bloomburg Philanthropies. At this point I would like to thank former New York mayor Michael Rubens Bloomburg  (“Hey Mike, do they also constantly miss-spell your family name? So annoying!”) for ending his 2020 bid for the Democratic presidential nomination.

This project stands on the shoulders of many many giants such as Howard Butler and his team who worked for years to create Entwine, Markus Schuetz who came out of nowhere and suddenly dropped Potree upon our community, and the many contributors in person at OSGeo hackfests and codesprints or remotely via patches and pull requests on the software infrastructure behind this effort. Don’t forget to support these folks!

Above other visualizations of the the Matarazzo Building where the City Hall of Sao Paulo is located. Below an overview of the city after zooming out a bit.

Surprise Release of Airborne LiDAR in Germany’s “Closed Data State” Bavaria

You have guessed correctly. This is mostly fake news as our “Freistaat” (read “Free State“) of Bavaria continues to tightly guard all of its tax-payer funded geospatial basis data for no good reason. Our other “Free State” – that of Thuringia – has become one poster child of open data in Germany with North Rhine-Westphalia being the original one. But there is “some” open LiDAR in Bavaria now.

The authors of a recent paper on change detection in urban areas have published two interesting airborne LiDAR data sets from 2008 and 2009 for the town of Abenberg in [Hebel, Arens, Stilla, 2013]. What is interesting about these data sets is that they (a) were flown with a forward looking laser scanner with flight trajectories from 4 different directions (as illustrated in the image below) and (b) were surveyed again the same way in the following year for temporal change detection.


Our lasview is able to visualize how the four different flight lines scan this house from four different directions, once we have reconstructed a properly populated LAZ file with flight line information, return numbers, and unique GPS time stamps.

The data is provided for download here as zipped ASCII files that have one line per returns containing 11 comma-separated values. Below is a sample of the first 10 lines of the 2008 data:

1, 290.243, 28.663, -11.787, 0.060, -0.052, 0.997, 517.3170, -58.6934, 313.0817, 52
1, 290.208, 28.203, -11.825, 0.062, -0.056, 0.996, 517.3167, -58.6934, 313.0817, 49
1, 290.182, 27.739, -11.852, 0.063, -0.055, 0.997, 517.3164, -58.6935, 313.0817, 53
1, 290.165, 27.272, -11.866, 0.061, -0.058, 0.996, 517.3161, -58.6935, 313.0817, 53
1, 290.163, 26.800, -11.858, 0.061, -0.053, 0.997, 517.3157, -58.6935, 313.0817, 68
1, 290.152, 26.334, -11.864, 0.059, -0.054, 0.997, 517.3154, -58.6936, 313.0817, 57
1, 290.092, 25.882, -11.938, 0.050, -0.057, 0.997, 517.3151, -58.6936, 313.0817, 57
1, 290.103, 25.406, -11.911, 0.043, -0.058, 0.997, 517.3147, -58.6937, 313.0817, 63
1, 290.067, 24.947, -11.952, 0.043, -0.061, 0.997, 517.3144, -58.6937, 313.0817, 63
1, 290.034, 24.488, -11.989, 0.044, -0.063, 0.997, 517.3141, -58.6937, 313.0817, 56

The first number is either a classification into ground, vegetation, or other surface, or represents an identifier for a planar shape that the return is part of. The next three numbers in red are the x, y, and z coordinate of the LiDAR point in a local coordinate system. The next three numbers in green are the x, y, and z coordinates of an estimated surface normal. The next three numbers in blue are the x, y, and z coordinates of the sensor position. The last number is the intensity of the LiDAR return.

This textual representation makes it difficult to efficiently load the data into most LiDAR processing software. Also several attributes such as return number, number of returns, flight line ID, and GPS time stamps are missing.

We use this as an opportunity for a little exercise in converting ASCII to LAZ while preserving any “additional attributes” using the “extra bytes” functionality available since the LAS 1.4 specification. This is a timely experiment as the LAS Working Group of the ASPRS is currently contemplating how to standardize some useful “additional attributes”. Here you can download the resulting files:

In order to replicate these steps, please get your hands of the very latest version of LAStools. First we use txt2las to convert the TXT file to LAZ format as follows:

txt2las -i abenberg_data_2008.txt ^
        -set_point_type 1 ^
        -parse 0xyz123456i ^
        -set_scale 0.001 0.001 0.001 ^
        -add_attribute 3 "planar shape ID" "preliminary classification" ^
        -add_attribute 4 "normal x coord" "local normal direction estimate" 0.001 ^
        -add_attribute 4 "normal y coord" "local normal direction estimate" 0.001 ^
        -add_attribute 4 "normal z coord" "local normal direction estimate" 0.001 ^
        -add_attribute 6 "sensor x coord" "sensor position" 0.0001 ^
        -add_attribute 6 "sensor y coord" "sensor position" 0.0001 ^
        -add_attribute 6 "sensor z coord" "sensor position" 0.0001 ^
        -odix _temp1 -olaz

We use a back and forth of lasview and las2las with option ‘-subseq 12345 67890’ to interactively find the exact index of the return where each flightline is ending and the next one is starting. The command below allows you to visualize the the trajectories.

lasview -i abenberg_data_2008_temp1.laz ^
        -copy_attribute_into_x 4 ^
        -copy_attribute_into_y 5 ^
        -copy_attribute_into_z 6 ^
        -points 10000000 ^
        -point_size 5


By hovering with the mouse over a point where the trajectory end and pressing ‘i’ like info we can query the coordinates and attributes of a point. The console output also lists the index of the point in the file. We use this index as the start index for the manual search of the exact index where the flight lines really ends by piping the coordinates as text to stdout and looking for a “jump” in coordinates that indicates the start of a new flightline as shown below.

las2las -i abenberg_data_2008_temp1.laz ^
        -subseq 1213485 1213505 ^
        -oparse xyz -otxt -stdout

-279.846 -98.442 -11.973
-279.984 -98.900 -12.123
-280.150 -99.357 -12.311
-280.195 -99.825 -12.332
-280.229 -100.294 -12.340
-280.275 -100.763 -12.363
-280.302 -101.233 -12.361
-280.344 -101.700 -12.379
-280.371 -102.171 -12.376
-280.156 -102.661 -12.042
110.259 304.077 3.873
110.646 304.118 3.925
111.036 304.154 3.970
111.424 304.167 3.982
111.811 304.211 4.037
112.201 304.252 4.088
112.588 304.278 4.118
112.976 304.315 4.164
113.364 304.336 4.186
113.752 304.344 4.192

Then we use the found index to seperate the first flight line form the rest with with two more runs of las2las:

las2las -i abenberg_data_2008_temp1.laz ^
         -subseq 0 1213495 ^
         -odix _strip1 -olaz

las2las -i abenberg_data_2008_temp1.laz ^
        -subseq 1213495 100000000 ^
        -odix _rest234 -olaz

We then repeat this procedure for the rest until we have four individual strips. Next we use las2las to change the point type from 0 to 1 to have a GPS time attribute that we will then populate with “fake” but useful GPS time stamps:

las2las -i abenberg_data_2008_temp1_strip*.laz ^
        -set_point_type 1 ^
        -odix _pt1 -olaz

We noticed in the original text file that subsequent groups of points often have the exact same value for the three numbers in blue that are the x, y, and z coordinates of the sensor position. This suggests that those points are multiple returns from the same laser shot. We wrote a small tool using LASlib that exploits this observed pattern to recover the “return number” and the “number of returns” attribute for each point and to set a “fake” but unique GPS time for each such group of returns. You can download the source code for this tool here. We run this tool for each strip with a different start GPS time:

lasrecover -i abenberg_data_2008_temp1_strip1_pt1.laz ^
           -gpstime_start 1000000 ^
           -odix _rec -olaz

lasrecover -i abenberg_data_2008_temp1_strip2_pt1.laz ^
           -gpstime_start 2000000 ^
           -odix _rec -olaz

lasrecover -i abenberg_data_2008_temp1_strip3_pt1.laz ^
           -gpstime_start 3000000 ^
           -odix _rec -olaz

lasrecover -i abenberg_data_2008_temp1_strip4_pt1.laz ^
           -gpstime_start 4000000 ^
           -odix _rec -olaz

Finally we merged the four strips back into one file using lasmerge:

lasmerge -i abenberg_data_2008_temp1_strip1_pt1_rec.laz ^
         -i abenberg_data_2008_temp1_strip2_pt1_rec.laz ^
         -i abenberg_data_2008_temp1_strip3_pt1_rec.laz ^
         -i abenberg_data_2008_temp1_strip4_pt1_rec.laz ^
         -faf ^
         -o abenberg_data_2008_temp2.laz

With lasview we are now able to visualize not only the multiple returns per shot but also the different angles from which the laser scanner has observed the scene.

Finally we use las2las and a “filtered transform” with the custom classfication code provided in the first “additional attribute” to populate the official LAS classification codes as ratified by the ASPRS. First we turn code 1 (“Ground level”) into ground points.

las2las -i abenberg_data_2008_temp2.laz ^
        -keep_attribute_between 0 1 1 ^
        -filtered_transform ^
        -set_classification 2 ^
        -o abenberg_data_2008_temp3.laz


Then we turn code 5 (“Vegetation”) into high vegetation points, code 6 (“Other Surface”) into keypoints, and code 9 or higher (“Planar shape #”) into building points.

las2las -i abenberg_data_2008_temp3.laz ^
        -keep_attribute_between 0 5 5 ^
        -filtered_transform ^
        -set_classification 5 ^
        -o abenberg_data_2008_temp4.laz

las2las -i abenberg_data_2008_temp4.laz ^
        -keep_attribute_between 0 6 6 ^
        -filtered_transform ^
        -set_classification 8 ^
        -o abenberg_data_2008_temp5.laz

las2las -i abenberg_data_2008_temp5.laz ^
         -keep_attribute_above 0 8 ^
         -filtered_transform ^
         -set_classification 6 ^
         -o abenberg_data_2008.laz

Then we are done. Here you can download the resulting files:


Hebel, M., Arens, M., Stilla, U., 2013. Change detection in urban areas by object-based analysis and on-the-fly comparison of multi-view ALS data. ISPRS Journal of Photogrammetry and Remote Sensing 86, pp. 52-64. [DOI: 10.1016/j.isprsjprs.2013.09.005] [PDF]

Philippines use Taal Vulcano Eruption as Opportunity to become Very First Asian Country with Open LiDAR

UPDATE: As of January 30th also orthophotos and classified LAZ tiles are available for download.

It took just a few years of nagging, a vulcanic eruption, and then a few more weeks of nagging but now it has happened. The Philippines have become the first country in Asia to offering LiDAR as open data for free and unencumbered download. The portal created by the UP Training Center for Applied Geodesy and Photogrammetry (UP TCAGP) and their DREAM and PHIL LiDAR program already offers LiDAR-derived 1 meter DTM and DSM data flown between 2013 and 2017 as part of a national mission to aquire flood mapping data for a certain area around the Taal Vulcano. In the coming days orthophotos and the classified LiDAR point cloud will be added (at the moment the data is still undergoing another quality assurance review process).

As a quick test we went to the new online portal and downloaded the 34 DTM raster tiles that cover the Taal Vulcano Lake as seen in the screenshot below.


Downloading the area-of-interest is easy with LiPAD’s nice download portal.

The downloaded 1 meter DTM tiles are in TIF format and each cover an area of 1000 by 1000 meter. However, they are overlapping because they have a 50 meter buffer, so that each raster contains elevation samples organized in 1100 columns by 1100 rows plus “no data” values. We use two LAStools commands to remove the buffers. First we use our new demzip to turn the TIF to RasterLAZ format. Use demzip from version 200131 of LAStools (or newer) as older releases did not handle “no data” values correctly.

demzip -i Taal\DTM\*.tif ^

The conversion from TIF to RasterLAZ also reduces the total file size for the 34 files from 157 MB to 27 MB. Next we remove the buffers using a new functionality in lasgrid (make sure you have the latest LAStools version 200112 or newer).

lasgrid -i Taal\DTM\*.laz ^
        -step 1 ^
        -use_tile_size 1000 ^
        -odir Taal\DTM_unbuffered ^

Without buffers the total file size in RasterLAZ format shrinks to 22 MB. Now we have the data in a format that can either be treated as a raster or as a point cloud. Hence we can use laspublish and quickly create a visualization of the Taal Vulcano Island with Potree which we then copied onto our university Web space for you to play with.  This was he are able to instantly create an 3D visualization portal that lets anyone do various simple and also more complex measurements.

laspublish -i Taal\DTM_unbuffered ^
           -elevation ^
           -odir Taal\DTM_portal ^
           -o TaalVulcanoIsland.html ^
           -title "DTM of Taal Vulcano Island" ^
           -description "DTM of Taal Vulcano Island" ^
           -olaz -overwrite

Below we see the result visualized with the Desktop version of Potree. You can access the interactive portal we have created here with any Web browser.


Visualizing the 1 meter DTM of Taal Vulcano Island as RasterLAZ point cloud with Potree to instantly create interactive portal allowing simple measurements that give an intuition about the height and the size of the vulcanic formation that makes up Taal Vulcano Island.

We would like to acknowledge the UP Training Center for Applied Geodesy and Photogrammetry (UP TCAGP) and their DREAM and PHIL LiDAR program for providing easy and unencumbered open access to this data with a license that encourages data reuse and repurposing. Kudos for being first in Asia to make open LiDAR happen!!!