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!

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.

Using Open LiDAR to Remove Low Noise from Photogrammetric UAV Point Clouds

We collected drone imagery of the¬†restored “Kance” tavern¬†during the lunch stop of the UAV Tartu summer school field trip (actually the organizer Marko Kohv¬†did that, as I was busy introducing students to SUP boarding). With Agisoft PhotoScan we then processed the images into point clouds below the deck of the historical “J√Ķmmu”¬†barge on the way home (actually Marko did that, because I was busy enjoying the view of the wetlands in the afternoon sun). The resulting data set with 7,855,699 points is shown below and can be downloaded here.

7,855,699 points produced with Agisoft Photoscan

Generating points using photogrammetric techniques in scenes containing water bodies tends to be problematic as dense blotches of noise points above and below the water surface are common as you can see in the picture below. Especially the low points are troublesome as they adversely affect ground classification which results in poor Digital Elevation Models (DTMs).

Clusters of low noise points nearly 2 meters below the actual surface in water areas.

In a previous article we have described a LAStools workflow that can remove excessive low noise. In this article here we use external information about the topography of the area to clean our photogrammetry points. How convenient that the Estonian Land Board has just released their entire LiDAR archives as open data.

Following these instructions here you can download the available open LiDAR for this area, which has the map sheet index 475681. Alternatively you can download the four currently available data sets here flown in spring 2010, in summer 2013, in spring 2014, and in summer 2017. In the following we will use the one flown in spring 2014.

We can view both data sets simultaneously in lasview. By adding ‘-faf’ to the command-line we can switch back and forth between the two data sets by pressing ‘0’ and ‘1’.

lasview -i Kantsi.laz ^
        -i 475681_2014_tava.laz ^
        -points 10000000 ^
        -faf

We find cut the 1 km by 1 km LiDAR tile down to a 250 m by 250 m tile that nicely surrounds our photogrammetric point set using the following las2las command-line:

las2las -i 475681_2014_tava.laz ^
        -inside_tile 681485 6475375 250 ^
        -o LiDAR_Kantsi.laz

lasview -i Kantsi.laz ^
        -i LiDAR_Kantsi.laz ^
        -points 10000000 ^
        -faf

Scrutinizing the two data sets we quickly find that there is a miss-alignment between the dense imagery-derived and the comparatively sparse LiDAR point clouds. With lasview we investigate the difference between the two point clouds by hovering over a point from one point cloud and pressing <i> and then hovering over a somewhat corresponding point from the other point cloud and pressing <SHIFT>+<i>. We measure displacements of around 2 meters vertically and of around 3 to 3.5 meter in total.

Before we can use the LiDAR points to remove the low noise from the photogrammetric points we must align them properly. For simple translation errors this can be done with a new feature that was recently added to lasview. Make sure to download the latest version (190404 or newer) of LAStools to follow the steps shown in the image sequence below.

las2las -i Kantsi.laz ^
        -translate_xyz 0.89 -1.90 2.51 ^
        -o Kantsi_shifted.laz

lasview -i Kantsi_shifted.laz ^
        -i LiDAR_Kantsi.laz ^
        -points 10000000 ^
        -faf

The result looks good in the sense that both sides of the photogrammetric roof are reasonably well aligned with the LiDAR. But there is still a shift along the roof so we repeat the same thing once more as shown in the next image sequence:

We use a suitable displacement vector and apply it to the photogrammetry points, shifting them again:

las2las -i Kantsi_shifted.laz ^
        -translate_xyz -1.98 -0.95 0.01 ^
        -o Kantsi_shifted_again.laz

lasview -i Kantsi_shifted_again.laz ^
        -i LiDAR_Kantsi.laz ^
        -points 10000000 ^
        -faf

The result is still not perfect as there is also some rotational error and you may find another software such as Cloud Compare more suited to align the two point clouds, but for this exercise the alignment shall suffice. Below you see the match between the photogrammetry points and the LiDAR TIN before and after shifting the photogrammetry points with the two (interactively determined) displacement vectors.

The final steps of this exercise use las2dem and the already ground-classified LiDAR compute a 1 meter DTM, which we then use as input to lasheight. We classify the photogrammetry points using their height above this set of ground points with 1 meter spacing: points that are 40 centimeter or more below the LiDAR DTM are classified as noise (7), points that are between 40 below to 1 meter above the LiDAR DTM are classified to a temporary class (here we choose 8) that has those points that could potentially be ground points. This will help, for example, with subsequent ground classification as large parts of the photogrammetry points – namely those on top of buildings and in higher vegetation – can be ignored from the start by a ground classification algorithm such as lasground.

las2dem -i LiDAR_Kantsi.laz ^
        -keep_class 2 ^
        -kill 1000 ^
        -o LiDAR_Kantsi_dtm_1m.bil
lasheight -i Kantsi_shifted_again.laz ^
          -ground_points LiDAR_Kantsi_dtm_1m.bil ^
          -classify_below -0.4 7 ^
          -classify_between -0.4 1.0 8 ^
          -o Kantsi_cleaned.laz

Below the results we have achieved after “roughly” aligning the two point clouds with some new lasview tricks and then using the LiDAR elevations to classify the photogrammetry points into “low noise”, “potential ground”, and “all else”.

We thank the Estonian Land Board for providing¬†open data with a permissive¬†license. Special thanks also go to the organizers of the UAV Summer School in Tartu, Estonia and the European Regional Development Fund for funding this event. Especially fun was the fabulous excursion to the Emaj√Ķe-Suursoo Nature Reserve and through to Lake Peipus aboard, overboard and aboveboard the historical barge “J√Ķmmu”. If you look carefully you can also find the barge in the photogrammetry point cloud. The photogrammetry data used here was acquired during our lunch stop.

Fun aboard and overboard the historical barge “J√Ķmmu”.

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

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

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.