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.

National Open LiDAR Strategy of Latvia humiliates Germany, Austria, and other European “Closed Data” States

Latvia, officially the Republic of Latvia, is a country in the Baltic region of Northern Europe has around 2 million inhabitants, a territory of 65 thousand square kilometers and Рsince recently Рalso a fabulous open LiDAR policy. Here is a list of 65939 tiles in LAS format available for free download that cover the entire country with airborne LiDAR with a density from 4 to 6 pulses per square meters. The data is classified into ground, building, vegetation, water, low noise, and a few other classifications. It is licensed Creative Commons CC0 1.0 Рmeaning that you can copy, modify, and distribute the data, even for commercial purposes, all without asking permission. And there is a simple and  functional interactive download portal where you can easily download individual tiles.

latvia_open_data_portal_01

Interactive open LiDAR download portal of Latvia.

We downloaded the 5 by 5 block of¬†square kilometer¬†tiles matching “4311-32-XX.las” for checking the quality and creating a 1m DTM and a 1m DSM raster. You can follow along after downloading the latest version of LAStools.

Quality Checking

We first run lasvalidate and lasinfo on the downloaded LAS files and then immediately compress them with laszip because multi-core processing of uncompressed LAS files will quickly overwhelm our file system, make processing I/O bound, and result in overall longer processing times with CPUs waiting idly for data to be loaded from the drives.

lasinfo -i 00_tiles_raw\*.las ^
        -compute_density ^
        -histo z 5 ^
        -histo intensity 256 ^
        -histo user_data 1 ^
        -histo scan_angle 1 ^
        -histo point_source 1 ^
        -histo gps_time 10 ^
        -odir 01_quality -odix _info -otxt ^
        -cores 3
lasvalidate -i 00_tiles_raw\*.las ^
            -no_CRS_fail ^
            -o 01_quality\report.xml

Despite already excluding a missing Coordinate Reference System (CRS) from being a reason to fail (the lasinfo reports show that the downloaded LAS files do not have any geo-referencing information) lasvalidate still reports a few failing files, but scrutinizing the resulting XML file ‘report.xml’ shows only minor issues.

Usually during¬†laszip compression we do not alter the contents of a file, but here we also add the EPSG code 3059 for CRS “LKS92 / Latvia TM” as we turn bulky LAS files into slim LAZ files so we don’t have to specify it in all future processing steps.

laszip -i 00_tiles_raw\*.las ^
       -epsg 3059 ^
       -cores 2

Compression reduces the total size of the 25 tiles from over 4.1 GB to below 0.6 GB.

Next we use lasgrid to visualize the last return density which corresponds to the pulse density of the LiDAR survey. We map each 2 by 2 meter pixel where the last return density is 2 or less to blue and each 2 by 2 meter pixel it is 8 or more to red.

lasgrid -i 00_tiles_raw\*.laz ^
        -keep_last ^
        -step 2 ^
        -density_16bit ^
        -false -set_min_max 2 8 ^
        -odir 01_quality -odix _d_2_8 -opng ^
        -cores 3

This we follow by the mandatory lasoverlap check for flight line overlap and alignment where we map the number of overlapping swaths as well as the worst vertical difference between overlapping swaths to a color that allows for quick visual quality checking.

lasoverlap -i 00_tiles_raw\*.laz ^
           -step 2 ^
           -min_diff 0.1 -max_diff 0.2 ^
           -odir 01_quality -opng ^
           -cores 3

The results of the quality checks with lasgrid and lasoverlap are shown below.

Raster Derivative Generation

Now we use first las2dem to create a Digital Terrain Model (DTM) and a Digital Surface Model (DSM) in RasterLAZ format and then use blast2dem to create merged and hill-shaded versions of both. Because we will use on-the-fly buffering to avoid edge effects along tile boundaries we first spatially index the data using lasindex for more efficient access to the points from neighboring tiles.

lasindex -i 00_tiles_raw\*.laz ^
         -cores 3

