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 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 too much red and blue in areas that are open and along the roof tops. We can inspect this in further detail by inspecting some of the worst areas. For this we first spatially index the strips with lasindex so that future 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 options.

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

lasview ^
-i 00_raw_laz\*.laz ^
-gui

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.

TO BE CONTINUED

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

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

Goal:
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

Data:
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

Background:
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.

lasmoons_zak_kus_0

Test print of San Francisco’s Golden Gate Park.

lasmoons_zak_kus_1

Test print of San Francisco’s Golden Gate Park.

Goal:
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.

Data:
+
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:
1)
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

Background:
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.

DCIM100MEDIADJI_0507.JPG

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

Goal:
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.

Data:
+
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:
1)
Removing Excessive Low Noise from Dense-Matching Point Clouds
2)
Digital Pothole Removal: Clean Road Surface from Noisy Pix4D Point Cloud
3)
Creating DTMs from dense-matched points of UAV imagery from SenseFly’s eBee

LASmoons: Volga Lipwoni

Volga Lipwoni (recipient of three LASmoons)
Department of Geography, School of Earth and Environment
University of Canterbury, NEW ZEALAND

Background:
Structure from motion (SfM) photogrammetry, has emerged as an effective tool to accurately extract three-dimensional (3D) structures from a series of overlapping two-dimensional (2D) Unmanned aerial vehicles (UAVs) images. The bid to switch from the current labour-intensive, and time consuming forestry inventory practices has seen a lot of interest geared towards understanding the use of SfM photogrammetry to derive forest metrics (Iglhaut et al., 2019). There are a range of commercial, free and open source SfM photogrammetric software packages that can be used to process UAV images into 3D point clouds. Selection of the most appropriate package has become an important issue for most projects (Turner, Lucieer, & Wallace, 2013). A comparison of software performance in terms of accuracy, processing times and related costs would help foresters in deciding the best tool for the job.

lasmoons_Volga_Lipwoni

Typical point cloud derived with SfM software from UAV imagery.

Goal:
The study will generate 3D point clouds of images of a young forest trial and LAStools will be used to derive canopy height models (CHM) for computing tree heights. Tree heights from LiDAR data will serve as a baseline for accuracy assessment of heights derived from the point clouds.

Data:
+
422 UAV images processed into 3D point clouds using ten (10) different commercial and open source SfM software packages

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 Modelsand Digital Surface Models [lasthin, las2dem]
5) produce Canopy Height Models for computing tree heights [lasheight, las2dem]

References:
Iglhaut, J., Cabo, C., Puliti, S., Piermattei, L., O’Connor, J., & Rosette, J. (2019). Structure from motion photogrammetry in forestry: A review. Current Forestry Reports, 5(3), 155-168. doi:https://doi.org/10.1007/s40725-019-00094-3
Turner, D., Lucieer, A., & Wallace, L. (2013). Direct georeferencing of ultrahigh-resolution UAV imagery. EEE Transactions on Geoscience and Remote Sensing, 52(5), 2738-2745. doi:10.1109/TGRS.2013.2265295

Removing Noise from Single Photon LiDAR to Generate a Smooth DTM

A while back we had a first look at the Single Photon LiDAR from Leica’s SPL100 sensor (that eventually turned out just to be an SPL99 because one beamlet or one receiver in the 10 by 10 array was broken and did not produce any returns). Today we are taking a closer look at a strategy to remove the excessive noise in the raw Single Photon LiDAR data from a “proper” SPL100 sensor (where all of the 100 beamlets are firing) that was flown in 2017 in Navarra, Spain.

navarra_spl_teaser

Profile through original points on top of generated DTM.

The data was provided as open data by the cartography section of Navarra’s Government and is available via a simple download FTP portal. We describe the LAStools processing steps that were used to eliminate the excessive noise and to generate a smooth DTM. In the following we are using the originally released version of the data, that we obtained shortly after the portal went online that seems to be a bit more “raw” than the current files available now. One starndard quality check with lasinfo was done with:

lasinfo -i 0_raw\*.laz ^
        -cd ^
        -histo intensity 1 ^
        -histo user_data 1 ^
        -histo point_source 1 ^
        -histo gps_time 10 ^
        -odir 1_quality -odix _info -otxt

Upon inspecting the lasinfo report we suggest a few changes in how to store this Single Photon LiDAR data for more efficient hosting via an online portal. We perform these changes here before starting the actual processing. First we use the las2las call shown below to fix an error in the global encoding bits, remove an irrelevant VLR, re-scale the coordinates from millimeter to centimeters, re-offset the coordinates to nice numbers, and – what is by far the most crucial change for better compression – remap the beamlet ID stored in the ‘user data’ field as described in an earlier article.

las2las -i 0_raw\*.laz ^
        -rescale 0.01 0.01 0.01 ^
        -auto_reoffset ^
        -set_global_encoding_gps_bit 1 ^
        -remove_vlr 1 ^
        -map_user_data beamlet_ID_map.txt ^
        -odir 2_fix_rescale_reoffset_remap -olaz ^
        -cores 3

Then we use two lassort calls, one to maximize compression and one to improve spatial coherence. One lassort call rearranges the points in increasing order first based on the GPS time stamps, then breaks ties based on the user data field (that stores the beamlet ID), and finally stores the returns of every beamlet ordered by return number. We also add spatial reference information in this step. The other lassort call rearranges the points into a spatially coherent layout. It uses a Z-order sort with the granularity of 50 meter by 50 meter buckets of points. Within each bucket the point order from the prior sort is kept.

lassort -i 2_fix_rescale_reoffset_remap\*.laz ^
        -epsg 25830 ^
        -gps_time ^
        -user_data ^
        -return_number ^
        -odir 2_maximum_compression -olaz ^
        -cores 3

lassort -i 2_maximum_compression\*.laz ^
        -bucket_size 50 ^
        -odir 2_spatial_coherence -olaz ^
        -cores 3

The resulting optimized nine tiles are around 200 MB each and can be downloaded as one file here or as individual tiles here:

Now we start the usual processing workflow by tiling the data with lastile into smaller 500 meter by 500 meter tiles with a 25 meter buffer. We also set the pre-existing point classification in the data to zero as we will compute our own later.

lastile -i 2_spatial_coherence\*.laz ^
        -set_classification 0 ^
        -tile_size 500 -buffer 25 -flag_as_withheld ^
        -odir 3_buffered -o yecora.laz

We notice that a large amount of the noise has intensity values below 1000. We are still a bit puzzled where those intensity values come from and what exactly they mean in a Single Photon LiDAR system. But it works. We run las2las with a “filtered transform” to set classification of all points whose intensity value is 1000 or less to the classification code 7 (aka “noise”).

las2las -i 3_buffered\*.laz ^
        -keep_intensity_below 1000 ^
        -filtered_transform ^
        -set_classification 7 ^
        -odir 4_intensity_denoised -olaz ^
        -cores 3

We then ignore this “easy-to-identify” noise and go after the remaining one with lasnoise by ignoring classification code 7 and setting the newly identified noise to classification code 9 – not because it’s “water” (the usual meaning of class 9) but because these points are drawn with a distinct blue color when checking the result with lasview.

 lasnoise -i 4_intensity_denoised\*.laz ^
         -ignore_class 7 ^
         -step_xy 1.0 -step_z 0.2 ^
         -isolated 5 ^
         -classify_as 9 ^
         -odir 4_isolation_denoised -olaz ^
         -cores 3

