NRW Open LiDAR: Merging Points into Proper LAS Files

In the first part of this series we downloaded, compressed, and viewed some of the newly released open LiDAR data for the state of North Rhine-Westphalia. In the second part we look at how to merge the multiple point clouds provided back into single LAS or LAZ files that are as proper as possible. Follow along with a recent version of LAStools and a pair of DGM and DOM files for your area of interest. For downloading the LiDAR we suggest using the wget command line tool with option ‘-c’ that after interruption in transmission will restart where it left off.

In the first part of this series we downloaded the pair of DGM and DOM files for the City of Bonn. The DGM file and the DOM file are zipped archives that contain the points in 1km by 1km tiles stored as x, y, z coordinates with centimeter resolution. We had already converted these textual *.xyz files into binary *.laz files. We did this with the open source LASzip compressor that is distributed with LAStools as described in that blog post. We continue now with the assumption that you have converted all of the *.xyz files to *.laz files as described here.

Mapping from tile names of DGM and DOM archives to classification and return type of points.

The mapping from tile names in DGM and DOM archives to the classification and return type of points: lp = last return. fp = first return, ab,aw,ag = synthetic points

There are multiple tiles for each square kilometer as the LiDAR has been split into different files based on classification and return type. Furthermore there are also synthetic points that were created by the land survey department to replace LiDAR under bridges and along buildings for generating higher quality rasters. We want to combine all points of a square kilometer into a single LAZ tile as it is usually expected. Simply merging the multiple files per tile is not an option as this would result in loosing point classifications and return type information as well as in duplicating all single returns that are stored in more than one file. The folks at OpenNRW offer this helpful graphic to know what the acronyms above mean:

Illustration of how acronyms used in tile names correspond to point classification and type.

Illustration of how acronyms used in tile names correspond to point classification and type.

In the following we’ll be looking at the set of files corresponding to the UTM tile 32366 / 5622. We wanted an interesting area with large buildings, a bridge, and water. We were looking for a suitable area using the KML overlays generated in part one. The tile along the Rhine river selected in the picture below covers the old city hall, the opera house, and the “Kennedy Bridge” has a complete set of DGM and DOM files:

      3,501 dgm1l-ab_32366_5622_1_nw.laz
     16,061 dgm1l-ag_32366_5622_1_nw.laz
      3,269 dgm1l-aw_32366_5622_1_nw.laz
    497,008 dgm1l-brk_32366_5622_1_nw.laz
  7,667,715 dgm1l-lpb_32366_5622_1_nw.laz
 12,096,856 dgm1l-lpnb_32366_5622_1_nw.laz
     15,856 dgm1l-lpub_32366_5622_1_nw.laz

      3,269 dom1l-aw_32366_5622_1_nw.laz
 21,381,106 dom1l-fp_32366_5622_1_nw.laz
We find the name of the tiles that cover the "Kennedy Bridge" using the KML overlays generated in part one.

We find the name of the tile that covers the “Kennedy Bridge” using the KML overlays generated in part one.

We now assign classification codes and flags to the returns from the different files using las2las, merge them together with lasmerge, and recover single, first, and last return information with lasduplicate. We set classifications to bridge deck (17), ground (2), to unclassified (1), and to noise (7) for all returns in the files with the acronym ‘brk’ (= bridge points), the acronym ‘lpb’ (= last return ground), the acronym ‘lpnb’ (= last return non-ground), and the acronym ‘lpub’ (= last return under ground). with las2las and store the resulting files to a temporary folder.

las2las -i dgm1l-brk_32366_5622_1_nw.laz ^
        -set_classification 17 ^
        -odir temp -olaz

las2las -i dgm1l-lpb_32366_5622_1_nw.laz ^
        -set_classification 2 ^
        -odir temp -olaz

las2las -i dgm1l-lpnb_32366_5622_1_nw.laz ^
        -set_classification 1 ^
        -odir temp -olaz

las2las -i dgm1l-lpub_32366_5622_1_nw.laz ^
        -set_classification 7 ^
        -odir temp -olaz

