We're really interesting in comparing the forest structure as measured by LiDAR to field measurements by ecologists. Effectively, can we take a slice through the point cloud and extract a vertical column. Our current plots record information across the forest strata: 0-0.5m, 0.5-1m. 1-2m, 2-5m, 5-10m, 10-15m, 15-20m, 20-35m, and >35m.
Each of the Pennsylvania LiDAR tiles covers is four kilometers to a side and thus the file size is usually about 350mb per tile. We have about 13,500 of them in the state, so we're only working with one for this test and we'll scale up later. It was relatively straightforward to load a single LAS file into R using lidR.
lidar = readLAS("50001990PAN.las")
Next, I generated a digital terrain model (DTM):
dtm = grid_terrain(lidar, res = 10, method = "knnidw", k = 6L)
This was a pretty quick proof-of-concept so the parameters for interpolating the grid are just placeholders. However, in later versions of this, we'll be able to use the DEM created by Pennsylvania TopoGeo as part of the original LiDAR collection.
The next step was to create a normalized point cloud that removes the ground contours and leaves us with just the canopy heights.
lidar_norm = lidar - dtm
From here, the next step was to calculate point density in each of the strata layers we measured in the field. First, we pulled out a slice of the point cloud (here's the zero to half meter extract:
ht0_0.5m = lidar_norm %>% lasfilter(Z>0.33,Z<1.64)
Then we calculate the point density:
den0m_0.5m <- grid_density(ht0_0.5m, res=5)
Here's a plot of the density for eight strata l