Category Archives: photogrammetry

LASmoons: Olumese Efeovbokhan

Posted on by

Olumese Efeovbokhan (recipient of another three LASmoons)
Geosciences, School of Geography
University of Nottingham, UK

Background:
Hydrological models require various input data for flood vulnerability mapping. An important input data for flood vulnerability mapping is the DTM over which flow is being routed. DTMs are generated using cartography, ground surveying, digital aerial photogrammetry, interferometric SAR (InSAR), LiDAR amongst other means. The accuracy of high resolution DTMs minimize errors that may emanate from input data when conducting hydrological modelling, especially in small built-up catchment areas. This research involves the application of digital aerial photogrammetry to generate point clouds which can subsequently be utilized for flood vulnerability mapping.

photogrammetry point cloud
photogrammetic DSM

Goal:
To consolidate on previous gains in using LAStools to generate DTMs required for flood vulnerability mapping. The suitability of these DTMs will be subsequently validated for flood vulnerability analysis. These results will be compared with other DTMs in order to determine the uncertainty associated with the use of such DTMs for flood vulnerability mapping.

photogrammetry point cloud
photogrammetric DSM

Data:
+ high-resolution photogrammetry point cloud and DSM for Lagos Island, Ikorodu and Ajah Nigeria
– – – imagery obtained with an Ebee Sensefly drone flight
– – – photogrammetry point cloud generated with Photoscan by AgiSoft 
+ rainfall data
+ classified LiDAR point cloud with a resolution of 1 pulse per square meter obtained for the study area from the Lagos State Government

Photogrammetry point cloud
Photogrammetric DSM

LAStools processing:
1) tile large photogrammetry point cloud into tiles with buffer [lastile]
2) mark set of points whose z coordinate is a certain percentile of that of their neighbors [lasthin]
3) remove isolated low points from the set of marked points [lasnoise]
4) classify marked points into ground and non-ground [lasground]
5) pull in points close above and below the ground [lasheight]
6) create Digital Terrain Model (DTM) from ground points [las2dem]
7) merge and hillshade individual raster DTMs [blast2dem]

LASmoons: Martin Romain

Posted on by

Martin Romain (recipient of three LASmoons)
Marshall Islands Conservation Society
Majuro, Republic of the MARSHALL ISLANDS

Background:
As a low-lying coastal nation, the Republic of the Marshall Islands (RMI) is at the forefront of exposure to climate change impacts. RMI has a strong dependence on natural resources and biodiversity not only for food and income but also for culture and livelihood. However, these resources are threatened by rising sea levels and associated coastal hazards (king tides, storm surges, wave run-up, saltwater intrusion, erosion). This project aims at addressing the lack of technical capacity and available data to implement effective risk reduction and adaptation measures, with a particular focus on inundation mapping and local evacuation planning in population centers.

DCIM100MEDIADJI_0507.JPG

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

Goal:
This project intends to use LAStools to generate a DEM of the inhabited sections of 3 remote atolls (Aur, Ebon, Likiep) and 1 island (Mejit). The resulting DEM will be used to produce an inundation exposure model (and map) under variable sea level rise projections for each site. The ultimate goal is to integrate the results into each site’s disaster risk reduction strategy (long-term outcome) and present it through community consultations in schools, community centers, and council houses.

Data:
+ Aerial imagery of 11.5 square kilometers of land (6.3% of total national landmass) using DJI Matrice 200 V2 & DJI Zenmuse X5S with a minimum overlap of 75/75 and maximum altitude of 120m.

LAStools processing:
1) tile large point cloud into tiles with buffer [lastile]
2) remove noise points [lasthin, lasnoise]
3) classify points into ground and non-ground [lasground]
4) create Digital Terrain Models and Digital Surface Models [lasthin, las2dem]

Potential LAStools pipelines:
1) Removing Excessive Low Noise from Dense-Matching Point Clouds
2) Digital Pothole Removal: Clean Road Surface from Noisy Pix4D Point Cloud
3) Creating DTMs from dense-matched points of UAV imagery from SenseFly’s eBee

LASmoons: Volga Lipwoni

Posted on by

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

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

lasmoons_Volga_Lipwoni

Typical point cloud derived with SfM software from UAV imagery.

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

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

LAStools processing:
1) tile large point cloud into tiles with buffer [lastile]
2) remove noise points [lasthin, lasnoise]
3) classify points into ground and non-ground [lasground]
4) create Digital Terrain Modelsand Digital Surface Models [lasthin, las2dem]
5) produce Canopy Height Models for computing tree heights [lasheight, las2dem]

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

LASmoons: Gabriele Garnero

Posted on by