Next we use the synthetic flag of the LAS format specification to flag any additional points that were added (no measured) by the survey department to generate better raster products. We set classifications to ground (2) and the synthetic flag for all points of the files with the acronym ‘ab’ (= additional ground) and the acronym ‘ag’ (= additional building footprint). We set classifications to water (9) and the synthetic flag for all points of the files with the acronym ‘aw’ (= additional water bodies). Files with acronym ‘aw’ appear both in the DGM and DOM archive. Obviously we need to keep only one copy.

las2las -i dgm1l-ab_32366_5622_1_nw.laz ^
        -set_classification 2 ^
        -set_synthetic_flag 1 ^
        -odir temp -olaz

las2las -i dgm1l-ag_32366_5622_1_nw.laz ^
        -set_classification 2 ^
        -set_synthetic_flag 1 ^
        -odir temp -olaz

las2las -i dgm1l-aw_32366_5622_1_nw.laz ^
        -set_classification 9 ^
        -set_synthetic_flag 1 ^
        -odir temp -olaz

Using lasmerge we merge all returns from files with acronyms ‘brk’ (= bridge points), ‘lpb’ (= last return ground),  ‘lpnb’ (= last return non-ground), and ‘lpub’ (= last return under ground) into a single file that will then contain all of the (classified) last returns for this tile.

lasmerge -i temp\dgm1l-brk_32366_5622_1_nw.laz ^
         -i temp\dgm1l-lpb_32366_5622_1_nw.laz ^
         -i temp\dgm1l-lpnb_32366_5622_1_nw.laz ^
         -i temp\dgm1l-lpub_32366_5622_1_nw.laz ^
         -o temp\dgm1l-lp_32366_5622_1_nw.laz

Next we run lasduplicate three times to recover which points are single returns and which points are the first and the last return of a pair of points generated by the same laser shot. First we run lasduplicate with option ‘-unique_xyz’ to remove any xyz duplicates from the last return file. We also mark all surviving returns as the second of two returns. Similarly, we remove any xyz duplicates from the first return file and mark all survivors as the first of two returns. Finally we run lasduplicate with option ‘-single_returns’ with the unique last and the unique first return files as ‘-merged’ input. If a return with the exact same xyz coordinates appears in both files only the first copy is kept and marked as a single return. In order to keep the flags and classifications from the last return file, the order in which the last and first return files are listed in the command line is important.

lasduplicate -i temp\dgm1l-lp_32366_5622_1_nw.laz ^
             -set_return_number 2 -set_number_of_returns 2 ^
             -unique_xyz ^
             -o temp\last_32366_5622_1_nw.laz

lasduplicate -i dom1l-fp_32366_5622_1_nw.laz ^
             -set_return_number 1 -set_number_of_returns 2 ^
             -unique_xyz ^
             -o temp\first_32366_5622_1_nw.laz

lasduplicate -i temp\last_32366_5622_1_nw.laz ^
             -i temp\first_32366_5622_1_nw.laz ^
             -merged ^
             -single_returns ^
             -o temp\all_32366_5622_1_nw.laz

We then add the synthetic points with another call to lasmerge to obtain a LAZ file containing all points of the tile correctly classified, flagged, and return-numbered.

lasmerge -i temp\dgm1l-ab_32366_5622_1_nw.laz ^
         -i temp\dgm1l-ag_32366_5622_1_nw.laz ^
         -i temp\dgm1l-aw_32366_5622_1_nw.laz ^
         -i temp\all_32366_5622_1_nw.laz ^
         -o temp\merged_32366_5622_1_nw.laz

Optional: For more efficient use of this file in subsequent processing – and especially to accelerate area-of-interest queries with lasindex – it is often of great advantage to reorder the points in a spatially coherent manner. A simple call to lassort will rearrange the points along a space-filling curve such as a Hilbert curve or a Z-order curve.

lassort -i temp\merged_32366_5622_1_nw.laz ^
        -o bonn_32366_5622_1_nw.laz

Note that we also renamed the file because a good name can be useful if you find that file again in two years from now. Let’s have a look at the result with lasview.