las2dem -i 00_tiles_raw\*.laz ^
        -keep_class 2 9 ^
        -buffered 25 ^
        -step 1 ^
        -use_orig_bb ^
        -odir Latvia\02_dtm_1m -olaz ^
        -cores 3

blast2dem -i 02_dtm_1m\*.laz ^
          -merged ^
          -hillshade ^
          -step 1 ^
          -o dtm_1m.png

las2dem -i 00_tiles_raw\*.laz ^
        -drop_class 1 7 ^
        -buffered 10 ^
        -spike_free 1.5 ^
        -step 1 ^
        -use_orig_bb ^
        -odir 03_dsm_1m -olaz ^
        -cores 3

blast2dem -i 03_dsm_1m\*.laz ^
          -merged ^
          -hillshade ^
          -step 1 ^
          -o dsm_1m.png

Because the overlaid imagery does not look as nice in our new Google Earth installation, below are the DTM and DSM at versions down-sampled to 25% of their original size.

Many thanks to SunGIS from Latvia who tweeted us about the Open LiDAR after we chatted about it during the Foss4G 2019 gala dinner. Kudos to the Latvian Geospatial Information Agency (LGIA) for implementing a modern national geospatial policy that created opportunity for maximal return of investment by opening the expensive tax-payer funded LiDAR data for re-purposing and innovation without barriers. Kudos!

Removing Low Noise in LiDAR Points with Median Ground Surface

Recently a user of LAStools asked a question in our user forum about how to classify LiDAR data that contains lots of low noise. A sample screen shot of the user’s failed attempt to correctly classify the noise using lasnoise and the ground with lasground is shown below: red points are noise, brown points are ground, and grey points are unclassified. In this article we show how to remove this low noise using a temporary ground surface that we construct from a subset of points at a certain elevation percentile. You can follow along by downloading the data and the sequence of command lines used.

example of miss-classified low noise points: ground points (brown) below ground

Download the LiDAR data set that was apparently flown with a RIEGL “crossfire” Q1560. You can also download the command line sequence here. We first run lasinfo with option ‘-compute_density’ (or ‘-cd’ for short) to get a rough idea about the last return density which is quite high with an average of over 31 last returns per square meter. We then use lasthin to classify one last return per square meter with the temporary classification code 8, namely the one whose elevation is closest to the 20th percentile per 1 meter by 1 meter grid cell. We then repeat this command line for the 30th, 40th, 50th percentile modifying the command line accordingly. You must use this version of lasthin that will part of a future LAStools release as options ‘-ignore_first_of_many’ and ‘-ignore_intermediate’ were just added this weekend.

lasthin -i crossfire.laz ^
        -ignore_first_of_many -ignore_intermediate ^
        -step 1 ^
        -percentile 20 15 ^
        -classify_as 8 ^
        -odix _p20 -olaz

Below you see the resulting subset of points marked with the temporary classification code 8 for the four different percentiles 20th, 30th, 40th, and 50th triangulated into a surface and hill-shaded.

Next we reclassify only those points marked with the temporary classification code 8 into ground (2) and unclassified (1) points using lasground by ignoring all points that still have the original classification code 0.

lasground -i crossfire_p20.laz ^
          -ignore_class 0 ^
          -wilderness ^
          -odix g -olaz

Below you see the resulting ground points computed from the subsets of points at four different percentiles 20th, 30th, 40th, and 50th triangulated into a surface and hill-shaded.

Both the ground classification of the 40th and the 50th percentile look reasonable. Only a few down spikes remain in the 40th percentile surface and a few additional bumps appear in the 50th percentile surface. Next we use lasheight with those two reasonable-looking ground surfaces to classify all points that are 20 centimeter below the triangulated ground surface into the noise classification code 7.

lasheight -i crossfire_p40g.laz ^
          -classify_below -0.2 7 ^
          -do_not_store_in_user_data ^
          -odix h -olaz

Now that the low noise points were removed (or rather classified as noise) we start the actual ground classification process. In this example we want to create a 50 cm DTM, hence it is more than sufficient to find one ground point per 25 cm cell. Therefore we first move all lowest non-noise last return per 25 cm cell to the temporary classification code 8.