Gabriele Garnero (recipient of three LASmoons)
Interuniversity Department of Regional and Urban Studies and Planning
Politecnico e Università degli Studi, Torino
ITALY

Background:
Last spring, the LARTU research group produced a laser scanner survey of the Abbey of Sant’Andrea in Vercelli, on the occasion of the VIII centenary of the dedication (1219). The database produced with a topographic tool that integrates the potential of a total station with laser scanner and photogrammetric sensors (Trimble SX 10), has been used to produce representations that can be consulted in interactive mode, navigating within the point clouds and producing a consultation platform that can also be accessed by non-specialist users such as art historians or archaeologists.

lasmoons_gabriele_garnero_0

Goal:
The LAStools software will be used to improve both the point cloud produced by eliminating the remaining noises, and check other ways of publishing the data, so as to make it usable from outside, to the community of researchers.

Data:
+ laser scanner and photogrammetric acquisitions of the interior of the building (150 millions of points)
+ laser scanner and photogrammetric acquisitions of the outside of the building (210 millions of points)
+ drone-based shooting of outdoor areas processed with Pix4D (23 millions of points)

LAStools processing:
1) tile large point cloud into tiles with buffer [lastile]
2) mark set of points whose z coordinate is a certain percentile of that of their neighbors [lasthin]
3) remove isolated low points from the set of marked points [lasnoise]
4) classify marked points into ground and non-ground [lasground]
5) creates a LiDAR portal for 3D visualization of LAS files [laspublish]

LASmoons: Olumese Efeovbokhan

Posted on by

Olumese Efeovbokhan (recipient of three LASmoons)
Geosciences, School of Geography
University of Nottingham, UK

Background:
One of the vital requirements to successfully drive and justify favorable flood risk management policies is the availability of reliable data for hydrological modelling. Unfortunately, this poses a big challenge in data-sparse regions and has resulted in uncoordinated and ineffective flood risk management policies with some areas left at the mercy of the floods they are exposed to. This research is focused on the ability to successfully generate data required for hydrological modelling using affordable and easy-to-replicate methods. The research will utilize unmanned aerial vehicles (UAVs) for the generation of bare earth models (DTMs) from photogrammetry points, which will be subsequently used for flood vulnerability mapping.

Photogrammetry point cloud of Tafawa Balewa Square in Lagos Island, Nigeria

Goal:
Generate a bare earth model using a combination of Agisoft Photoscan and LAStools and then validate its suitability for hydrological modelling. Should the generated model prove to be suitable we will use it to conduct flood sensitivity analysis and inundation modelling in other data-sparse regions using high resolution bare earth models generated the same way.

Data:
+ high-resolution photogrammetry point cloud for a portion of the study area
– – – imagery obtained with an Ebee Sensefly drone flight
– – – photogrammetry point cloud generated with Photoscan by AgiSoft 
+ classified LiDAR point cloud with a resolution of 1 pulse per square meter obtained for the study area from the Lagos State Government

LAStools processing:
1) tile large photogrammetry point cloud into tiles with buffer [lastile]
2) mark set of points whose z coordinate is a certain percentile of that of their neighbors [lasthin]
3) remove isolated low points from the set of marked points [lasnoise]
4) classify marked points into ground and non-ground [lasground]
5) pull in points close above and below the ground [lasheight]
6) create Digital Terrain Model (DTM) from ground points [las2dem]
7) merge and hillshade individual raster DTMs [blast2dem]

Clean DTM from Agisoft Photogrammetric Points of Urban Scene

Posted on by

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

Posted on by

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.

The miss-alignment between sparse LiDAR and dense photogrammetry points is clearly visible in open terrain areas.
The miss-alignment error between sparse LiDAR and dense photogrammetry points is easily measured along the roof of the building.
We can measure a displacement of around 2 meters vertically and 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.

Press ‘x’ once, select an area-of-interest around the house, and press ‘x’ again.
Use pop-up menu via right-click to select “reclassify points as” some unique and non-existing classification such as “road surface”.
Draw a polygon on a flat part of the roof.
Press ‘r’ to register the edit.
Use pop-up menu via right-click to change to “color by classification”.
Many wrong points. Use pop-up menu to render only “road surface”.
So many wrong points. Use pop-up menu to select “reclassify points as never classified”.
Draw a polygon around the wrong points. above.
Press ‘r’ to register the edit.
Press ‘e’ to continue editing and draw a polygon around the wrong points below.
Press ‘r’ to register the edit. Repeat until only points forming a planar area remain.
Press ‘p’ to compute a plane through these points.
Press ‘a’ to turn on all points.
Press ‘1’ to turn on the LiDAR points only (only works if lasview was run with option ‘-faf’) and press ‘t’ to triangulate them. Hover over a point and press ‘i’.
Hover over another point and press ‘i’ again. The distance from the plane to the point is displayed in the black control window.
Hover over another point and press ‘i’ again. The distances are similar.
Use pop-up menu to select “reclassify points as keypoint” and draw a polygon around the flat part of the triangulated surface.
Some wrong points again. Use pop-up menu to select “reclassify points as never classified”.
Draw a polygon around the wrong points and press ‘r’ to register the edit.
Once only the LiDAR points from the same roof that also form a plane are left press ‘l’.
A list of distances and an average distance / translation vector are output in the control window. We use this to shift the photogrammetry points closer to the LiDAR. 
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:

Comparing a wireframe TIN of the LiDAR with the photogrammetry points shows another misalignment.
Using the same technique as earlier we classify a subset of points from the building facade to a unique class and fit a plane.
By picking various different LiDAR points shows the perpendicular displacement vectors to the plane.
By picking various different LiDAR points shows the perpendicular displacement vectors to the plane.
By picking various different LiDAR points shows the perpendicular displacement vectors to the plane.
By picking various different LiDAR points shows the perpendicular displacement vectors to the plane.

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.

Original alignment.
After manual alignment with two interactively determined translation 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”.

Photogrammetry points colored by RGB.
Photogrammetry points colored by elevation.
Photogrammetry points colored by classification.
Points 40 centimeter or more below the LiDAR DTM classified as low noise.
Same as before but colored by RGB.
Points between 40 centimeter below and 100 centimeter above the LiDAR DTM classified as potential ground.
Points higher than 100 centimeter above the LiDAR DTM.

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

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

Posted on by

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.

all points.
after lasthin: all 70th elevation percentile points per 20 cm by 20 cm cell.
same. all 70th elevation percentile points per 20 cm by 20 cm cell.
after lasnoise with isolated points in blue.
after lasnoise without isolated points.
after lasground.
after lasheight: using surface defined by ground points to classify low noise.
low noise only.
points without low noise.
all points.
temporary ground surface for classifying low noise.
low noise points.
points without low noise.

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.

final DSM.
final DTM.

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

Posted on by

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.

 

Clumps of low noise points typical for photogrammetry point clouds.
The thinned 20th elevation percentile points.
The thinned 20th elevation percentile points classified with code 8.
In yellow the points classified to code 12 by lasnoise.
The ground points used to create a temporary surface.
The pink points below the temporary ground surface are noise.
The points classified as noise,
The clean photogrammetric point clouds.

 

Clumps of low noise points typical for photogrammetry point clouds.
The thinned 20th elevation percentile points.
The thinned 20th elevation percentile points classified with code 8.
In yellow the points classified to code 12 by lasnoise.
The ground points used to create a temporary surface.
The pink points below the temporary ground surface are noise.
The points classified as noise,
The clean photogrammetric point clouds.

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.

 

final 25 cm Digital Surface Model
final 25 cm Digital Terrain Model

Removing Excessive Low Noise from Dense-Matching Point Clouds

Posted on by

Point clouds produced with dense-matching by photogrammetry software such as SURE, Pix4D, or Photoscan can include a fair amount of the kind of “low noise” as seen below. Low noise causes trouble when attempting to construct a Digital Terrain Model (DTM) from the points as common algorithm for classifying points into ground and non-ground points – such as lasground – tend to “latch onto” those low points, thereby producing a poor representation of the terrain. This blog post describes one possible LAStools workflow for eliminating excessive low noise. It was developed after a question in the LAStools user forum by LASmoons holder Muriel Lavy who was able to share her noisy data with us. See this, this, this, this, this, and this blog post for further reading on this topic.

Lots of low noise in tile ‘panca_388500_4963000.laz’.
Cross section of noisy area of tile ‘panca_388500_4963000.laz’.
Lots of very low noise in tile ‘panca_388500_4963000.laz’.
Cross section of very noisy area of tile ‘panca_388500_4963000.laz’.

Here you can download the dense matching point cloud that we are using in the following work flow:

We leave the usual inspection of the content with lasinfo, lasview, and lasvalidate that we always recommend on newly obtained data as an exercise to the reader. Note that a check for proper alignment of flightlines with lasoverlap that we consider mandatory for LiDAR data is not applicable for dense-matching points.

With lastile we turn the original file with 87,261,083 points into many smaller 500 by 500 meter tiles for efficient multi-core processing. Each tile is given a 25 meter buffer to avoid edge artifacts. The buffer points are marked as withheld for easier on-the-fly removal. We add a (terser) description of the WGS84 UTM zone 32N to each tile via the corresponding EPSG code 32632:
lastile -i muriel\20161127_Pancalieri_UTM.laz ^
        -tile_size 500 -buffer 25 -flag_as_withheld ^
        -epsg 32632 ^
        -odir muriel\tiles_raw -o panca.laz