lasview -i bonn_32366_5622_1_nw.laz

In lasview you can press <c> repeatedly to switch through all available coloring modes until you see the yellow (single) / red (first) / last (blue) coloring of the returns. The recovered return types are especially evident under vegetation, in the presence of wires, and along building edges. Press <x> to select an area of interest and press <x> again to inspect it more closely. Press <i> while hovering above a point to show its coordinates, classification, and return type.

We did each processing in separate steps to illustrate the overall workflow. The above sequence of LAStools command line calls can be shortened by combining multiple processing steps into one operation. This is left as an exercise for the advanced LAStools user … (-;

Acknowledgement: The LiDAR data of OpenNRW comes with a very permissible license. It is called “Datenlizenz Deutschland – Namensnennung – Version 2.0” or “dl-de/by-2-0” and allows data and derivative sharing as well as commercial use. It only requires us to name the source. We need to cite the “Land NRW (2017)” with the year of the download in brackets and specify the Universal Resource Identification (URI) for both the DOM and the DGM. Done. So easy. Thank you, OpenNRW … (-:

LASmoons: Jesús García Sánchez

Jesús García Sánchez (recipient of three LASmoons)
Landscapes of Early Roman Colonization (LERC) project
Faculty of Archaeology, Leiden University, The Netherlands

Background:
Our project Landscapes of Early Roman Colonization (LERC) has been studying the hinterland of the Latin colony of Aesernia (Molise region, Italy) using several non-destructive techniques, chiefly artefactual survey, geophysics, and interpretation of aerial photographs. Nevertheless large areas of the territory are covered by the dense forests of the Matese mountains, a ridge belonging the Apennine chain, or covered by bushes due to the abandonment of the countryside. The project won’t be complete without integrating the marginal, remote and forested areas into our study of the Roman hinterland. Besides, it’s also relevant to discuss the feasibility of LiDAR data sets in the study of Mediterranean landscapes and its role within contemporary Landscape Archaeology.

some clever caption

LiDAR coverage in Molise region, Italy.

Goal:
+ to study in detail forested areas in the colonial hinterland of Aesernia.
+ to found the correct parameters of the classification algorithm to be able to locate possible archaeological structures or to document appropriately those we already known.
+ to document and create new visualization of hill-top fortified sites that belong to the indigenous population and are currently poorly studied due to inaccessibility and forest coverage (Monte San Paolo, Civitalla, Castelriporso, etc.)
+ to demonstrate the archaeological potential of LiDAR data in Italy and help other scholars to work with that kind of data, explaining basic information about data quality, where and how to acquire imagery and examples of application in archaeology. A paper entitled “Working with ALS – LiDAR data in Central South Italy. Tips and experiences”, will be presented in the International Mediterranean Survey Workshop by the end of February in Athens.

Civitella hillfort (Longano, IS) and its local context: ridges and forest belonging to the Materse mountains and the Appenines.

Data:
Recently the LERC project has acquired a large LiDAR dataset created by the Italian Geoportale Nazionale and the Minisstero dell’Ambiente e della Tutella del Territorio e del Mare. The data was produced originally to monitor land-slides and erosive risk.
The average point resolution is 1 meter.
+ The data sets were cropped originally in 1 sq km. tiles by the Geoportale Nazionale for distribution purposes.

LAStools processing:
1) data is provided in *.txt files thus the first step is to create appropriate LAS files to work with [txt2las]
2) combine areas of circa 16 sq km (still fewer than 20 million points to be processed in one piece with LAStools) in the surroundings of the colony of Aesernia and in the Matese mountains [lasmerge]
3) assign the correct projection to the data [lasmerge or las2las]
4) extract the care-earth with the best-fitting parameters [lasground or lasground_new]
5) create bare-earth terrain rasters as a first step to visualize and analyze the area [lasdem]

Pre-Processing Mobile Rail LiDAR with LAStools

