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]
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.
Side by side comparison of DTM and DSM generated with LAStools from photogrammetry points.
Side by side comparison of DTM and DSM generated with LAStools from photogrammetry points.
Side by side comparison of DTM and DSM generated with LAStools from photogrammetry points.
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.
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.
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).
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).
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.
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.
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.
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.
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.
Recently we attempted to do a small LiDAR survey by drone for a pet project of our CEO in our “code and surf camp” here in Samara, Costa Rica. But surveying is difficult when you are a novice and we ran into a trajectory issue. The dramatic “wobbles” were entirely our fault, but fortunately our mistakes also led to something useful: We found some LAS export bugs. Our laser scanner was a Velodyne HDL-32E integrated with a NovAtel INS into the Snoopy Series A HD made by LiDARUSA. The system was carried by a DJI Matrice 600 (M600) drone. We processed the trajectory with NovAtel Inertial Explorer (here we made the “wobbles” error) and finally exported the LAS and LAZ files with ScanLook PC (version 1.0.182) from LiDARUSA.
While we were investigating our “wobbles” (which clearly were our mistake) we also found five different LAS export bugs in ScanLook PC that seem to have started sometime after version 1.0.171 and will likely end with version 1.0.193. Below an illustration of a correct export from version 1.0.129 and a buggy export from version 1.0.182. In both instances you see the returns from one revolution of the Velodyne HDL-32E scanner head ordered by their GPS time stamps and colored to distinguish the 32 separate beams. In the buggy version, groups of around seven non-adjacent returns are given the same time stamp. This bug will only affect you, if correct GPS time stamps are important for your subsequent LiDAR processing or if your client explicitly asked for ASPRS specification compliant LAS files. We plan to publish another blog post detailing how to find this GPS time stamping bug (and the other four bugs we found).
First we generate a lasinfo report that includes a number of histograms for on-the-fly merged flight lines with lasinfo and then use the z coordinate histogram from the lasinfo report to set reasonable min/max values for the elevation color ramp of lasview:
The lasinfo report shows no information about the coordinate reference system. We found out experimentally that the horizontal coordinates seem to be EPSG code 2236 and that the vertical units are most likely be US survey feet. The warnings you will see in the lasinfo report have to do with the fact that the double-precision bounding box stored in the LAS header was populated with numbers that have many more decimal digits than the coordinates in the file, which only have millifeet resolution as all three scale factors are 0.001 (meaning coordinates have three decimal digits). The information which of the 32 lasers was collecting which point is stored in both the ‘user data’ and the ‘point source ID’ field which is evident from the histograms in the lasinfo report. We need to be careful not to override both fields in later processing.
Point cloud colored by elevation with blue being 25 feet and red being 75 feet.
Point cloud colored by elevation with black being 25 feet and white being 75 feet.
Point cloud colored by intensity suddenly shows the landing strip that has low reflectivity.
Next we use lasoverlap to check how well the LiDAR points from the flight out and the flight back align vertically. This tool computes the difference of the lowest points for each square foot covered by multiple flight lines. Differences of less than a quarter of a foot are both times mapped to white, differences of more than one foot (more than half a foot) are mapped to saturated red or blue depending on whether the difference is positive or negative in the first run (in the second run):
Difference of the lowest points for each square foot covered by multiple flight lines: differences of less than a quarter of a foot are white, differences of more than one foot are saturated red or blue.
Difference of the lowest points for each square foot covered by multiple flight lines: differences of less than a quarter of a foot are white, differences of more than half a foot are saturated red or blue.
We use a new feature of the LAStools GUI (as of version 180429) to closer inspect large red or blue areas. With lasmerge we clip out regions that looks suspect for closer examination with lasview. First we spatially index the flight lines to make this process faster. With the ‘-gui’ switch we start the tool in GUI mode with flight lines already loaded. Using the new PNG overlay roll-out on the left we add the ‘overlap_025_050_diff.png’ image from the quality folder created in the last step and clip out three areas.
The reader may inspect the areas 939500_889860.laz, 940400_889620.laz, and 940500_890180.laz with lasview using profile views via hot keys ‘x’ and switching back and forth between the points from different flight lines via hot keys ‘0’, ‘1’, ‘2’, ‘3’, … for individual and ‘a’ for all flight lines as we have done it in previous tutorials [1,2,3]. Using drop-lines or rise-lines via the pop-up menu gives you a sense of scale. Removing points with lastrack that are horizontally too far from the trajectory could be one strategy to use fewer outliers. But as our surfaces are expected to be “fluffy” (because we have a Velodyne LiDAR system), we accept these flight line differences and continue processing.
Hillshaded DTM with half foot resolution generated via average ground computation with LAStools.
In the first step we lastile the six flight lines into 250 by 250 feet tiles with 25 feet buffer while preserving flight line information. The flight line information will be stored in the “point source ID” field of each point and therefore override the beam ID that is currently stored there. But the beam ID is also stored in the “user data” field as the lasinfo report had told us. We set all classifications to zero and add information about the horizontal coordinate reference system EPSG code 2236 and the vertical units (US Survey Feet).
Next we use lasthin to classify the point whose elevation is closest to the 5th elevation percentile among all points falling into its cell with classification code 8. We run lasthin multiple times and each time increase the cell size from 1, 2, 4, 8 to 16 foot. We do this because we have requested the 5th elevation percentile to only be computed when there are at least 20 points in the cell. Percentiles are statistical measures and need a reasonable sample size to be stable. Because drone flights are very dense in the center and more sparse at the edges this increase in cell size assures that we have a good selection of points classified with classification code 8 across the entire survey area.
Then we let lasground_new run on only the points classified with classification code 8 (i.e. by ignoring the points still classified with code 0) which classifies them into ground (code 2) and non-ground (code 1).
The ground points we have computed form somewhat of a lower envelope of the “fluffy” points of a Velodyne scanner. With lasheight we now draw all the points near the ground – namely those from 0.1 foot below to 0.4 foot above the ground – into a new classification code 6 that we term “thick ground”. The ‘-do_not_store_in_user_data’ switch prevent the default behavior of lasheight from happening, which would override the beam ID information that it stored in the ‘user data’ field with approximate height value.
A few close-up shots of the resulting “thick ground” are shown in the picture gallery below.
Thick ground points in orange. Unclassified points in black.
Thick ground only.
Unclassified points only. The z distance between black points above and below the thick ground is just above half a feet.
Thick ground points in orange. Unclassified points in black.
Thick ground only.
Unclassified points only.
Thick ground points in orange. Unclassified points in black.
Thick ground only.
Unclassified points only. The z distance between black points above and below the thick ground is just above half a feet.
We then use lasgrid to average the (orange) thick ground points onto a regular grid with a cell spacing of half a foot. We do not grid the tile buffers by adding the ‘-use_tile_bb’ switch.
Finally we use blast2dem to merge all the averaged ground point grids into one file, interpolate across open areas without ground points, and compute the hillshaded DTM shown above. All command lines used are summarized in this text file.
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.
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’.
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:
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.
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:
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.
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”.
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
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.
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.
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.
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:
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.
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.
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.
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).
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.
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.
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.
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.
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.
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.
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.
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.
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).
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.
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.