Of the surviving non-noise points we then use lasthin to reclassify the point closest to the 20th elevation percentile per 50 cm by 50 cm area with classification code 8 (for all areas that have more than 5 non-noise points per 50 cm by 50 cm area. We repeat the same for every 1 meter by 1 meter area.

lasthin -i 4_isolation_denoised\*.laz ^
        -ignore_class 7 9 ^
        -step 0.5 -percentile 20 5 ^
        -classify_as 8 ^
        -odir 5_thinned_p20_050cm -olaz ^
        -cores 3

lasthin -i 5_thinned_p20_050cm\*.laz ^
        -ignore_class 7 9 ^
        -step 1.0 -percentile 20 5 ^
        -classify_as 8 ^
        -odir 5_thinned_p20_100cm -olaz ^
        -cores 3

We then perform a more agressive second noise removal step one with lasnoise using only those points with classification code 8, namely those non-noise points that were the 20th elevation percentile in either a 50 cm by 50 cm cell or a 1 meter by 1 meter cell. This can be done by ignoring classification code 0, 7, and 9. We mark those noise points as 6 so they appear orange in the point cloud with lasview.

lasnoise -i 5_thinned_p20_100cm\*.laz ^
         -ignore_class 0 7 9 ^
         -step_xy 2.0 -step_z 0.2 ^
         -isolated 1 ^
         -classify_as 6 ^
         -odir 5_thinned_p20_100cm_denoised -olaz ^
         -cores 3

The 20th elevation percentile points that survive the last noise removal are then classified into ground (2) and non-ground (1) points with lasground_new by ignoring all other points, namely those with classification codes 0, 6, 7, and 9.

lasground_new -i 5_thinned_p20_100cm_denoised\*.laz ^
              -ignore_class 0 6 7 9 ^
              -town ^
              -odir 5_tiles_ground_050cm -olaz ^
              -cores 3

These images below illustrate the steps we took. They also show that not all data was used and might give you ideas where to tweak our workflow for even better results.

Finally we raster the ground points into 1 meter Digital Terrain Model (DTM) rasters with las2dem and store the result (without buffers) to the RasterLAZ format.

las2dem -i 5_tiles_ground_050cm\*.laz ^
        -keep_class 2 ^
        -step 1.0 ^
        -use_tile_bb ^
        -odir 6_tiles_dtm_100cm -olaz ^
        -cores 3

Finally we merged all RasterLAZ tiles into one and compute the final hillshaded DTM with blast2dem.

blast2dem -i 6_tiles_dtm_100cm\*.laz -merged ^
          -step 1.0 ^
          -hillshade ^
          -o yecora_dtm_100cm.png

The hillshaded DTM that is result of the entire sequence of processing steps described above is shown below.

DTM from ground classification created with LAStools

For comparison we generate the same DTM using the originally provided classification. According to the README file the original ground points are classified with code 22 in areas of flight line overlap and as the usual code 2 elsewhere. Hence we must use both classification codes to construct the DTM. We do this analogue to the earlier processing steps with the three LAStools commands lastile, las2dem, and blast2dem below.

lastile -i 2_spatial_coherence\*.laz ^
        -tile_size 500 -buffer 25 -flag_as_withheld ^
        -odir 3_tiles_buffered_orig -o yecora.laz

las2dem -i 3_tiles_buffered_orig\*.laz ^
        -keep_class 2 22 ^
        -step 1.0 ^
        -use_tile_bb ^
        -odir 6_tiles_dtm_100cm_orig -olaz ^
        -cores 3

blast2dem -i 6_tiles_dtm_100cm_orig\*.laz -merged ^
          -step 1.0 ^
          -hillshade ^
          -o yecora_dtm_100cm_orig.png

Below the hillshaded DTM generated from the ground classification that was provided with the LiDAR when it was originally released as open data.

DTM from ground classification of originally released data.

In the meantime Andorra’s SPL data have been updated with a newer version in the open data portal. The new version of the data contains a much better ground classification that might have been improved manually as the new files now have the the string ‘cam’ instead of ‘ca’ in the file name, which probably means ‘classified automatically and manually’ instead of the original ‘classified automatically’. We decided not to switch to the new data release as it seemed less “raw” than the original release. For example there are suddenly points with GPS times and returns counts and numbers of zero in the file that seem synthetic. But we also computed the hillshaded DTM for the new release which is shown below.

DTM from ground classification of newly released data.

We thank the cartography section of Navarra’s Government for providing their LiDAR as open data. This not only allows re-purposing expensive data paid for by public taxes but also generates additional value, encourages citizen science, and provides educational opportunity and insights such as this blog article.

Another European Country Opens LiDAR: Welcome to the Party, Slovakia!

We got a little note from Vítězslav Moudrý from CULS pointing out that the Geodesy, Cartography and Cadastre Authority of the Slovak Republic has started releasing LiDAR as open data on their interactive Web portal. Congratulations, Slovakia!!! Welcome to the Open Data Party!!! We managed to download some data starting from this Web portal link and describe the process of obtaining LiDAR data from the Low Tatras mountain range in central Slovakia with pictures below.

open_data_portal_slovakia_01

(1) click the new “data export” link

open_data_portal_slovakia_02

(2) change the export selection to “Shape”

open_data_portal_slovakia_03

(3) change the file format to “LAZ”

open_data_portal_slovakia_04

(4) zoom to a colored area-of-interest

open_data_portal_slovakia_05

(5) zoom further and draw a nice polygon

open_data_portal_slovakia_06

(6) edit polygon into nice shape and realize heart is red because area is too big

open_data_portal_slovakia_07

(7) zoom further and draw polygon smaller than 2 square kilometer

open_data_portal_slovakia_08

(8) when polygon turns green, accept license, enter email address and export

open_data_portal_slovakia_09

(9) short wait and you get download link to such an archive

open_data_portal_slovakia_10

(10) license conditions: PDF auto-translated from Slovak to English

 

open_data_portal_slovakia_11

(11) LiDAR are spatially indexed flight lines clipped to area-of-interest

open_data_portal_slovakia_12_density_all_returns_20_50

(12) all return density: blue = 20 and red = 50 returns per square meter

lasgrid -i LowTatras\*.laz -merged ^
        -step 2 -point_density_16bit ^
        -false -set_min_max 20 50 ^
        -o LowTatras\density_all_returns_20_50.png

open_data_portal_slovakia_13_density_last_returns_4_40

(13) last return density: blue = 4 and red = 40 last returns per square meter

lasgrid -i LowTatras\*.laz -merged ^
        -keep_last ^
        -step 2 -point_density_16bit ^
        -false -set_min_max 4 40 ^
        -o LowTatras\density_last_returns_4_40.png

open_data_portal_slovakia_14_density_ground_returns_4_40

(14) ground return density: blue = 4 and red = 40 ground returns per square meter

lasgrid -i LowTatras\*.laz -merged ^
        -keep_classification 2 ^
        -step 2 -point_density_16bit ^
        -false -set_min_max 4 40 ^
        -o LowTatras\density_ground_returns_4_40.png

open_data_portal_slovakia_14_overlap_10cm_20cm_diff

(15) flight line difference image: white <= +/- 10 cm and red/blue >= +/- 20 cm

lasoverlap -i LowTatras\*.laz -faf ^
           -drop_classification 7 18 ^
           -min_diff 0.1 -max_diff 0.2 ^
           -o LowTatras\overlap_10cm_20cm.png

Finally we compute a DSM and a corresponding DTM using the already existing ground classification with BLAST using the command sequence shown below.

 

lasthin -i LowTatras\*.laz -merged ^
        -drop_classification 7 18 ^
        -step 0.5 -highest ^
        -o LowTatras\highest_50cm.laz

blast2dem -i LowTatras\highest_50cm.laz ^
          -hillshade ^
          -o LowTatras -o dsm_1m_hillshaded.png

blast2dem -i LowTatras\*.laz -merged ^
          -keep_classification 2 ^
          -thin_with_grid 0.5 ^
          -hillshade ^
          -o LowTatras\dtm_1m_hillshaded.png

We thank the Geodesy, Cartography and Cadastre Authority of the Slovak Republic for providing their LiDAR as open data for both commercial and non-commercial purposes and name the source of the data used above (as the license requires) as the Office of Geodesy, Cartography and Cadastre of the Slovak Republic (GCCA SR) or – in Slovak – the Úrad geodézie, kartografie a katastra Slovenskej republiky (ÚGKK SR).

Which European country goes next? Czech Republic? Poland? Hungary? Switzerland?

 

 

Completeness and Correctness of Discrete LiDAR Returns per Laser Pulse fired

Again and again we have preached about the importance of quality checking when you first get your expensive LiDAR data from the vendor or your free LiDAR data from an open data portal. The minimal quality check we usually advocate consists of lasinfo, lasvalidate, lasoverlap, and lasgrid. The information computed by these LAStools can reassure you that the data contains the right information, is specification conform, has properly aligned flight lines, and has the density distribution you expect. For deliveries or downloads in LAZ format we in addition recommend running laszip with the option ‘-check’ to find the rare file that might have gotten bit-corrupted or truncated during the transfer or the download. Today we learn about a more advanced quality check that can be done by running lassort followed by lasreturn.

For every laser shot fired there are usually between one to five discrete LiDAR returns and some full-waveform systems may even deliver up to fifteen returns. Each of these one to fifteen returns is then given the exact same GPS time stamp that corresponds to the moment in time the laser pulse was fired. By having these unique GPS time stamps we can always recover the set of returns that come from the same laser shot. This makes it possible to check completeness (are all the returns in the file) and correctness (is the returns numbering correct) for the discrete returns of each laser pulse.

optech_galaxy_issue

Showing all sets of returns in the file that do not have an unique GPS time stamp because the set has one or more duplicate returns (e.g. two first returns, two second returns, … ).

With LAStools we can do this by running lassort followed by lasreturn for any LiDAR that comes from a single beam system. For LiDAR that comes from some multi-beam system, such as the Velodyne 16, 32, 64, or 128, the Optech Pegasus, the RIEGL LMS 1560 (aka “crossfire”), or the Leica ALS70 or ALS80 we first need to seperate the files into one file per beam, which can be done with lassplit.  In the following we investigate data coming from an Optech Galaxy single-beam system. First we sort the returns by GPS time stamp using lassort (this step can be omitted if the data is already sorted in acquisition order (aka by increasing GPS time stamps)) and then we check the return numbering with lasreturn:

lassort -i L001-1-M01-S1-C1_r.laz -gps_time -odix _sorted -olaz

lasreturn -i L001-1-M01-S1-C1_r_sorted.laz -check_return_numbering
checked returns of 11809046 multi and 8585573 single return pulses. took 26.278 secs
missing: 0 duplicate: 560717 too large: 0 zero: 0
duplicate
========
200543 returns with n = 1 and r = 1 are duplicate
80548 returns with n = 2 and r = 1 are duplicate
80548 returns with n = 2 and r = 2 are duplicate
41962 returns with n = 3 and r = 1 are duplicate
41962 returns with n = 3 and r = 2 are duplicate
41962 returns with n = 3 and r = 3 are duplicate
13753 returns with n = 4 and r = 1 are duplicate
13753 returns with n = 4 and r = 2 are duplicate
13753 returns with n = 4 and r = 3 are duplicate
13753 returns with n = 4 and r = 4 are duplicate
3636 returns with n = 5 and r = 1 are duplicate
3636 returns with n = 5 and r = 2 are duplicate
3636 returns with n = 5 and r = 3 are duplicate
3636 returns with n = 5 and r = 4 are duplicate
3636 returns with n = 5 and r = 5 are duplicate
WARNING: there are 59462 GPS time stamps that have returns with different number of returns

The output we see above indicates a problem in the return numbering. A recently added new options to lasreturn that allow to reclassify those returns that seem to be part of a problematic set of returns that either contains missing returns, duplicate returns, or returns with different values for the “numbers of returns of given pulse” attribute. This allows us to visualize the issue with lasview. All returns whose are part of a problematic set is shown in the image above.

lasreturn -i L001-1-M01-S1-C1_r_sorted.laz ^
          -check_return_numbering ^
          -classify_as 8 ^
          -classify_duplicate_as 9 ^
          -classify_violation_as 7 ^
          -odix _marked -olaz

This command will mark all sets of returns (i.e. returns that have the exact same GPS time stamp) that have missing returns as 8, that have duplicate returns as 9, and that have returns which different “number of returns per pulse” attribute as 7. The data we have here has no missing returns (no returns are classified as 8) but we have duplicate (9) and violating (7) returns. We look at them closely in single scan lines to conclude.

It immediately becomes obvious that the same GPS time stamp was assigned to the returns of pair of subsequent shots. If the subsequent shots have the same number of returns per shot they are classified as duplicate (9 or blue). If the subsequent shots have different number of returns per shot they are marked as violating (7 or violett) but the reason for the issue is the same. We can look at a few of these return sets in ASCII. Here two subsequent four return shots that have the same GPS time stamp.

237881.011730 4 1 691602.736 5878246.425 141.992 6 79
237881.011730 4 2 691602.822 5878246.415 141.173 6 89
237881.011730 4 3 691603.051 5878246.389 138.993 6 44
237881.011730 4 4 691603.350 5878246.356 136.150 6 169
237881.011730 4 1 691602.793 5878246.439 142.037 6 114
237881.011730 4 2 691602.883 5878246.429 141.185 6 96
237881.011730 4 3 691603.109 5878246.404 139.033 6 50
237881.011730 4 4 691603.414 5878246.370 136.129 6 137

Here a four return shot followed by a three return shot that have the same GPS time stamp.

237881.047753 4 1 691603.387 5878244.501 140.187 6 50
237881.047753 4 2 691603.602 5878244.476 138.141 6 114
237881.047753 4 3 691603.776 5878244.456 136.490 6 60
237881.047753 4 4 691603.957 5878244.436 134.767 6 116
237881.047753 3 1 691603.676 5878244.492 138.132 6 97
237881.047753 3 2 691603.845 5878244.473 136.534 6 90
237881.047753 3 3 691604.034 5878244.452 134.739 6 99

It appears the GPS time counter in the LMS export software did not store the GPS time with sufficient resolution to always distinguish subsequent shots. The issue was confirmed by Optech and was already fixed a few months ago.

We should point out that these double-used GPS time stamps have zero impact on the geometric quality of the point cloud or the distribution of returns. The drawback is that not all returns can easily be grouped into one unique set per laser shot and that the files are not entirely specification conform. Any software that relies on accurate and unique GPS time stamps (such as flight line alignment software) may potentially struggle as well. The bug of the twice-used GPS time stamps was a discovery that is probably of such low consequence that no user of Optech Galaxy data had noticed it in the 4 years that Galaxy had been sold … until we really really scrutinized some data from one of our clients. Optech reports that the issue has been fixed now. But there are other vendors out there with even more serious GPS time and return numbering issues … to be continued.

Another German State Goes Open LiDAR: Saxony

Finally some really good news out of Saxony. 😊 After North Rhine-Westphalia and Thuringia released the first significant amounts of open geospatial data in Germany in a one-two punch in January 2017, we now have a third German state opening their entire tax-payer-funded geospatial data holdings to the tax-paying public via a simple and very easy-to-use online download portal. Welcome to the open data party, Saxony!!!

Currently available via the online portal are the LiDAR-derived raster Digital Terrain Model (DTM) at 1 meter resolution (DGM 1m) for everything flown since 2015 and and at 2 meter resolution (DGM 2m) or 20 meter resolution (DGM 20m) for the entire state. The horizontal coordinates use UTM zone 33 with ETRS89 (aka EPSG code 25833) and the vertical coordinate uses the “Deutsche Haupthöhennetz 2016” or “DHHN2016” (aka EPSG code 7837). Also available are orthophotos at 20 cm (!!!) resolution (DOP 20cm).

dgm_1000_rdax_87

Overview of current LiDAR holdings. Areas flown 2015 or later have LAS files and 1 meter rasters. Others have LiDAR as ASCII files and lower resolution rasters.

Offline – by ordering through either this online form or that online form – you can also get the 5 meter DTM and the 10 meter DTM, the raw LiDAR point clouds, LiDAR intensity rasters, hill-shaded DTM rasters, as well as the 1 meter and the 2 meter Digital Surface Model (DSM) for a small administrative fee that ranges between 25 EUR and 500 EUR depending on the effort involved.

Our immediate thought is to get a copy on the entire raw LiDAR points clouds (available as LAS 1.2 files for all  data acquired since 2015 and as ASCII text for earlier acquisitions) and find some portal willing to hosts this data online. We are already in contact with the land survey of Saxony to discuss this option and/or alternate plans.

Let’s have a look at the data. First we download four 2 km by 2 km tiles of the 1 meter DTM raster for an area surrounding the so called “Greifensteine” using the interactive map of the download portal, which are provided as simple XYZ text. Here a look at the contents of one ot these tiles:

more Greifensteine\333525612_dgm1.xyz
352000 5613999 636.26
352001 5613999 636.27
352002 5613999 636.28
352003 5613999 636.27
352004 5613999 636.24
[...]

Note that the elevation are not sampled in the center of every 1 meter by 1 meter cell but exactly on the full meter coordinate pair, which seems especially common  in German-speaking countries. Using txt2las we convert these XYZ rasters to LAZ format and add geo-referencing information for more efficient subsequent processing.

txt2las -i greifensteine\333*_dgm1.xyz ^
        -set_scale 1 1 0.01 ^
        -epsg 25833 ^
        -olaz

Below you see that going from XYZ to LAZ reduces the amount of  data from 366 MB to 10.4 MB, meaning that the data on disk becomes over 35 times smaller. The ability of LASzip to compress elevation rasters was first noted during the search for missing airliner MH370 and resulted in our new LAZ-based compressor for height grid called DEMzip.  The resulting LAZ files now also include geo-referencing information.

96,000,000 333525610_dgm1.xyz
96,000,000 333525612_dgm1.xyz
96,000,000 333545610_dgm1.xyz
96,000,000 333545612_dgm1.xyz
384,000,000 bytes

2,684,820 333525610_dgm1.laz
2,590,516 333525612_dgm1.laz
2,853,851 333545610_dgm1.laz
2,795,430 333545612_dgm1.laz
10,924,617 bytes

Using blast2dem we then create a hill-shaded version of the 1 meter DTM in order to overlay a visual representation of the DTM onto Google Earth.

blast2dem -i greifensteine\333*_dgm1.laz ^
          -merged ^
          -step 1 ^
          -hillshade ^
          -o greifensteine.png

Below the result that nicely shows how the penetrating laser of the LiDAR allows us to strip away the forest to see interesting geological features in the bare-earth terrain.

In a second exercise we use the available RGB orthophoto images to color one of the DTM tiles and explore it using lasview. For this we download the image for the top left of the four tiles that covers the area containing the “Greifensteine” from the interactive download portal for orthophotos. As the resolution of the TIF image is 20 cm and that of the DTM is only 1 meter, we first down-sample the TIF using gdalwarp of GDAL.

gdalwarp -tr 1 1 ^
         -r cubic ^
         greifensteine\dop20c_33352_5612.tif ^
         greifensteine\dop1m_33352_5612.tif

If you are not yet using GDAL today is a good day to start. It nicely complements the point cloud processing functionality of LAStools for raster inputs. Next we use lascolor to give each elevation pixel of the DTM stored in LAZ format its corresponding color from the orthophoto.

lascolor -i greifensteine\333525612_dgm1.laz ^
         -image greifensteine\dop1m_33352_5612.tif ^
         -odix _rgb -olaz

Now we can view the colored DTM in LAZ format interactively with lasview or any other LiDAR viewing software and turn on the RGB colors from the orthophoto as needed to understand the scene.

lasview -i greifensteine\333525612_dgm1_rgb.laz

We thank the “Staatsbetrieb Geobasisinformation und Vermessung Sachsen (GeoSN)” for giving us easy access to the 1 meter DTM and the 20 cm orthophoto that we have used in this article through their new open geodata portal as open data under the user-friendly license “Datenlizenz Deutschland – Namensnennung – Version 2.0.