The majority of LAStools users are processing airborne LiDAR. That should not surprise as airborne is by far the most common form of LiDAR in terms of square kilometers covered. The availability of LiDAR as “open data” is also pretty much restricted to airborne surveys, which are often tax-payer funded and then distributed freely to achieve maximum return of investment.

But folks are increasingly using our software to do some of the “heavy lifting” for mobile LiDAR, either mounted on a truck for scanning cities or on a train for capturing railroad infrastructure. The LiDAR collected for the cities of Budapest and Singapore, for example, was pre-processed by multi-core scripted LAStools when the scanning trucks returned with their daily trajectories worth of point clouds captured by a RIEGL VMX-450 mobile mapping system.

One customer who was recently scanning railroad infrastructure wanted to do automatic ground classification as a first step prior to further segmentation of the data. We were asked for advice because on such data the standard settings of lasground left too many patches of ground unclassified. Also the uniform tiling lastile generates by default is not a good way to break such data into manageable pieces given the drastically varying point densities in mobile scanning.

We obtained a 217 MB file in LAZ format with 40 million points corresponding to a 2.7 km stretch of railway track. We first run a quick lasindex (with the options for ‘mobile’) on the file that creates a spatial indexing LAX file with maximally 10 meter resolution. This not only allows faster area-of-interest queries but also gives us a more detailed preview than just the bounding box of where the LiDAR points actually are in the GUI of LAStools.

mobile_rail_lidar_01

Presence of LAX files results in actual extend of LiDAR being shown in GUI.

lasindex -i segment.laz -tile_size 10 -maximum -100

We then run lastile four times to create an adaptive tiling in which no tile has more than 6 million points. The first call creates the initial 1000 by 1000 meter tiles. The following three calls refine all those tiles that still have more than 6 million points first into 500 by 500 meter, then 250 by 250 meter, and finally 125 by 125 meter tiles in parallel on 4 cores. Note the ‘-refine_tiling’ option is used in the first call to lastile and the ‘-refine_tiles’ option in all subsequent calls.

lastile -i segment.laz ^
        -tile_size 1000 ^
        -buffer 10 -flag_as_withheld ^
        -refine_tiling 6000000 ^
        -odir tiles_raw -o rail.laz
lastile -i tiles_raw\*_1000.laz ^
        -flag_as_withheld ^
        -refine_tiles 6000000 ^
        -olaz ^
        -cores 4
lastile -i tiles_raw\*_500.laz ^
        -flag_as_withheld ^
        -refine_tiles 6000000 ^
        -olaz ^
        -cores 4
lastile -i tiles_raw\*_250.laz ^
        -flag_as_withheld ^
        -refine_tiles 6000000 ^
        -olaz ^
        -cores 4

The resulting tiles all have fewer than 6 million points but still have the initial 10 meter buffer that was specified by the first call to lastile. Two tiles were sufficiently small after the 1st call, three tiles after the 2nd call, eleven tiles after 3rd call, and three tiles after the 4th.

contents of tile shown in blue in adaptive tiling below

points of adaptive tile (high-lighted in blue below) colored by intensity

Adaptive tiling created with four calls to lastile.

Adaptive tiling created with four calls to lastile. Scale factors of 0.00025 (see mouse cursor) implies that point coordinates are stored with quarter millimeter resolution. Lowering them to 0.001 would result in better compression and lower I/O.

Noise in the data – especially low noise – can lead lasground into choosing the wrong points during ground classification by latching on to those low noise points. We first classify the noise points into a different class (7) using lasnoise so we can later ignore them. These particular settings were found by experimenting on a few tiles with different values (see the README file) until visual inspection showed that most low points had been classified as noise.

lasnoise -i tiles_raw\*.laz ^
         -step_xy 0.5 -step_z 0.1 ^
         -odir tiles_denoised -olaz ^
         -cores 4
noise points shown in violett

noise points shown in violett

The points classified as noise will not be considered as ground points during the next step. For this it matters little that lamp posts, wires, or vegetation are wrongly marked as noise now. We can always undo their noise classification once the ground points were classified. Important is that those pointed to by the mouse cursor, which are below the desired ground, are excluded from consideration during the ground classification step. Here those low points are not actually noise but returns generated wherever the laser was able to “peek” through an opening to a lower surface.