Because dense-matching points often have a poor point order in the files they get delivered in we use lassort to rearrange them into a space-filling curve order as this will speed up most following processing steps:
lassort -i muriel\tiles_raw\panca*.laz ^
        -odir muriel\tiles_sorted -olaz ^
        -cores 7
We then run lasthin to reclassify the highest point of every 2.5 by 2.5 meter grid cell with classification code 8. As the spacing of the dense-matched points is around 40 cm in both x and y, around 40 points will fall into each such grid cell from which the highest is then classified as 8:
lasthin -i muriel\tiles_sorted\panca*.laz ^
        -step 2.5 ^
        -highest -classify_as 8 ^
        -odir muriel\tiles_thinned -olaz ^
        -cores 7
Considering only those points classified as 8 in the last step we then run lasnoise to find points that are highly isolated in wide and flat neighborhoods that are then reclassified as 7. See the README file of lasnoise for a detailed explanation of the different parameters:
lasnoise -i muriel\tiles_thinned\panca*.laz ^
         -ignore_class 0 ^
         -step_xy 5 -step_z 0.1 -isolated 4 ^
         -classify_as 7 ^
         -odir muriel\tiles_isolated -olaz ^
         -cores 7
Now we run a temporary ground classification of only (!!!) on those points that are still classified as 8 using the default parameters of lasground. Hence we only use the points that were the highest points on the 2.5 by 2.5 meter grid and that were not classified as noise in the previous step. See the README file of lasground for a detailed explanation of the different parameters:
lasground -i muriel\tiles_isolated\panca*.laz ^
          -city -ultra_fine -ignore_class 0 7 ^
          -odir muriel\tiles_temp_ground -olaz ^
          -cores 7
The result of this temporary ground filtering is then merely used to mark all points that are 0.5 meter below the triangulated TIN of these temporary ground points with classification code 12 using lasheight. See the README file of lasheight for a detailed explanation of the different parameters:
lasheight -i muriel\tiles_temp_ground\panca*.laz ^
          -do_not_store_in_user_data ^
          -classify_below -0.5 12 ^
          -odir muriel\tiles_temp_denoised -olaz ^
          -cores 7
In the resulting tiles the low noise (but also many points above the ground) are now marked and in a final step we produce properly classified denoised tiles by re-mapping the temporary classification codes to conventions that are more consistent with the ASPRS LAS specification using las2las:
las2las -i muriel\tiles_temp_denoised\panca*.laz ^
        -change_classification_from_to 1 0 ^
        -change_classification_from_to 2 0 ^
        -change_classification_from_to 7 0 ^
        -change_classification_from_to 12 7 ^
        -odir muriel\tiles_denoised -olaz ^
        -cores 7
Let us visually check what each of the above steps has produced by zooming in on a 300 meter by 100 meter strip of points with the bounding box (388500,4963125) to (388800,4963225) in tile ‘panca_388500_4963000.laz’:
Raw colored point cloud.
Highest point per 2.5 by 2.5 meter grid cell classified as 8.
Marking isolated points as noise with classification 7 (pink).
Point surviving thinning and noise removal are input for temporary ground classification.
Temporary ground points used for low noise removal.
All raw points falling 0.5 meter under TIN of temporary ground points become low noise (pink).
Points that are not low noise are used as input for final ground classification.
Final ground classification.
Triangulation of final ground points.
The final classification of all points that are not already classified as noise (7) into ground (2) or non-ground (1) was done with a final run of lasground. See the README file of lasground for a detailed explanation of the different parameters:
lasground -i muriel\tiles_denoised\panca*.laz ^
          -ignore_class 7 ^
          -city -ultra_fine ^
          -odir muriel\tiles_ground -olaz ^
          -cores 7
Then we create a seamless hill-shaded DTM tiles by triangulating all the points classified as ground into a temporary TIN (including those in the 25 meter buffer) and then rasterizing only the inner 500 meter by 500 meter of each tile with option ‘-use_tile_bb’ of las2dem. For more details on the importance of buffers in tile-based processing see this blog post here.
las2dem -i muriel\tiles_ground\panca*.laz ^
        -keep_class 2 ^
        -step 1 -hillshade ^
        -use_tile_bb ^
        -odir muriel\tiles_dtm -opng ^
        -cores 7

And here the original DSM side-by-side with resulting DTM after low noise removal. One dense forested area near the center of the data was not entirely removed due to the lack of ground points in this area. Integrating external ground points or manual editing with lasview are two possible way to rectify these few remaining errors …

Original DSM.
Classified DTM.