Side note: One might also consider to modify the following workflow to run the ground classification on more than just the last returns by omitting ‘-ignore_first_of_many’ and ‘-ignore_intermediate’ from the lasthin call and by adding ‘-all_returns’ to the lasground call. Why? Because for all laser shots that resulted in a low noise point, this noise point will usually be the last return, so that the true ground hit could be the second to last return.

lasthin -i crossfire_p40gh.laz ^
        -ignore_first_of_many -ignore_intermediate ^
        -ignore_class 7 ^
        -step 0.25 ^
        -lowest ^
        -classify_as 8 ^
        -odix _low25 -olaz

The final ground classification is obtained by running lasground only on the points with temporary classification code 8 by ignoring all others, namely the noise points (7) and the unclassified points (0 and 1).

lasground -i crossfire_p40gh_low25.laz ^
          -ignore_class 0 1 7 ^
          -wilderness ^
          -odix g -olaz

We then use las2dem to create the 50 cm DTM from the points classified as ground. We store this DTM raster to the LAZ format which has shown to be the most efficient format for storing elevation or height rasters. We have started calling this format RasterLAZ. It is supported by all LAStools and the new DEMzip tool. One advantage is that we can feed RasterLAZ directly back into LAStools, for example as done below, for a second call to las2dem that computes a hill-shaded DTM.

las2dem -i crossfire_p40gh_low25g.laz ^
        -keep_class 2 ^
        -step 0.5 ^
        -ocut 9 -odix _dtm50 -olaz

las2dem -i crossfire_p40_dtm50.laz ^
        -step 0.5 ^
        -hillshade ^
        -odix _hill -opng

Below the resulting hill-shaded DTMs computed for the 40th and the 50th elevation percentile – as well as for the 45th elevation percentile that we’ve added for comparison.

Below we finally take a closer look at an example 1 meter profile line through the LiDAR classified by the 45th percentile workflow. There is a small stretch of ground points that was incorrectly classified as noise points (find the mouse cursor) so it might be worthwhile to change parameters slightly to make the noise classification less aggressive.

Side note follow-up: The return coloring shows there are indeed some ‘intermediate’ as well some ‘first of many returns’ just where we expect the bare terrain to be. However, there are not so many that the results can be expected to drastically change by including them into the ground finding process.

Clean DTM from Agisoft Photogrammetric Points of Urban Scene