lasground -i tiles_denoised\*.laz ^
          -ignore_class 7 ^
          -step 1 -sub 3 -bulge 0.1 -spike 0.1 -offset 0.02 ^
          -odir tiles_ground -olaz ^
          -cores 4

For classification with lasground there are a number of options to play with  (see the README file) but the most important is the correct step size. It is terrain along the railway track bed that is supposed to get represented well. The usual step of 5 to 40 meter for lasground aim at the removal of vegetation and man-made structures from airborne LiDAR. They are not the right choice here. A step of 1 and the parameters shown above gives us the ground shown below.

Classification of terrain along railway track using lasground with '-step 1'

Classification of terrain along railway track bed using lasground with ‘-step 1’

The new ‘-flag_as_withheld’ option in lastile that flags each point in the buffer with the withheld flag is useful in case we want to remove all buffer points on-the-fly, for example, in order to create a DTM hillshade of 25 cm resolution for a visual quality check of the entire 2.7 km track using blast2dem from the BLAST extension of LAStools.

blast2dem -i tiles_ground\*.laz -merged ^
          -drop_withheld -keep_class 2 ^
          -hillshade -step 0.25 ^
          -o dtm_hillshaded.png
Small 600 x 600 pixel detail of hill-shaded 5663 x 9619 pixel DTM raster generated by blast2dem

Small 600 x 600 pixel detail of hill-shaded 5663 x 9619 pixel DTM raster generated by blast2dem.

LASmoons: Rachel Opitz

Rachel Opitz (recipient of three LASmoons)
Center for Virtualization and Applied Spatial Technologies
Department of Anthropology, University of South Florida, USA

Background:
In Spring 2017 Rachel Opitz will be teaching a course on Remote Sensing for Human Ecology and Archaeology at the University of South Florida. The aim of the course is to provide students with the practical skills and knowledge needed to work with contemporary remote sensing data. The course focuses on airborne laser scanning and hyper-spectral data and their application in Human Ecology and Archaeology. Through the course students will be introduced to a number of software packages commonly used to process and interpret these data, with an emphasis on free and/or open source tools.

Classification parameters and the resolution at which the DTM is interpolated both have a significant impact on our ability to recognize anthropogenic features in the landscape. Here we see a small quarry. More aggressive filtering and a coarser DTM resolution (left) makes it difficult to recognize that this is a quarry. Less aggressive filtering and a higher resolution (right) leaves some vegetation behind, but makes the edges of the quarry and some in-situ blocks clearly visible.

Goal:
The students will develop practical skills in applied remote sensing through hands-on exercises. Learning to assess, manage and process large data sets is essential. In particular, the students in the course will learn to:
+ Identify the set of techniques needed to solve a problem in applied remote sensing
+ Find public imagery and specify acquisitions
+ Assess data quality
+ Process airborne LiDAR data
+ Combine complementary remote sensing data sources
+ Create effective data visualizations
+ Analyze digital topographic and spectral data to answer questions in human ecology and archaeology

Data:
The course will include case studies that draw on a variety of publicly available data sets that will all be used in the exercises:
+ the PNOA data from Spain
+ data held by NOAA
+ data collected using NASA’s GLiHT platform

LAStools processing:
LAStools will be used throughout the course, as students learn to assess the quality of LiDAR data, classify raw LiDAR point clouds, create raster terrain and canopy models, and produce visualizations. The online tutorials and videos available via the company website and the over 50 hours of video on YouTube as well as the LAStools user forum will be used as resources during the course.

LASmoons: Alen Berta

Alen Berta (recipient of three LASmoons)
Department of Terrestrial Ecosystems and Landscape, Faculty of Forestry
University of Zagreb and Oikon Ltd Institute for Applied Ecology, CROATIA

