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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

LASmoons: David Bandrowski

David Bandrowski (recipient of three LASmoons)
Yurok Tribe
Native American Indian Tribe in Northern California, USA

Background:
Wild spring-run Chinook salmon populations on the South Fork Trinity River in Northern California are near the brink of extinction. The South Fork Trinity River is the most remote and the largest un-dammed river in the State of California, federally designated as a wild and scenic river, and is a keystone watershed within the Klamath River basin supporting one of the last remaining populations of wild spring-run Chinook salmon. Ecosystem restoration is urgently needed to improve watershed health in the face of climate change, land use, and water diversions. This drastic decline of the wild salmon species motivated the Yurok Tribe and its partners to take action and implement this project as a last opportunity to save this species before extinction. Spring-run Chinook are extremely important for the Yurok people culturally, spiritually, and for a subsistence food source.

sample of the available photogrammetry data

Goal:
Due to budgetary constraints, airborne LiDAR is not available; therefore the Yurok Tribe has been using aerial drones and Structure for Motion (SfM) photogrammetry to develop DTM models that can be used in determining available salmon habitat and to develop prioritized locations for restoration. The watershed has extremely heavy vegetation, and obtaining bare-earth surfaces for hydraulic modeling is difficult without the proper tools. The goal is to use LAStools to further restoration science and create efficient workflows for DTM development.

Data:
+
 length of river mapped: 8 Kilometers
+ number of points: 150,856,819
+ horizontal datum: North American Datum 83 – California State Plane – Zone 1 (usft)
+ vertical datum: North American Vertical Datum 88

LAStools processing:
1) data quality checking [lasinfo, lasview, lasgrid]
2) classify ground and non-ground points [lasground and lasground_new]
3) remove low and high outliers [lasheight, lasnoise]
4) create DTM tiles at appropriate resolution [las2dem]
5) create a normalized point cloud [lasheight]

Digital Pothole Removal: Clean Road Surface from Noisy Pix4D Point Cloud

How to generate a clean Digital Terrain Model (DTM) from point clouds that were generated with the image matching techniques implemented in various photogrammetry software packages like those from Pix4D, AgiSoft, nframes, DroneDeploy and others has become an ever more frequent inquiry. There are many other blog posts on the topic that you can peruse as well [1,2,3,4,5,6,7,8]. In the following we go step by step through the process of removing low noise from a high-density point cloud that was generated with Pix4D software. A composite of the resulting DTM and DSM is shown below.

Final DSM and DTM created with LAStools for a photogrammetric point cloud of a road generated by Pix4D.

After downloading the data it is useful to familiarize oneself with the number of points, the density of points and their geo-location. This can be done with lasview, lasinfo, and lasgrid using the command lines shown below. There are around 19 million points in the file and their density averages around 2300 points per square meter. Because the RGB values have a 16 bit range (as evident in the lasinfo report) we need to add the option ‘-scale_rgb_down’ to the command line when producing the RGB raster with lasgrid.

lasview -i 0_photogrammetry\densified_point_cloud.laz

lasinfo -i 0_photogrammetry\densified_point_cloud.laz ^
        -cd ^
        -o 1_quality\densified_point_cloud.txt

lasgrid -i 0_photogrammetry\densified_point_cloud.laz ^
        -scale_rgb_down ^
        -step 0.10 ^
        -rgb ^
        -fill 1 ^
        -o 1_quality\densified_point_cloud.png

The first step is to use lastile and create smaller and buffered tiles for these 19 million photogrammetry points. We use a tile size of 100 meters, request a buffer of 10 meters around every tile, and flag buffer points as withheld so they can be easily be dropped later. We also make sure that all classification codes are reset to 0.

lastile -i 0_photogrammetry\points.laz ^
        -set_classification 0 ^
        -tile_size 100 -buffer 10 -flag_as_withheld ^
        -o 2_tiles_raw\seoul.laz -olaz

We start with lassort as a pre-processing step that rearranges the points into a more coherent spatial order which often accelerates subsequent processing steps.

lassort -i 2_tiles_raw\seoul_*.laz ^
        -odir 3_tiles_temp0 -olaz ^
        -cores 4

Next we use a sequence of four modules, namely lasthin, lasnoiselasground, and lasheight with fine-tuned parameters to remove the low noise points that are typical for point clouds generated from imagery by photogrammetry software. A typical example for such noise points are shown in the image below generated with this call to lasview:

lasview -i 3_tiles_temp0\seoul_210400_542900.laz ^
        -inside 210406 542894 210421 542921 ^
        -points 20000000 ^
        -kamera 0 -95 90 0 -0.3 1.6 ^
        -point_size 4

Ground surface noise (exaggerated by pressing <]> in lasview which doubles the scale in z).

As always, the idea is to construct a surface that is close to the ground but always above the noise so that it can be used to declare all points beneath it as noise. Below is a processing pipeline whose parameters work well for this data and that you can fine tune for the point density and the noise profile of your own data.

