Flying LiDAR in regions with frequent cloud cover presents a significant challenge. If flight plan constraints do not allow to stay below all of the clouds then some of them will be scanned from above. For denser clouds this often means that all of the laser’s energy gets reflected or absorbed by the cloud and no returns on the terrain are generated. Clouds of points that are all cloud points can spell trouble in subsequent processing steps like automated classification with lasground. This is especially true for large dense clouds that start to look like features.
Airborne LiDAR surveys are often carried out to create an improved Digital Terrain Model (DTM) with higher resolution than previous elevation products. Hence, usually there is already some lower resolution model that – as we will see in the following – can be used to robustly remove or mark all cloud points, at least those sufficiently high above this older ground approximation. We use data from the DREAM LiDAR Project in the Philippines who often acquire LiDAR in areas with a lot of cloud cover.
Our input are the 4 very short LiDAR strips in LAZ format shown above with lots of cloud returns and a coarse Aster DTM in ASC format at 50 m resolution. We will classify all LiDAR points far above the Aster DTM because they correspond to cloud returns. We need to be very conservative (because the Aster DTM is coarse and inaccurate) and only remove points that are really far above the Aster DTM whose ASC header is shown below.
ncols 2109 nrows 2162 xllcorner 512546.427859 yllcorner 1290961.682335 cellsize 50 NODATA_value -9999 -9999 -9999 -9999 -9999 -9999 [...]
First we convert the Aster DTM from the inefficient ASC format to the efficient LAZ format using lassort. Why lassort? Because that puts the rasters that are on-the-fly converted to points on a grid into a spatially coherent order. This will allow efficient area-of-interest (AOI) queries once we have LAXxed the file with lasindex.
>> lassort -i masbate.asc -o masbate.laz >> dir masbate* 10,993,886 masbate.asc 2,120,442 masbate.laz
Next we create a spatial indexing for ‘masbate.laz’ with lasindex. The default granularity of lasindex is 100 by 100 meter because it is set for airborne LiDAR. Given that there is only one point every 50 by 50 meters in the point cloud grid of the Aster DTM we increase the granularity to 1000 by 1000 meters.
>> lasindex -i masbate.laz -tile_size 1000 -append >> dir masbate* 10,993,886 masbate.asc 2,132,006 masbate.laz
The ‘masbate.laz’ file is a little bit bigger than before because the tiny ‘masbate.lax’ file (that can also be stored separately) was appended to the end of ‘masbate.laz’ via the ‘-append’ option. Spatial indexing is realized with an underlying quadtree that can be visualized with lasview by pressing ‘Q’ (or with the pop-up menu) and the points corresponding to each quadtree cell can be turned blue by hovering above a cell and pressing ‘q’.
Before we use the sorted points from the Aster DTM to classify, flag, or delete the LiDAR returns far above the ground with lasheight we do a visual sanity check with lasview using the GUI settings shown below.
>> lasview -i raw_strips\*.laz ^ -i masbate.laz ^ -gui
By triangulating only the Aster DTM points we can visually confirm that the LiDAR points on the terrain are also close to the Aster DTM whereas the cloud returns are far up with a clear seperation between. The pink “drop lines” of length 50 meters that are dangling off each point give us a sense of scale (enabled via the pop-up menu that appears with a right click). Note how many last returns get stuck in the clouds. Those would likely cause troubles in subsequent ground classification with lasground.
Finally we run lasheight on 4 cores to classify all points that are 150 meter or more above the Aster DTM as noise (class 7).
>> lasheight -i raw_strips\*.laz ^ -ground_points masbate.laz ^ -classify_above 150 7 ^ -odix _cloud -olaz ^ -cores 4
Below the result with the cloud returns in pink (i.e. the points classified as 7). This process scales well even for larger ground points files, because lasheight uses an area-of-index (AOI) query to load only those ground points from the spatially indexed ‘masbate.laz’ file, that fall into the (slightly enlarged) bounding box of the raw LiDAR strip that is being processed.