Background:
After becoming the EU member state, Croatia is obliged to fulfill the obligation risen from the Kyoto protocol: National Inventory Report (NIR) of the Green House Gasses according to UNFCCC. One of the most important things during the creation of the NIR is to know how many forested areas there are and their wood stock and increment. This is needed to calculate the size of the existing carbon pool and its potential for sequestration. Since in Croatia, according to legislative, it is not mandatory to calculate the wood stock and yield of the degraded forest areas (shrubbery and thickets) during the creation of the usual forest management plans, this data is missing. So far, only a rough approximation of the wood stock and increment is used during the creation of NIR. However, these areas are expanding every year due to depopulation of the rural areas and the cessation of traditional farming.

very diverse stand structure of degraded forest areas (shrubbery and thickets)

Goal:
This study will focus on two things: (1) Developing regression models for biomass volume estimation in continental shrubberies and thickets based on airborne LiDAR data. To correlate LiDAR data with biomass volume, over 70 field plots with a radius of 12 meters have been established in more than 550 ha of the hilly and lowland shrubberies in Central Croatia and all trees and shrubberies above 1 cm Diameter at Breast Height (DBH) were recorded with information about tree species, DBH and height. Precise locations of the field plots are measured with survey GNNS and biomass is calculated with parameters from literature. For regression modeling, various statistics from the point clouds matching the field plots will be used (i.e. height percentiles, standard deviation, skewness, kurtosis, …). 2) Testing the developed models for different laser pulse densities to find out if there is a significant deviation from results if the LiDAR point cloud is thinner. This will be helpful for planning of the later scanning for the change detection (increment or degradation).

Data:
+
641 square km of discrete returns LiDAR data around the City of Zagreb, the capitol of Croatia (but since it is highly populated area, only the outskirts of the area will be used)
+ raw geo-referenced LAS files with up to 3 returns and an average last return point density of 1 pts/m².

LAStools processing:
1)
extract area of interest [lasclip or las2las]
2) create differently dense versions (for goal no. 2) [lasthin]
3) remove isolated noise points [lasnoise]
4) classify point clouds into ground and non-ground [lasground]
5) create a Digital Terrain Model (DTM) [las2dem]
6) compute height of points above the ground [lasheight]
7) classify point clouds into vegetation and other [lasclassify]
8) normalize height of the vegetation points [lasheight]
9) extract the areas of the field plots [lasclip]
10) compute various metrics for each plot [lascanopy]
11) convert LAZ to TXT for regression modeling in R [las2txt]

LASmoons: Stéphane Henriod

Stéphane Henriod (recipient of three LASmoons)
National Statistical Committee of the Kyrgyz Republic
Bishkek, Kyrgyzstan

This pilot study is part of the International Climate Initiative project called “Ecosystem based Adaptation to Climate change in the high mountainous regions of Central Asia” that is funded by the Federal Ministry for the Environment, Nature Conservation, Building and Nuclear Safety (BMU) of Germany and implemented by the Deutsche Gesellschaft für Internationale Zusammenarbeit (GIZ) GmbH in Kyrgyzstan, Tajikistan and Kazakhstan.

lasmoons_Stephane_Henriod_1

Background:
The ecosystems in high mountainous regions of Central Asia are characterized by a unique diversity of flora and fauna. In addition, they are the foundation of the livelihoods of the local population. Specific benefits include clean water, pasture, forest products, protection against floods and landslides, maintenance of soil fertility, and ecotourism. However, the consequences of climate change such as melting glaciers, changing river runoff regimes, and weather anomalies including sharp temperature fluctuations and non-typical precipitation result in negative impacts on these ecosystems. Coupled with unwise land use, these events damage fragile mountain ecosystems and reduce their regeneration ability undermining the local population’s livelihoods. Therefore, people living in rural areas and directly depending on natural resources must adapt to adverse impacts of climate change. This can be done through a set of measures, known in the world practice as ecosystem-based adaptation (EbA) approach. It promotes the sustainable use of natural resources to sustain and enhance the livelihood of the population depending on those resources.