First we use lasthin to give those points the classification code 8 that are closest to the 70th percentile in elevation within every 20 cm by 20 cm cell. As statistics like percentiles are only stable for a sufficient number of points we only do this for cells that contain 25 points or more. Given that we have an average of 2300 points per square meter this should easily be the case for all relevant cells.

lasthin -i 3_tiles_temp0\seoul_*.laz ^
        -step 0.20 ^
        -percentile 70 25 ^
        -classify_as 8 ^
        -odir 3_tiles_temp1 -olaz ^
        -cores 4

The we run lasnoise only points on the points with classification code 8 and reclassify all “overly isolated” points with code 9. The check for isolation uses cells of size 20 cm by 20 cm by 5 cm and reclassifies the points in the center cell when the surrounding neighborhood of 27 cells has only 3 or fewer points in total. Changing the parameters for ‘-step_xy 0.20 -step_z 0.05 -isolated 3’ will remove isolated points more or less aggressive.

lasnoise -i 3_tiles_temp1\seoul_*.laz ^
         -ignore_class 0 ^
         -step_xy 0.20 -step_z 0.05 -isolated 3 ^
         -classify_as 9 ^
         -odir 3_tiles_temp2 -olaz ^
         -cores 4

Next we use lasground to ground-classify only the surviving points (that still have classification code 8) by ignoring those with classification codes 0 or 9. This sets their classification code to either ground (2) or non-ground (1). The temporary surface defined by the resulting ground points will be used to classify low points as noise in the next step.

lasground -i 3_tiles_temp2\seoul_*.laz ^
          -ignore_class 0 9 ^
          -town -ultra_fine -bulge 0.1 ^
          -odir 3_tiles_temp3 -olaz ^
          -cores 4

Then we use lasheight to classify all points that are 2.5 cm or more below the triangulated surface of temporary ground points as points as noise (7) and all others as unclassified (1).

lasheight -i 3_tiles_temp3\seoul_*.laz ^
          -classify_below -0.025 7 ^
          -classify_above -0.225 1 ^
          -odir 4_tiles_denoised -olaz ^
          -cores 4

The progress of each step is illustrated visually in the two image sequences shown below.

Now that all noise points are classified we start a standard processing pipeline, but always ignore the low noise points that are now classified with classification code 7.

The processing steps below create a 10 cm DTM raster. We first use lasthin to classify the lowest non-noise point per 10 cm by 10 cm cell. Considering only those lowest points we use lasground with options ‘-town’, ‘-extra_fine’, ‘-bulge 0.05’, and ‘-spike 0.05’. Using las2dem the resulting ground points are interpolated into a TIN and rasterized into a 10 cm DTM cutting out only the center 100 meter by 100 meter tile. We store the DTM raster as a gridded LAZ for maximal compression and finally merge these gridded LAZ files to create a hillshaded raster in PNG format with blast2dem.

lasthin -i 4_tiles_denoised\seoul_*.laz ^
        -ignore_class 7 ^
        -step 0.10 ^
        -lowest ^
        -classify_as 8 ^
        -odir 5_tiles_thinned_lowest -olaz ^
        -cores 4

lasground -i 5_tiles_thinned_lowest\seoul_*.laz ^
          -ignore_class 1 7 ^
          -town -extra_fine ^
          -bulge 0.05 -spike 0.05 ^
          -odir 6_tiles_ground -olaz ^
          -cores 4

las2dem -i 6_tiles_ground\seoul_*.laz ^
        -keep_class 2 ^
        -step 0.10 ^
        -use_tile_bb ^
        -odir 7_tiles_dtm -olaz ^
        -cores 4

blast2dem -i 7_tiles_dtm\seoul_*.laz -merged ^
          -hillshade ^
          -step 0.10 ^
          -o dtm_hillshaded.png

The processing steps below create a 10 cm DSM raster. We first use lasthin to classify the highest non-noise point per 10 cm by 10 cm cell. With las2dem the highest points are interpolated into a TIN and rasterized into a 10 cm DSM cutting out only the center 100 meter by 100 meter tile. Again we store the raster as gridded LAZ for maximal compression and merge these files to create a hillshaded raster in PNG format with blast2dem.

lasthin -i 4_tiles_denoised\seoul_*.laz ^
        -ignore_class 7 ^
        -step 0.10 ^
        -highest ^
        -classify_as 8 ^
        -odir 8_tiles_thinned_highest -olaz ^
        -cores 4

las2dem -i 8_tiles_thinned_highest\seoul_*.laz ^
        -keep_class 8 ^
        -step 0.10 ^
        -use_tile_bb ^
        -odir 9_tiles_dsm -olaz ^
        -cores 4

blast2dem -i 9_tiles_dsm\seoul_*.laz -merged ^
          -hillshade ^
          -step 0.10 ^
          -o dsm_hillshaded.png

