Calculating Grid Statistics with lidR

So far I have been working only with the LiDAR-derived DEM for calculating my topographic indices (TWI and TPI). For the Grid Statistics, I will be using the raw LiDAR point cloud (.las or in my case .laz files since they are compressed). There are many ways we can summarize this data and what I need to figure out is which are useful for mapping wetlands. Here is a list of grid statistics previously used for mapping wetlands or characterizing vegetation (note: z = elevation, i = intensity, std = standard deviation, max = maximum, and min = minimum):

Corcoran et al. (2015) used z-std, z-max, z-min and z-mean as well as i-std, i-max, i-min, and i-mean for their wetland classification. Intensity can been used as a correlation between wetness and soil moisture. Their results showed that z-std, z-max, i-mean, i-min and i-deviation ranked particularly high for variable importance.

In a study by Millard and Richardson (2013), they used standard deviation of all hits to approximate for height of vegetation; standard deviation of ground hits as a proxy for hummocky peatlands; the ratio between all hits and ground hits to indicate canopy density; ground hit intensity for flat surfaces/open water; and the density of all hits for their model. They found z-std, ratio of all hits to ground hits, and the density of all hits (d) to rank high in importance. They also concluded that intensity could not be used to distinguish between wetland classes.

In Hopkinson et al. study from 2006, these authors suggest that z-std provided an important metric for estimating vegetation height. In forest environments, there is a greater standard deviation due to distribution of returns from LiDAR (Millard and Richardson 2013). Millard and Richardson (2013) also found reduced standard deviation within bogs containing little/no tall vegetation.

The coefficient of variation (standard deviation divided by the mean) can be used as a measure of texture, which can reveal landscape structural patterns that are characteristic of different wetland types.

Calculating grid metrics from LiDAR data using lidR in R

Grid statistics are an area-based calculation for each pixel. I used the package lidR in R to calculate the following grid statistics: zmin, zmax, zmean, zstd, imin, imax, imean, istd. I chose to start with these simpler calculations to get a handle on how the package worked before calculating more complicated metrics such as the ratio of all hits to ground hits, coefficient of variation and density.

In order to run grid metrics using the lidR package, you must specify the resolution (in this case I used 10m), the data you are using (z, i, classification, return number, etc.), and the output file name (more on this later as it outputs as a data stack).

There are several important settings to understand when using the lidR package and working with a large las catalogue. My notes and R-code can be viewed here. In addition, some useful documents and online help pages can be found at the links below.

LiDAR basics:

https://www.neonscience.org/resources/learning-hub/tutorials/lidar-basics

LiDAR BC specs:

https://www2.gov.bc.ca/assets/gov/data/geographic/digital-imagery/geobc_lidar_specifications_v40.pdf

https://www.asprs.org/wp-content/uploads/2010/12/LAS_1_4_r13.pdf

LAS catalogue how to:

https://jean-romain.github.io/lidRbook/engine.html

https://rdrr.io/cran/lidR/man/catalog_retile.html

https://gis.stackexchange.com/questions/308843/calculating-standard-metrics-in-multiple-lidar-files-with-lidr-2-0

lidR package:

https://cran.r-project.org/web/packages/lidR/lidR.pdf

https://github.com/Jean-Romain/lidR/blob/master/vignettes/lidR-LAS-class.Rmd

https://jean-romain.github.io/lidRbook/engine2.html

https://cran.r-project.org/web/packages/lidR/vignettes/lidR-LAScatalog-class.html

https://cran.r-project.org/web/packages/lidR/vignettes/lidR-LAScatalog-engine.html#filter-points

https://jean-romain.github.io/lidRbook/aba.html

Overlapping scans:

https://gis.stackexchange.com/questions/270961/open-source-approach-to-classifying-and-removing-lidar-points-from-overlapping-s

Future package:

https://www.r-bloggers.com/2020/09/future-1-19-1-making-sure-proper-random-numbers-are-produced-in-parallel-processing/

Parallel computing:

https://rdrr.io/cran/lidR/man/lidR-parallelism.html

Trying to understand random numbers:

https://stackoverflow.com/questions/61781599/why-random-number-generation-using-parallel-processing-makes-different-results-e

Grid metrics: 

https://rdrr.io/cran/lidR/man/grid_metrics.html

https://www.rdocumentation.org/packages/lidR/versions/1.2.1/topics/grid_metrics

https://gis.stackexchange.com/questions/310834/filtering-lasfilter-on-huge-lidar-dataset/310847#310847

Naming each metric:

https://gis.stackexchange.com/questions/347064/keeping-name-of-metrics-as-band-names-when-using-grid-metrics-lidr

Working with rasters in R:

https://www.neonscience.org/resources/learning-hub/workshops/work-lidar-derived-rasters-r

https://www.neonscience.org/resources/learning-hub/workshops/hands-primer-working-big-data-r-introduction-hierarchical-data

References:

Corcoran, J., Knight, J., Pelletier, K., Rampi, L., & Wang, Y. (2015). The Effects of Point or Polygon Based Training Data on RandomForest Classification Accuracy of Wetlands. Remote Sensing, 7(4), 4002–4025. https://doi.org/10.3390/rs70404002.

Jean-Romain Roussel and David Auty (2021). Airborne LiDAR Data Manipulation and Visualization for Forestry Applications. R package version 3.1.0. https://cran.r-project.org/package=lidR

Millard, K., & Richardson, M. (2013). Wetland mapping with LiDAR derivatives, SAR polarimetric decompositions, and LiDAR–SAR fusion using a random forest classifier. Canadian Journal of Remote Sensing, 39(4), 290–307.

Roussel, J.R., Auty, D., Coops, N. C., Tompalski, P., Goodbody, T. R. H., Sánchez Meador, A., Bourdon, J.F., De Boissieu, F., Achim, A. (2020). lidR : An R package for analysis of Airborne Laser Scanning (ALS) data. Remote Sensing of Environment, 251 (August), 112061. doi:10.1016/j.rse.2020.112061.