lasmoons_Stephane_Henriod_2 Goal:
In two selected pilot regions of Kyrgyzstan and Tajikistan, planned measures will concentrate on climate-informed management of ecosystems in order to maintain their services for the rural population. EbA always starts with identifying the vulnerability of the local population. Besides analyzing the socio-economic situation of the local population, this includes (1) assessing the ecological conditions of the ecosystems in the watershed and the related ecosystem services people benefit from, (2) identifying potential disaster risks, and (3) analyzing glacier dynamics with its response to water runoff. As a baseline to achieve this and to get spatially explicit, remote sensing based techniques and mapping activities need to be utilized.

A first UAV (unmanned aerial vehicle) mission has taken place in the Darjomj watershed of the Bartang valley in July 2016. RGB-NIR images as well as a high-resolution Digital Surface Model have been produced that now need to be segmented and analysed in order to produce comprehensive information. The main processing that will take advantage of LAStools is the generation of a DTM from the DSM that will then be used for identifying risk areas (flood zones, landslides and avalanches, etc.). The results of this approach will ultimately be compared with lower-cost satellite images (RapidEye, Planet, Sentinel).

Data:
+ High-resolution RGB and NIR image (10 cm) from a SenseFly Ebee
+ High-resolution DSM (10 cm) from a SenseFly Ebee

LAStools processing:
1)
classify DSM points obtained via dense-matching photogrammetry into a SenseFly Ebee imagery into ground and non-ground points via processing pipelines as described here and here [lastile, lassort, lasnoise, lasground]
2) create a DTM [las2dem, lasgrid, blast2dem]
3) produce 3D visualisations to facilitate the communication around adaptation to climate change [lasview]
lasmoons_Stephane_Henriod_0

Creating a Better DTM from Photogrammetic Points by Avoiding Shadows

At INTERGEO 2016 in Hamburg, the guys from Aerowest GmbH shared with us a small photogrammetric point cloud from the city of Soest that had been generated with the SURE dense-matching software from nFrames. We want to test whether – using LAStools – we can generate a decent DTM from these points that are essentially a gridded DSM. If this interest you also see this, this, this, and this story.

soest_00_google_earth

Here you can download the four original tiles (tile1, tile2, tile3, tile4) that we are using in these experiments. We first re-tile them into smaller 100 meter by 100 meter tiles with a 20 meter buffer using lastile. We use option ‘-flag_as_withheld’ that flags all the points falling into the buffer using the withheld flag so they can easily be removed on-the-fly later with the ‘-drop_withheld’ filter (see the README for more on this). We also add the missing projection with ‘-epsg 32632’.

lastile -i original\*.laz ^
        -tile_size 100 -buffer 20 ^
        -flag_as_withheld ^
        -epsg 32632 ^
        -odir tiles_raw -o soest.laz

Below is a screenshot from one of the resulting 100 meter by 100 meter tiles (with 20 meter buffer) that we will be focusing on in the following experiments.

The tiles 'soest_437900_5713800.laz'

The tile ‘soest_437900_5713800.laz’ used in all experiments.

Next we use lassort to reorder the points ordered along a coherent space-filling curve as the existing scan-line order has the potential to cause our triangulation engine to slow down. We do this on 4 cores.

lassort -i tiles_raw\*.laz ^
        -odir tiles_sorted -olaz ^
        -cores 4

We then use lasthin to lower the number of points that we consider as ground points (see the README for more info on this tool). We do this because the original 5 cm spacing of the dense matching points is a bit excessive for generating a DTM with a resolution of, for example, 50 cm. Instead we only use the lowest point in each 20 cm by 20 cm cell as a candidate for becoming a ground point, which reduces the number of considered points by a factor of 16. We achieve this by classifying these lowest point to a unique classification code (here: 9) and later tell lasground to ignore all other classifications.

lasthin -i tiles_sorted\*.laz ^
        -step 0.2 -lowest -classify_as 9 ^
        -odir tiles_thinned -olaz ^
        -cores 4