The final result is below. The entire script is linked here. Simply download it, modify it as needed, and try it on this data or on your own data.

Scripting LAStools to Create a Clean DTM from Noisy Photogrammetric Point Cloud

A recent inquiry by Drone Deploy in the LAStools user forum gave us access to a nice photogrammetric point cloud for the village of Fillongley in the North Warwickshire district of England. They voted “Leave” with a whopping 66.9% according to the EU referendum results by the BBC. Before we say “Good riddance, Fillongley!” we EU-abuse this little village one last time and remove all their low noise points to create a nice Digital Terrain Model (DTM). The final result is shown below.

Side by side comparison of DTM and DSM generated with LAStools from photogrammetric point cloud by Drone Deploy.

After downloading the data it is useful to familiarize yourself with the point number, the point density and their geo-location, which can be done with lasview, lasinfo, and lasgrid using the command lines shown below. There are around 50 million points and their density averages close to 70 points per square meter.

lasview -i 0_photogrammetry\points.laz

lasinfo -i 0_photogrammetry\points.laz ^
        -cd ^
        -o 1_quality\fillongley.txt

lasgrid -i 0_photogrammetry\points.laz ^
        -step 0.50 ^
        -rgb ^
        -fill 1 ^
        -o 1_quality\fillongley.png

The first step is to use lastile and create smaller and buffered tiles for these 50 million photogrammetry points. We use a tile size of 200 meters, request a buffer of 25 meters around every tile, and flag buffer points as withheld so they can be easily be dropped later.

lastile -i 0_photogrammetry\points.laz ^
        -tile_size 200 -buffer 25 -flag_as_withheld ^
        -o 2_tiles_raw\fillongley.laz -olaz

Next we use a sequence of four modules, namely lasthin, lasnoiselasground, and lasheight with fine-tuned parameters to remove the low noise points that are typical for point clouds generated from imagery by photogrammetry software. A typical example for such noise points are shown in the image below.

lasview -i 2_tiles_raw\fillongley_596000_5815800.laz ^
        -inside 596050 5815775 596150 5815825 ^
        -kamera 0 -89 -1.75 0 0 1.5 ^
        -point_size 3

Clumps of low noise points typical for photogrammetry point clouds.

The idea to identify those clumps of noise is to construct a surface that is sufficiently close to the ground but always above the noise so that it can be used to classify all points beneath it as noise. However, preserving true ground features without latching onto low noise points often requires several iterations of fine-tuning the parameters. We did this interactively by repeatedly running the processing on only two representative tiles until a desired outcome was achieved.

First we use lasthin to give the point the classification code 8 that is closest to the 20th percentile in elevation within every 90 cm by 90 cm cell (but only if the cells containing at least 20 points). Choosing larger step sizes or higher percentiles resulted in missing ground features. Choosing smaller step sizes or lower percentiles resulted in low noise becoming part of the final ground model.

lasthin -i 2_tiles_raw\fillongley_*.laz ^
        -step 0.90 ^
        -percentile 20 20 ^
        -classify_as 8 ^
        -odir 3_tiles_temp1 -olaz ^
        -cores 4

The we run lasnoise only points on the points with classification code 8 (by ignoring those with classification code 0) and reclassify all “overly isolated” points with code 12. The check for isolation uses cells of size 200 cm by 200 cm by 50 cm and reclassifies the points in the center cell when the surrounding neighborhood of 27 cells has only 3 or fewer points in total. Changing the parameters for ‘-step_xy 2.00 -step_z 0.50 -isolated 3’ will remove noise more or less aggressive.

lasnoise -i 3_tiles_temp1\fillongley_*.laz ^
         -ignore_class 0 ^
         -step_xy 2.00 -step_z 0.50 -isolated 3 ^
         -classify_as 12 ^
         -odir 3_tiles_temp2 -olaz ^
         -cores 4

Next we use lasground to ground-classify only the surviving points (that still have classification code 8) by ignoring those with classification codes 0 or 12 and set their classification code to ground (2) or non-ground (1). The temporary surface defined by the resulting ground points will be used to classify low points as noise in the next step.

lasground -i 3_tiles_temp2\fillongley_*.laz ^
          -ignore_class 0 12 ^
          -town -ultra_fine ^
          -odir 3_tiles_temp3 -olaz ^
          -cores 4

Then we use lasheight to classify all points that are 20 cm or more below the triangulated surface of temporary ground points as points as noise (7) and all others as unclassified (1).

lasheight -i 3_tiles_temp3\fillongley_*.laz ^
          -classify_below -0.20 7 ^
          -classify_above -0.20 1 ^
          -odir 4_tiles_denoised -olaz ^
          -cores 4

The progress of each step is illustrated visually in the two image sequences shown below.

 

 

Now that all noise points are classified we start a standard processing pipeline, but always ignore the noise points that are now classified with classification code 7.