We occasionally get permission to distribute a nice data sets and blog about how to best process it with LAStools because this gets around having to pay our “outrageous” consulting fees. (-: This time we received a nice photogrammetric point cloud of the Tafawa Balewa Square in Lagos Island, Lagos, Nigeria. This area is part of the central business district of Lagos and characterized by high-rise buildings. The Tafawa Balewa Square was constructed in 1972 over the site of a defunct track for horse racing and is bounded by Awolowo road, Cable street, Force road, Catholic Mission street and the 26-story Independence House. We want to create a nice Digital Terrain Model from the dense-matching point cloud that was generated with Photoscan by AgiSoft and – as always with photogrammetry – we have to take special care of low noise points. The final result is shown below. All processing commands used are here.

After downloading the data it is useful to familiarize yourself with the file, which can be done with lasview, lasinfo, and lasgrid using the command lines shown below. According to the lasinfo report there are around 47 million points points with RGB colors in the file and their average density is around 100 points per square meter.

lasview -i 0_raw\TafawaBalewa.laz

lasinfo -i 0_raw\TafawaBalewa.laz ^
        -cd -histo intensity 256 ^
        -histo z 1 ^
        -odir 1_quality -odix _info -otxt

lasgrid -i 0_raw\TafawaBalewa.laz ^
        -step 1 ^
        -density ^
        -false -set_min_max 50 150 ^
        -odir 1_quality -odix _d_50_150 -opng

The average point density value of 100 from the lasinfo report suggests that 50 as the minimum and 150 as the maximum are good false color ramp values for a map showing how the point density per square meter is distributed.

Color-coded point density: blue equals 50 or less and red means 150 or more points per square meter.

We use lastile to create a buffered tiling for the 47 million points. We use a tile size of 200 meters and request a large buffer of 50 meters around every tile because there are large buildings in the survey areas. We also flag buffer points as withheld so they can be easily be dropped later.

lastile -i 0_raw\TafawaBalewa.laz ^
        -tile_size 200 -buffer 50 -flag_as_withheld ^
        -odir 2_tiles_raw -o tafawa.laz

If you inspect the resulting tiles – such as ‘tafawa_544000_712600.laz’ as shown here – with lasview you will see the kind of low noise that is shown below and that may cause a ground classification algorithm. While our lasground software is able to deal with a certain amount of low noise – if there are too many it will likely latch onto them. Therefore we will first generate a subset of points that has as few as possible of such low noise points.

Typical low noise in dense-matching photogrammetry points in urban scene.

Next we use a sequence of three LAStools modules, namely lasthin, lasground, and lasheight to classify this photogrammtric point cloud into ground and non-ground points. All processing commands used are here. First we use lasthin to give the point the classification code 8 that is closest to the 50th percentile in elevation within every 50 centimeter by 50 centimeter cell (but only if the cells containing at least 20 points).

lasthin -i 2_tiles_raw\tafawa*.laz ^
        -step 0.5 ^
        -percentile 50 20 ^
        -classify_as 8 ^
        -odir 3_tiles_median_50cm -olaz ^
        -cores 3

Next we use¬†lasground to ground-classify only the points that have classification code 8 (i.e. by ignoring those with classification codes 0) and set their classification code to ground (2) or non-ground (1). Because of the large buildings in this urban scene we use ‘-metro’ which uses a large step size of 50 meters for the pre-processing. This also sets the internally used bulge parameter to 5.0 which you can see if you run the tool in verbose ‘-v’ mode. In three different trial runs we determined that forcing the bulge parameter to be 0.5 instead gave better results. The bulge and the spike parameters can be useful to vary in order to improve ground classification results (also see the README file).

lasground -i 3_tiles_median_50cm\tafawa*.laz ^
          -ignore_class 0 ^
          -metro -bulge 0.5 ^
          -odir 4_tiles_ground_50cm -olaz ^
          -cores 3

The resulting ground points are a subset with a resolution of 50 centimeter that is good enough to create a DTM with meter resolution, which we do with las2dem command line shown below. We really like storing DTM elevation rasters to the LAZ point format because it is a more compressed way of storing elevation rasters compared to ASC, BIL, TIF, or IMG. It also makes the raster output a natural input to subsequent LAStools processing steps.

las2dem -i 4_tiles_ground_50cm\tafawa*.laz ^
        -keep_class 2 ^
        -step 1 -kill 100 ^
        -use_tile_bb ^
        -odir 5_tiles_dtm_1m -olaz ^
        -cores 3

Finally we use blast2dem to create a seamless hill-shaded version of our 1 meter DTM from on-the-fly merged elevation rasters. This is the DTM pictured at the beginning of this article.

blast2dem -i 5_tiles_dtm_1m\tafawa*.laz -merged ^
          -step 1 ^
          -hillshade ^
          -o dtm_1m.png

The corresponding DSM pictured at the beginning of this article was generated with the two command lines below by first keeping only the 95th percentile highest elevation of every 50 cm by 50 cm cell with lasthin (which remove spurious high noise points) and then by triangulating the surviving points with blast2dem into a seamless TIN that is also hill-shaded and rasterized with 1 meter resolution. Running the 64 bit version of lasthin (note the ‘-cpu64‘ switch) allows us to work on the entire data set (rather than its tiles version) at once, where the standard 32 bit version may run out of memory.

lasthin -i 0_raw\TafawaBalewa.laz ^
        -cpu64 ^
        -step 0.5 ^
        -percentile 95 20 ^
        -o 0_raw\TafawaBalewa_p95_50.laz

blast2dem -i 0_raw\TafawaBalewa_p95_50.laz ^
          -step 1 ^
          -hillshade ^
          -o dsm_1m.png

In order to generate the final DTM at higher resolution we use lasheight to pull all points into the ground class that lie within a 5 cm distance vertically below or a 15 cm distance vertically above the triangulated surface of ground points computed in the previous step. You could experiment with other values here to be less or more conservative about pulling detail into the ground class.

lasheight -i 4_tiles_ground_50cm\tafawa*.laz ^
          -classify_between -0.05 0.15 2 ^
          -odir 6_tiles_ground -olaz ^
          -cores 3

We repeat the same processing step as before las2dem to create the raster DTM tiles, but this time with a resolution of 25 cm.

las2dem -i 6_tiles_ground\tafawa*.laz ^
        -keep_class 2 ^
        -step 0.25 -kill 100 ^
        -use_tile_bb ^
        -odir 7_tiles_dtm_25cm -olaz ^
        -cores 3

And we again use blast2dem to create a seamless hill-shaded version of the DTM from on-the-fly merged elevation rasters, but this time with a resolution of 25 cm. This is the DTM shown below. All processing commands used are here.

blast2dem -i 7_tiles_dtm_25cm\tafawa*.laz -merged ^
          -step 0.25 ^
          -hillshade ^
          -o dtm_25cm.png

Hill-shade of final DTM with resolution of 25 cm.

Smooth DTM from Drone LiDAR off Velodyne HDL 32A mounted on DJI M600 UAV

Recently we attempted to do a small LiDAR survey by drone for a pet project of our CEO in our “code and surf camp” here in Samara, Costa Rica. But surveying is difficult when you are a novice and we ran into a trajectory issue. The dramatic “wobbles” were entirely our fault, but fortunately our mistakes also led to something useful: We found some LAS export bugs. Our laser scanner was a Velodyne HDL-32E integrated with a NovAtel INS into the Snoopy Series A HD made by LiDARUSA. The system was carried by a DJI Matrice 600 (M600) drone. We processed the trajectory with NovAtel Inertial Explorer (here we made the “wobbles” error) and finally exported the LAS and LAZ files with ScanLook PC (version 1.0.182) from LiDARUSA.

While we were investigating our “wobbles” (which clearly were our mistake) we also found five different LAS export bugs in ScanLook PC that seem to have started sometime after version 1.0.171 and will likely end with version 1.0.193. Below an illustration of a correct export from version 1.0.129 and a buggy export from version 1.0.182. In both instances you see the returns from one revolution of the Velodyne HDL-32E scanner head ordered by their GPS time stamps and colored to distinguish the 32 separate beams. In the buggy version, groups of around seven non-adjacent returns are given the same time stamp. This bug will only affect you, if correct GPS time stamps are important for your subsequent LiDAR processing or if your client explicitly asked for ASPRS specification compliant LAS files. We plan to publish another blog post detailing how to find this GPS time stamping bug (and the other four bugs we found).

During the many interactions we had working through “wobbles” and export bugs, we obtained a nice set of six flight lines from Seth Gulich of Bowman Consulting – a US American company based in Stuart, Florida – who flew an identical “Snoopy Series A HD” system also on a DJI Matrice 600 drone at approximately 100 feet above ground level above a model airplane airport in Palm Beach, Florida. You can download the data set here. In the following we will check the flight line alignment of this data set and then process it into a smooth DTM. All command lines used are summarized in this text file.

First we generate a lasinfo report that includes a number of histograms for on-the-fly merged flight lines with lasinfo and then use the z coordinate histogram from the lasinfo report to set reasonable min/max values for the elevation color ramp of lasview:

lasinfo -i 0_strips_raw\Velodyne*.laz -merged ^
        -cd ^
        -histo z 1 ^
        -histo user_data 1 ^
        -histo point_source 1 ^
        -o 1_quality\Velodyne_merged_info.txt

lasview -i 0_strips_raw\Velodyne*.laz ^
        -points 10000000 ^
        -set_min_max 25 75

The lasinfo report shows no information about the coordinate reference system. We found out experimentally that the horizontal coordinates seem to be EPSG code 2236 and that the vertical units are most likely be US survey feet. The warnings you will see in the lasinfo report have to do with the fact that the double-precision bounding box stored in the LAS header was populated with numbers that have many more decimal digits than the coordinates in the file, which only have millifeet resolution as all three scale factors are 0.001 (meaning coordinates have three decimal digits). The information which of the 32 lasers was collecting which point is stored in both the ‘user data’ and the ‘point source ID’ field which is evident from the histograms in the lasinfo report. We need to be careful not to override both fields in later processing.

Next we use lasoverlap to check how well the LiDAR points from the flight out and the flight back align vertically. This tool computes the difference of the lowest points for each square foot covered by multiple flight lines. Differences of less than a quarter of a foot are both times mapped to white, differences of more than one foot (more than half a foot) are mapped to saturated red or blue depending on whether the difference is positive or negative in the first run (in the second run):

lasoverlap -i 0_strips_raw\Velodyne*.laz ^
           -faf ^
           -min_diff 0.25 -max_diff 1.00 -step 1 ^
           -odir 1_quality -o overlap_025_100.png

lasoverlap -i 0_strips_raw\Velodyne*.laz ^
           -faf ^
           -min_diff 0.25 -max_diff 0.50 -step 1 ^
           -odir 1_quality -o overlap_025_050.png

We use a new feature of the LAStools GUI (as of version 180429) to closer inspect large red or blue areas. With lasmerge we clip out regions that looks suspect for closer examination with lasview. First we spatially index the flight lines to make this process faster. With the ‘-gui’ switch we start the tool in GUI mode with flight lines already loaded. Using the new PNG overlay roll-out on the left we add the ‘overlap_025_050_diff.png’ image from the quality folder created in the last step and clip out three areas.

lasindex -i 0_strips_raw\Velodyne*.laz
         -tile_size 10 -maximum -100 ^
         -cores 3

lasmerge -i 0_strips_raw\Velodyne*.laz -gui

You can also clip out these three areas using the command lines below:

lasmerge -i 0_strips_raw\Velodyne*.laz ^
         -faf ^
         -inside_tile 939500 889860 100 ^
         -o 1_quality\939500_889860.laz

lasmerge -i 0_strips_raw\Velodyne*.laz ^
         -faf ^
         -inside_tile 940400 889620 100 ^
         -o 1_quality\940400_889620.laz

lasmerge -i 0_strips_raw\Velodyne*.laz ^
         -faf ^
         -inside_tile 940500 890180 100 ^
         -o 1_quality\940500_890180.laz

The reader may inspect the areas 939500_889860.laz, 940400_889620.laz, and 940500_890180.laz with lasview using profile views via hot keys ‘x’ and switching back and forth between the points from different flight lines via hot keys ‘0’, ‘1’, ‘2’, ‘3’, … for individual and ‘a’ for all flight lines as we have done it in previous tutorials [1,2,3]. Using drop-lines or rise-lines via the pop-up menu gives you a sense of scale. Removing points with lastrack that are horizontally too far from the trajectory could be one strategy to use fewer outliers. But as our surfaces are expected to be “fluffy” (because we have a Velodyne LiDAR system), we accept these flight line differences and continue processing.

Here the complete LAStools processing pipeline for creating an average ground model from the set of six flight lines that results in the hillshaded DTM shown below. The workflow is similar to those we have developed in earlier blog posts for Velodyne Puck based systems like the Hovermap and the Yellowscan and in the other Snoopy tutorial. All command lines used are summarized in this text file.

Hillshaded DTM with half foot resolution generated via average ground computation with LAStools.

In the first step we lastile the six flight lines into 250 by 250 feet tiles with 25 feet buffer while preserving flight line information. The flight line information will be stored in the “point source ID” field of each point and therefore override the beam ID that is currently stored there. But the beam ID is also stored in the “user data” field as the¬† lasinfo report had told us. We set all classifications to zero and add information about the horizontal coordinate reference system EPSG code 2236 and the vertical units (US Survey Feet).

lastile -i 0_strips_raw\*.laz ^
        -faf ^
        -set_classification 0 ^
        -epsg 2236 -elevation_survey_feet ^
        -tile_size 250 -buffer 25 -flag_as_withheld ^
        -odir 2_tiles_raw -o pb.laz

On three cores in parallel we then lassort the points in the tiles into a space-filling curve order which will accelerate later operations.

lassort -i 2_tiles_raw\*.laz ^
        -odir 2_tiles_sorted -olaz ^
        -cores 3

Next we use lasthin to classify the point whose elevation is closest to the 5th elevation percentile among all points falling into its cell with classification code 8. We run lasthin multiple times and each time increase the cell size from 1, 2, 4, 8 to 16 foot. We do this because we have requested the 5th elevation percentile to only be computed when there are at least 20 points in the cell. Percentiles are statistical measures and need a reasonable sample size to be stable. Because drone flights are very dense in the center and more sparse at the edges this increase in cell size assures that we have a good selection of points classified with classification code 8 across the entire survey area.

lasthin -i 2_tiles_sorted\*.laz ^
        -step 1 -percentile 5 20 -classify_as 8 ^
        -odir 3_tiles_thinned_p05_step01 -olaz ^
        -cores 3

lasthin -i 3_tiles_thinned_p05_step01\*.laz ^
        -step 2 -percentile 5 20 -classify_as 8 ^
        -odir 3_tiles_thinned_p05_step02 -olaz ^
        -cores 3

lasthin -i 3_tiles_thinned_p05_step02\*.laz ^
        -step 4 -percentile 5 20 -classify_as 8 ^
        -odir 3_tiles_thinned_p05_step04 -olaz ^
        -cores 3

lasthin -i 3_tiles_thinned_p05_step04\*.laz ^
        -step 8 -percentile 5 20 -classify_as 8 ^
        -odir 3_tiles_thinned_p05_step08 -olaz ^
        -cores 3

lasthin -i 3_tiles_thinned_p05_step08\*.laz ^
        -step 16 -percentile 5 20 -classify_as 8 ^
        -odir 3_tiles_thinned_p05_step16 -olaz ^
        -cores 3

Then we let lasground_new run on only the points classified with classification code 8 (i.e. by ignoring the points still classified with code 0) which classifies them into ground (code 2) and non-ground (code 1).

lasground_new -i 3_tiles_thinned_p05_step16\*.laz ^
              -ignore_class 0 ^
              -town ^
              -odir 4_tiles_ground_low -olaz ^
              -cores 3

The ground points we have computed form somewhat of a lower envelope of the “fluffy” points of a Velodyne scanner. With lasheight we now draw all the points near the ground – namely those from 0.1 foot below to 0.4 foot above the ground – into a new classification code 6 that we term “thick ground”. The ‘-do_not_store_in_user_data’ switch prevent the default behavior of lasheight from happening, which would override the beam ID information that it stored in the ‘user data’ field with approximate height value.

lasheight -i 4_tiles_ground_low\*.laz ^
          -classify_between -0.1 0.4 6 ^
          -do_not_store_in_user_data ^
          -odir 4_tiles_ground_thick -olaz ^
          -cores 3

A few close-up shots of the resulting “thick ground” are shown in the picture gallery below.

We then use lasgrid to average the (orange) thick ground points onto a regular grid with a cell spacing of half a foot. We do not grid the tile buffers by adding the ‘-use_tile_bb’ switch.

lasgrid -i 4_tiles_ground_thick\*.laz ^
        -keep_class 6 ^
        -step 0.5 -average ^
        -use_tile_bb ^
        -odir 5_tiles_gridded_mean_ground -olaz &
        -cores 3

Finally we use blast2dem to merge all the averaged ground point grids into one file, interpolate across open areas without ground points, and compute the hillshaded DTM shown above. All command lines used are summarized in this text file.

blast2dem -i 5_tiles_gridded_mean_ground\*.laz ^
          -merged ^
          -step 0.5 ^
          -hillshade ^
          -o dtm.png

We thank Seth Gulich of Bowman Consulting for sharing this LiDAR data set with us. It was flown with a DJI Matrice 600 drone carrying a “Snoopy A series HD” LiDAR system from LidarUSA.

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.

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, lasnoise, lasground, 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.