Then we run lasground on 4 cores to classify the ground points with options ‘-step 20’, ‘-bulge 0.5’, ‘-spike 0.5’ and ‘-fine’ that we selected after two trials on a single tile. There are several other options in lasground to play with that may achieve better results on other data sets (see README file for more on this). The ‘-ignore_class 0’ option instructs lasground to ignore all points that are not classified so that only those points that lasthin was classifying as 9 in the previous step are considered.
lasground -i tiles_thinned\*.laz ^
          -step 20 -bulge 0.5 -spike 0.5 -fine ^
          -ignore_class 0 ^
          -odir tiles_ground -olaz ^
          -cores 4
However, when we scrutinize the resulting ground classification notice that there are bumps in the corresponding ground TIN that seem to correspond to areas where the original imagery has dark shadows that in turn seem to significantly affect the geometric accuracy of the points generated by the dense-matching algorithm.
Looking a the bump from below we identify the RGB colors of the points have that form the bump that seem to be of much lower accuracy than the rest. This is an effect that we have noticed in the past for data generated with other dense-matching software and maybe an approach similar to the one we take here could also help with this low noise. It seems that points that are generated from shadowed areas in the input images can have a lot lower accuracy than those from well-lit areas. We use this correlation between RGB color and geometric accuracy to simply exclude all points whose RGB colors indicate that they might be from shadow areas from becoming ground points.
The RGB colors of low-accuracy points are mostly from very dark shadowed areas.

The RGB colors of low-accuracy points are mostly from very dark shadowed areas.

We use las2las with the relatively new ‘-filtered_transform’ option to reclassify all points whose RGB color is close to zero to yet classification code 7 (see README file for more on this). We do this for all points whose red value is between 0 and 30, whose green value is between 0 and 35, and whose blue value is between 0 and 40. Because the RGB values were stored with 16 bits in these files we have to multiply these values with 256 when applying the filter.
las2las -i tiles_thinned\*.laz ^
        -keep_RGB_red 0 7680 ^
        -keep_RGB_green 0 8960 ^
        -keep_RGB_blue 0 10240 ^
        -filtered_transform ^
        -set_classification 7 ^
        -odir tiles_rgb_filtered -olaz ^
        -cores 4
Below you see all points that will no longer be considered because their classification was set to 7 by the command above.
Points whose RGB values indicate they might lie in the shadows.

Points whose RGB values indicate they might lie in the shadows.

We then re-run lasground with the same options ‘-step 20’, ‘-bulge 0.5’, ‘-spike 0.5’ and ‘-fine’ as before but now we ignore all points that are still have classification 0 because they were not classified as 9 by lasthin earlier and we also ignore all points that have been assigned classification 7 by las2las in the previous step.
lasground -i tiles_thinned\*.laz ^
          -step 20 -bulge 0.5 -spike 0.5 -fine ^
          -ignore_class 0 7 ^
          -odir tiles_ground -olaz ^
          -cores 4
The situation has significantly improved for the bumb we saw before as you can see in the images below.

Finally we create a DTM with blast2dem (see README) and a DSM with lasgrid (see README) by merging all points into one file but dropping the buffer points that were marked as withheld by the initial run of lastile (see README).

blast2dem -i tiles_ground\*.laz -merged ^
          -drop_withheld -keep_class 2 ^
          -step 0.5 ^
          -o dtm.bil

lasgrid -i tiles_ground\*.laz -merged ^
        -drop_withheld ^
        -step 0.5 -average ^
        -o dsm.bil
 As our venerable lasview (see README) can directly read BIL rasters as points (just like all the other LAStools), so we can compare the DTM and the DTM by loading them as two flight lines (‘-faf’) and then switch back and forth between the two by pressing ‘0’ and ‘1’ in the viewer.
lasview -i dtm.bil dsm.bil -faf

Above you see the final DTM and the original DSM. So yes, LAStools can definitely create a DTM from point clouds that are the result of dense-matching photogrammetry. We used the correlation between shadowed areas in the source image and geometric errors to remove those points from consideration for ground points that are likely to be too low and result in bumps. For dense-matching algorithms that also output an uncertainty value for each point there is the potential for even better results as our range of eliminated RGB colors may not cover all geometrically uncertain points while at the same time may be too conservative and also remove correct ground points.