The processing steps below create a 25 cm DTM raster. We first use lasthin to classify the lowest non-noise point per 25 cm by 25 cm cell. Considering only those lowest points we use lasground with options ‘-town’, ‘-extra_fine’, or ‘-bulge 0.1’. Using las2dem the resulting ground points are interpolated into a TIN and rasterized into a 25 cm DTM cutting out only the center 200 meter by 200 meter tile. We store the DTM raster as a gridded LAZ for maximal compression and finally merge these gridded LAZ files to create a hillshaded raster in PNG format with blast2dem.

lasthin -i 4_tiles_denoised\fillongley_*.laz ^
        -ignore_class 7 ^
        -step 0.25 ^
        -lowest ^
        -classify_as 8 ^
        -odir 5_tiles_thinned_lowest -olaz ^
        -cores 4

lasground -i 5_tiles_thinned_lowest\fillongley_*.laz ^
          -ignore_class 1 7 ^
          -town -extra_fine -bulge 0.1 ^
          -odir 6_tiles_ground -olaz ^
          -cores 4

las2dem -i 6_tiles_ground\fillongley_*.laz ^
        -keep_class 2 ^
        -step 0.25 ^
        -use_tile_bb ^
        -odir 7_tiles_dtm -olaz ^
        -cores 4

blast2dem -i 7_tiles_dtm\fillongley_*.laz -merged ^
          -hillshade ^
          -step 0.25 ^
          -o dtm_hillshaded.png

The processing steps below create a 25 cm DSM raster. We first use lasthin to classify the highest non-noise point per 25 cm by 25 cm cell. With las2dem the highest points are interpolated into a TIN and rasterized into a 25 cm DSM cutting out only the center 200 meter by 200 meter tile. Again we store the raster as gridded LAZ for maximal compression and merge these files to create a hillshaded raster in PNG format with blast2dem.

lasthin -i 4_tiles_denoised\fillongley_*.laz ^
        -ignore_class 7 ^
        -step 0.25 ^
        -highest ^
        -classify_as 8 ^
        -odir 8_tiles_thinned_highest -olaz ^
        -cores 4

las2dem -i 8_tiles_thinned_highest\fillongley_*.laz ^
        -keep_class 8 ^
        -step 0.25 ^
        -use_tile_bb ^
        -odir 9_tiles_dsm -olaz ^
        -cores 4

blast2dem -i 9_tiles_dsm\fillongley_*.laz -merged ^
          -hillshade ^
          -step 0.25 ^
          -o dsm_hillshaded.png

The final result is below. The entire script is linked here. Simply download it, modify it as needed, and try it on your data.

 

In Sweden, all they Wanted for Christmas was National LiDAR as Open Data

Let’s heat up some sweet, warm and spicy Glögg in celebration! They must have been good boys and girls up there in Sweden. Because “Jultomten” or simply ”Tomten” – how Sweden’s Santa Clause is called – is assuring a “God Jul” for all the Swedish LiDAR lovers this Christmas season.

Only a few weeks ago this tweet of ours had (mistakenly) included Sweden in a list of European countries that had released their national LiDAR archives as open data for public reuse over the past six years.

Turns out we were correct after all. Sweden has just opened their LiDAR data for free and unencumbered download. To get the data simply create a user account and browse to the ftp site for download as shown in the image sequence below.

The released LiDAR data was collected with a density of 1 to 2 pulses per square meter and is distributed in LASzip compressed LAZ tiles of 2500 by 2500 meters. The returns are classified into four classes: ground (2), water (9), low noise (7) and high noise (18). All items that can not be classified as any of the first four classes coded as left unclassified (1). The LAZ files do not contain CRS information, but this can easily be added with horizontal coordinates in SWERED99 TM (EPSG code 3006) and elevations in RH2000 height (EPSG code 5613).

Below a look with lasview at a 5 km by 5 km area that composed of the four tiles ‘18P001_67100_5800_25.laz‘, ‘18P001_67100_5825_25.laz‘, ‘18P001_67125_5800_25.laz‘ and ‘18P001_67125_5825_25.laz‘ with several of the different color modes available.

 

Some more details: The data was acquired at flying altitude of around 3000 meter with a maximum scan angle of ± 20º and a minimum side overlap of 10% between the flightlines. The laser footprint on ground is below 75 centimeters with slight variation based on the flying altitude. The laser scanning survey was performed with LiDAR instruments that can provide at least three returns from the same pulse. All LiDAR returns are preserved throughout the entire production chain.

The LiDAR data comes with the incredibly Creative Commons – CC0 license, which means that you can use, disseminate, modify and build on the data – even for commercial purposes – without any restrictions. You are free to acknowledge the source when you distribute the data further, but it is not required.

The LiDAR data will eventually cover approximately 75% of Sweden and new point clouds will continuously be added as additional scanning is performed according to the schedule shown below. The survey will be returning to scan every spot again after about 7 years.

2018-2022 LiDAR acquisition plan for Sweden

Below a lasinfo report for tile ‘18P001_67125_5825_25.laz‘. One noticeable oddity is the distribution of intensities. The histogram across all intensities with bins of size 256 shows two clearly distinct sets of intensities each with their own peak and a void of values between 3000 and 10000.

lasinfo -i 18P001_67125_5825_25.laz -cd -histo intensity 256
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            1
  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000
  version major.minor:        1.2
  system identifier:          ''
  generating software:        'TerraScan'
  file creation day/year:     303/2018
  header size:                227
  offset to point data:       227
  number var. length records: 0
  point data format:          1
  point data record length:   28
  number of point records:    20670652
  number of points by return: 13947228 4610837 1712043 358397 42147
  scale factor x y z:         0.01 0.01 0.01
  offset x y z:               0 0 0
  min x y z:                  582500.00 6712500.00 64.56
  max x y z:                  584999.99 6714999.99 136.59
LASzip compression (version 3.2r2 c2 50000): POINT10 2 GPSTIME11 2
reporting minimum and maximum for all LAS point record entries ...
  X            58250000   58499999
  Y           671250000  671499999
  Z                6456      13659
  intensity          32      61406
  return_number       1          5
  number_of_returns   1          5
  edge_of_flight_line 0          1
  scan_direction_flag 0          1
  classification      1         18
  scan_angle_rank   -19         19
  user_data           0          1
  point_source_ID  1802       1804
  gps_time 222241082.251248 222676871.876191
number of first returns:        13947228
number of intermediate returns: 2110980
number of last returns:         13952166
number of single returns:       9339722
covered area in square units/kilounits: 5923232/5.92
point density: all returns 3.49 last only 2.36 (per square units)
      spacing: all returns 0.54 last only 0.65 (in units)
overview over number of returns of given pulse: 9339722 5797676 4058773 1263967 210514 0 0
histogram of classification of points:
        10888520  unclassified (1)
         9620725  ground (2)
           22695  noise (7)
          138147  water (9)
             565  Reserved for ASPRS Definition (18)
intensity histogram with bin size 256.000000
  bin [0,256) has 1753205
  bin [256,512) has 3009640
  bin [512,768) has 2240861
  bin [768,1024) has 1970696
  bin [1024,1280) has 1610647
  bin [1280,1536) has 1285858
  bin [1536,1792) has 974475
  bin [1792,2048) has 790480
  bin [2048,2304) has 996926
  bin [2304,2560) has 892755
  bin [2560,2816) has 164142
  bin [2816,3072) has 57367
  bin [3072,3328) has 18
         [void]
  bin [10752,11008) has 589317
  bin [11008,11264) has 3760
  bin [11264,11520) has 99653
  bin [11520,11776) has 778739
  bin [11776,12032) has 1393569
  bin [12032,12288) has 1356850
  bin [12288,12544) has 533202
  bin [12544,12800) has 140223
  bin [12800,13056) has 16195
  bin [13056,13312) has 2319
  bin [13312,13568) has 977
  bin [13568,13824) has 765
  bin [13824,14080) has 648
  bin [14080,14336) has 289
  bin [14336,14592) has 513
  bin [14592,14848) has 383
  bin [14848,15104) has 178
  bin [15104,15360) has 526
  bin [15360,15616) has 108
  bin [15616,15872) has 263
  bin [15872,16128) has 289
  bin [16128,16384) has 69
  bin [16384,16640) has 390
  bin [16640,16896) has 51
  bin [16896,17152) has 186
  bin [17152,17408) has 239
  bin [17408,17664) has 169
  bin [17664,17920) has 58
  bin [17920,18176) has 227
  bin [18176,18432) has 169
  bin [18432,18688) has 40
  bin [18688,18944) has 401
  bin [18944,19200) has 30
  bin [19200,19456) has 411
  bin [19456,19712) has 34
  bin [19712,19968) has 34
  bin [19968,20224) has 398
  bin [20224,20480) has 24
  bin [20480,20736) has 108
  bin [20736,20992) has 267
  bin [20992,21248) has 29
  bin [21248,21504) has 318
  bin [21504,21760) has 26
  bin [21760,22016) has 59
  bin [22016,22272) has 184
  bin [22272,22528) has 52
  bin [22528,22784) has 18
  bin [22784,23040) has 116
  bin [23040,23296) has 55
  bin [23296,23552) has 89
  bin [23552,23808) has 250
  bin [23808,24064) has 24
  bin [24064,24320) has 52
  bin [24320,24576) has 14
  bin [24576,24832) has 29
  bin [24832,25088) has 71
  bin [25088,25344) has 74
  bin [25344,25600) has 2
  bin [25600,25856) has 17
  bin [25856,26112) has 2
  bin [26368,26624) has 9
  bin [26624,26880) has 1
  bin [26880,27136) has 1
  bin [27136,27392) has 1
  bin [27392,27648) has 1
  bin [27648,27904) has 3
  bin [28416,28672) has 2
  bin [29184,29440) has 4
  bin [30720,30976) has 1
  bin [30976,31232) has 2
  bin [31232,31488) has 1
  bin [32512,32768) has 1
  bin [36864,37120) has 1
  bin [58368,58624) has 1
  bin [61184,61440) has 1
  average intensity 3625.2240208968733 for 20670652 element(s)

CyArk partners with Google, takes over “Don’t be Evil” Mantra, opens LiDAR Archive

One of our most popular (and controversial) blog articles was “Can You Copyright LiDAR“. It was written after we saw the then chief executive director at CyArk commenting “Sweeeet use of CyArk data” on an article describing the creation of a sugary fudge replica of Guatemala’s Tikal temple promoting a series of sugars by multinational agribusiness Tate & Lyle. Yet just a few months earlier our CEO’s university was instructed to take down his Web pages that – using the same data set – were demonstrating how to realize efficient 3D content delivery across the Web. CyArk told university administrators in an email that he was “[…] hosting unauthorized content from CyArk […]”. The full story is here.

Back then, the digital preservation strategy of CyArk was to keep their archaeological scans safe through their partnership with Iron Mountain. In the comment section of “Can You Copyright LiDAR” you can find several entries that are critical of this approach. But that was five years ago. Earlier this year and just after Google removed the “Don’t be Evil” mantra from their code of conduct, CyArk stepped up to take it over and completely changed their tight data control policies. Through their “Open Heritage initiative” CyArk released for the first time their raw LiDAR and imagery with an open license. Here in their own words:

In 2018, CyArk launched the Open Heritage initiative, a
collaboration with Google Arts and Culture to make available
our archive to a broader audience. This was the first time
CyArk has made available primary data sets, including lidar
scans, photogrammetric imagery and corresponding metadata
in a standardized format on a self-serve platform. We are
committed to opening up our archive further as we collect
new data and publishing existing projects where permissions
allow. The data is made available for education, research
and other non-commercial uses via a a Creative Commons
Attribution-NonCommercial 4.0 International License.

This is a HUGE change from the situation in 2013 that resulted in the deletion of our CEO’s Web pages. So we went to download Guatemala’s Tikal temple – the one that got him into trouble back then. It is provided as a single E57 file called ‘Tikal.e57’ with a size of 1074 MB that contains 35,551,759 points in 118 individual scan positions. Using the e572las.exe tool that is part of LAStools we converted this into a single LAZ file ‘Tikal.laz’ with a size of 164 MB.

C:\LAStools\bin>e572las -i c:\data\Tikal\Tikal.e57 ^
                        -o c:\data\Tikal\Tikal.laz

We were not able to find information about the Coordinate Reference System (CRS), but after looking at the coordinate bounding box (see lasinfo report at the end of the article) and the set of projections covering Guatemala, one can make an educated guess that it might be UTM 16 north. Generating a false-colored highest-return 0.5 meter raster with lasgrid and loading it into Google Earth quickly confirms that this is correct.

lasgrid -i c:\data\Tikal\Tikal.laz ^
        -step 0.5 ^
        -highest ^
        -false ^
        -utm 16north ^
        -odix _elev -opng

Now we can laspublish the file with the command line below to create an interactive 3D Web portal using Potree. Unlike five years ago we should now be permitted to create an online portal without the headaches of last time. The CC BY-NC 4.0 license allows to copy and redistribute the material in any medium or format.

laspublish -i c:\data\Tikal\Tikal.laz ^
           -rgb ^
           -utm 16north ^
           -o tikal.html ^
           -title "CyArk's LiDAR Scan of Tikal" ^
           -description "35,551,759 points from 118 individual scans (licensed CC BY-NC 4.0)" ^
           -odir C:\data\Tikal\Tikal -olaz ^
           -overwrite

Below are two screenshots of the online portal that we have just created including some quick distance measurements. This is amazing data. Wow!

Looking at “Templo del Gran Jaguar” from “La Gran Plaza” after taking two measurements.

Overlooking “La Gran Plaza” out of the upper opening of “Templo del Gran Jaguar” with “Templo del las Mascaras” in the back.

We congratulate CyArk to their new Open Heritage initiative and thank them for providing easy access to the Tikal temple LiDAR scans as open data with a useful Creative Commons Attribution-NonCommercial 4.0 International license. Thank you, CyArk, for your contribution to open data and open science. Kudos!

C:\LAStools\bin>lasinfo -i c:\data\Tikal\Tikal.laz
lasinfo (181119) report for 'c:\data\Tikal\Tikal.laz'
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            0
  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000
  version major.minor:        1.2
  system identifier:          'LAStools (c) by Martin Isenburg'
  generating software:        'e572las.exe (version 180919)'
  file creation day/year:     0/0
  header size:                227
  offset to point data:       227
  number var. length records: 0
  point data format:          2
  point data record length:   26
  number of point records:    35551759
  number of points by return: 35551759 0 0 0 0
  scale factor x y z:         0.001 0.001 0.001
  offset x y z:               220000 1900000 0
  min x y z:                  220854.951 1905881.781 291.967
  max x y z:                  221115.921 1906154.829 341.540
LASzip compression (version 3.2r4 c2 50000): POINT10 2 RGB12 2
reporting minimum and maximum for all LAS point record entries ...
  X              854951    1115921
  Y             5881781    6154829
  Z              291967     341540
  intensity       24832      44800
  return_number       1          1
  number_of_returns   1          1
  edge_of_flight_line 0          0
  scan_direction_flag 0          0
  classification      0          0
  scan_angle_rank     0          0
  user_data           0          0
  point_source_ID     1        118
  Color R 0 65280
        G 0 65280
        B 0 65280
number of first returns:        35551759
number of intermediate returns: 0
number of last returns:         35551759
number of single returns:       35551759
overview over number of returns of given pulse: 35551759 0 0 0 0 0 0
histogram of classification of points:
        35551759  never classified (0)

LASmoons: Sebastian Flachmeier

Sebastian Flachmeier (recipient of three LASmoons)
UniGIS Master of Science, University of Salzburg, AUSTRIA
Bavarian Forest National Park, administration, Grafenau, GERMANY

Background:
The Bavarian Forest National Park is located in South-Eastern Germany, along the border with the Czech Republic. It has a total area of 240 km² and its elevation ranges from 600 to 1453 m. In 2002 a project called “High-Tech-Offensive Bayern” was started and a few first/last return LiDAR transects were flown to compute some forest metrics. The results showed that LiDAR has an advantage over other methods, because the laser was able to get readings from below the canopy. New full waveform scanner were developed that produced many more returns in the lower canopy. The National Park experimented with this technology in several projects and improved their algorithms for single tree detection. In 2012 the whole park was flown with full waveform and strategies for LiDAR based forest inventory for the whole National Park were developed. This is the data that is used in the following workflow description.

The whole Bavarian Forest National Park (black line), 1000 meter tiles (black dotted lines), the coverage of the recovered flight lines (light blue). In the area marked yellow within the red frame there are gaps in some of the flightlines. The corresponding imagery in Google Earth shows that this area contains a water reservoir.

Goal:
Several versions of the LiDAR existed on the server of the administration that didn’t have the attributes we needed to reconstruct the original flight lines. The number of returns per pulse, the flight line IDs, and the GPS time stamps were missing. The goal was a workflow to create a LAStools workflow to convert the LiDAR from the original ASCII text files provided by the flight company into LAS or compressed LAZ files with all fields properly populated.

Data:
+
 ALS data flown in 2012 by Milan Geoservice GmbH 650 m above ground with overlap.
+ full waveform sensor (RIEGL 560 / Q680i S) with up to 7 returns per shot
+ total of 11.080.835.164 returns
+ in 1102 ASCII files with *.asc extension (changed to *.txt to avoid confusion with ASC raster)
+ covered area of 1.25 kilometers
+ last return density of 17.37 returns per square meter

This data is provided by the administration of Bavarian Forest National Park. The workflow was part of a Master’s thesis to get the academic degree UniGIS Master of Science at the University of Salzburg.

LAStools processing:

The LiDAR was provided as 1102 ASCII text files named ‘spur000001.txt’ to ‘spur001102.txt’ that looked like this:

more spur000001.txt
4589319.747 5436773.357 685.837 49 106 1 215248.851500
4589320.051 5436773.751 683.155 46 24 2 215248.851500
4589320.101 5436773.772 686.183 66 87 1 215248.851503
[…]

Positions 1 to 3 store the x, y, and z coordinate in meter [m]. Position 4 stores the “echo width” in 0.1 nanoseconds [ns], position 5 stores the intensity, position 6 stores the return number, position 7 stores the GPS time stamp in seconds [s] of the current GPS week. The “number of returns (of given pulse)” information is not explicitly stored and will need to be reconstructed in order, for example, to identify which returns are last returns. The conversion from ASCII text to LAZ was done with the txt2las command line shown below that incorporates these rationals:

  • Although the ASCII files list the three coordinates with millimeter resolution (three decimal digits), we store only centimeter resolution which is sufficient to capture all the precision in a typical airborne LiDAR survey.
  • After computing histograms of the “return number” and the “echo width” for all points with lasinfo and determining their maximal ranges it was decided to use point type 1 which can store up to 7 returns per shot and store the “echo width” as an additional attribute of type 3 (“unsigned short”) using “extra bytes”.
  • The conversion from GPS time stamp in GPS week time to Adjusted Standard time was done by finding out the exact week during which Milan Geoservice GmbH carried out the survey and looking up the corresponding GPS week 1698 using this online GPS time calculator.
  • Information about the Coordinate Reference System “DHDN / 3-degree Gauss-Kruger zone 4” as reported in the meta data is added in form of EPSG code 31468 to each LAS file.
txt2las -i ascii\spur*.txt ^
        -parse xyz0irt ^
        -set_scale 0.01 0.01 0.01 ^
        -week_to_adjusted 1698 ^
        -add_attribute 3 "echo width" "of returning waveform [ns]" 0.1 0 0.1 ^
        -epsg 31468 ^
        -odir spur_raw -olaz ^
        -cores 4

The 1102 ASCII files are now 1102 LAZ files. Because we switched from GPS week time to Adjusted Standard GPS time stamps we also need to set the “global encoding” flag in the LAS header from 0 to 1 (see ASPRS LAS specification). We can do this in-place (i.e. without creating another set of files) using the following lasinfo command:

lasinfo -i spur_raw\spur*.laz ^
        -nh -nv -nc ^
        -set_global_encoding 1

To reconstruct the missing flight line information we look for gaps in the sequence of GPS time stamps by computing GPS time histograms with lasinfo and bins of 10 seconds in size:

lasinfo -i spur_raw\spur*.laz -merged ^
        -histo gps_time 10 ^
        -o spur_raw_all.txt

The resulting histogram exhibits the expected gaps in the GPS time stamps that happen when the survey plane leaves the target area and turns around to approach the next flight line. The subsequent histogram entries marked in red show gaps of 120 and 90 seconds respectively.

more spur_raw_all.txt
[...]
bin [27165909.595196404,27165919.595196255) has 3878890
bin [27165919.595196255,27165929.595196106) has 4314401
bin [27165929.595196106,27165939.595195957) has 435788
bin [27166049.595194317,27166059.595194168) has 1317998
bin [27166059.595194168,27166069.595194019) has 4432534
bin [27166069.595194019,27166079.59519387) has 4261732
[...]
bin [27166239.595191486,27166249.595191337) has 3289819
bin [27166249.595191337,27166259.595191188) has 3865892
bin [27166259.595191188,27166269.595191039) has 1989794
bin [27166349.595189847,27166359.595189698) has 2539936
bin [27166359.595189698,27166369.595189549) has 3948358
bin [27166369.595189549,27166379.5951894) has 3955071
[...]

Now that we validated their existence, we use these gaps in the GPS time stamps to split the LiDAR back into the original flightlines it was collected in. Using lassplit we produce one file per flightline as follows:

lassplit -i spur_raw\spur*.laz -merged ^
         -recover_flightlines_interval 10 ^
         -odir strips_raw -o strip.laz

In the next step we repair the missing “number of returns (per pulse)” field that was not provided in the ASCII file. This can be done with lasreturn assuming that the point records in each file are sorted by increasing GPS time stamp. This happens to be true already in our case as the original ASCII files where storing the LiDAR returns in acquisition order and we have not changed this order. If the point records are not yet in this order it can be created with lassort as follows. As these strips can have many points per file it may be necessary to run the new 64 bit executables by adding ‘-cpu64’ to the command line in order to avoid running out of memory.

lassort -i strips_raw\strips*.laz ^
        -gpstime -return_number ^
        -odir strips_sorted -olaz ^
        -cores 4 -cpu64

An order sorted by GPS time stamp is necessary as lasreturn expects point records with the same GPS time stamp (i.e. returns generated by the same laser pulse) to be back to back in the input file. To ‘-repair_number_of_returns’ the tool will load all returns with the same GPS time stamp  and update the “number of returns (per pulse)” attribute of each return to the highest “return number” of the loaded set.

lasreturn -i strips_sorted\strips*.laz ^
          -repair_number_of_returns ^
          -odir strips_repaired -olaz ^
          -cores 4

In a final step we use las2las with the ‘-files_are_flightlines’ option (or short ‘-faf’) to set the “file source ID” field in the LAS header and the “point source ID” attribute of every point record in the file to the same unique value per strip. The first file in the folder will have all its field set to 1, the next file will have all its field set to 2, the next file to 3 and so on. Please do not run this on multiple cores for the time being.

las2las -i strips_repaired\strips*.laz ^
        -files_are_flightlines ^
        -odir strips_final -olaz

It’s always useful to run a final validation of the files using lasvalidate to reassure yourself and the people you will be sharing the data with that nothing funky has happened during any of these conversion steps.

lasvalidate -i strips_final\strip*.laz ^
            -o strips_final\report.xml

And it can also be useful to add an overview in SHP or KML format to the delivery that can be created with lasboundary as follows:

lasboundary -i strips_final\strip*.laz ^
            -overview -labels ^
            -o strips_final\overview.kml

The result was 89 LAZ files (each containing one complete flightline) totaling 54 GB compared to 1102 ASCII files (each containing a slice of a flightline) totaling 574 GB.