In adam-erickson/gapfraction: Functions for processing airborne LiDAR scans of forests

"The science of today is the technology of tomorrow." Edward Teller

The gapfraction package for R is designed for processing airborne laser scanning (ALS) light-detection-and-ranging (LiDAR) data of forests. The package was produced for a chapter of my doctoral thesis at University of British Columbia^[Erickson, A. (2017) A hybrid modeling approach to simulating past-century understory solar irradiation in Alberta, Canada. University of British Columbia.]. The package is designed to be used with LiDAR data pre-processed with LAStools^[LAStools: http://rapidlasso.com/lastools/], USDA Fusion^[USDA Fusion: http://forsys.cfr.washington.edu/fusion/fusionlatest.html], or the new lidR package from Jean-Romain Roussel^[lidR: https://github.com/Jean-Romain/lidR]. The main input to functions are LAS or LAZ format height-normalized point clouds, typically LiDAR plots corresponding to field plots. The functions are designed to accept either paths to height-normalized LAS/LAZ files or data.frame/data.table objects created with lidR or rLiDAR from Carlos Silva^[rLiDAR: https://cran.r-project.org/web/packages/rLiDAR/index.html] packages. Functions automatically detect whether a path or object is input for the las parameter.

Functions included

The gapfraction package implements a new barycentric interpolation spike-free canopy height model (CHM) algorithm based on Khosravipour et al. (2013)^[Khosravipour et al. (2013) Development of an algorithm to generate a LiDAR pit-free canopy height model. http://www.riegl.com/uploads/tx_pxpriegldownloads/khosravipour_SilviLaser2013.pdf] and (2016)^[Khosravipour et al. (2016) Generating spike-free digital surface models using LiDAR raw point clouds: A new approach for forestry applications. http://www.sciencedirect.com/science/article/pii/S0303243416300873], in addition to two new LiDAR metrics of canopy gap fraction ($P_o$), several individual tree crown (ITC) detection algorithms, canopy distance and direction metrics, effective leaf area index ($L_e$) and apparent clumping index ($ACI$) estimation methods. Four mathematical fisheye (hemispherical) lens models are implemented: equi-angular, equi-distant, stereographic, and orthographic. Finally, a new radial.grid function is included for polar plots, while a modified sun path plotting function is also provided. An alphabetical list of functions contained in the gapfraction package is provided below.

• chm - Simple canopy height model
• chm.sf - Barycentric interpolation spike-free canopy height model
• dd.canopy - Euclidean distance and direction to nearest canopy pixel from plot center
• dd.crown - Euclidean distance and direction to nearest tree crown from plot centers
• itc.mw - Variable-window individual tree crown detection
• itc.mw.h - Hierarchical variable-window individual tree crown detection
• itc.wat - Watershed segmentation individual tree crown detection
• itc.wat.h - Hierarchical watershed segmentation individual tree crown detection
• Le.bl - Ground-to-total-return ratio with a spherical leaf angle distribution
• Le.n - Contact frequency and vertical canopy cover-based effective LAI
• radial.grid.hemi - Modified radial.grid function supporting hemispherical lens geometries
• P.hv - Hemipsherical Voronoi canopy gap fraction
• P.hv.par - Parallel hemispherical Voronoi canopy gap fraction with SOCKS
• P.pdn - Point-density-normalized canopy gap fraction, effective LAI (Le), and ACI
• sun.path - Modified solar position plots of Thomas Steiner
• vcc.aci - Above-height cover index of vertical canopy cover
• vcc.bl - Beer-Lambert-Law-modified intensity-return ratio of vertical canopy cover
• vcc.cv - 2-D Cartesian Voronoi vertical canopy cover
• vcc.fci - First-echo cover index of vertical canopy cover
• vcc.fr - Canopy-to-first-return ration of vertical canopy cover
• vcc.ir - Intesity-return ratio of vertical canopy cover
• vcc.p - Canopy-to-total-pixel ratio of vertical canopy cover
• vcc.r - Canopy-to-total-return ratio of vertical canopy cover
• vcc.sci - Solberg's cover index of vertical canopy cover

Canopy light transmission equations implemented

Point-density normalized gap fraction:

$$\begin{align} P_\text{pdn} &= \sum_{i=1}^{n} \sum_{j=1}^{n} \left(\frac{{n_\text{first return}}{i,j}} {\rho\text{first return} / {{A_\text{sector}}{i,j}}} \times \frac{A{\text{sector}{i,j}}}{A\text{hemisphere}}\right) \end{align}$$

Hemispherical Voronoi gap fraction above a canopy height threshold with an equidistant lens projection:

$$\begin{align} x_\text{equidistant} &= \left( \frac{x}{\sqrt{x^2+y^2}} \right) \times atan2 \left( \sqrt{x^2+y^2},z \right) \ y_\text{equidistant} &= \left( \frac{y}{\sqrt{x^2+y^2}} \right) \times atan2 \left( \sqrt{x^2+y^2},z \right) \[4pt] p &= \left( x_\text{equidistant}, y_\text{equidistant} \right) \[4pt] V_a(A_i) &= p|d_a (p,A_i) \le d_a(p,A_j), j \ne i, j=1, \cdots, n \ P_\text{hv} &= 1 - \sum_{i=1}^{n} V_a(A_i) \text{ where } h_i \ge h_\text{threshold} \end{align}$$

Beer-Lambert Law canopy gap fraction:

$$\begin{align} L_e &= 2\sum_{i=1}^n - \ln\overline{P(\theta_i)}\cos{\theta} \frac{\sin\theta_i}{\sum_{j=1}^n\sin\theta_j} \ P_o &= \sum_{i=1}^n \exp \left (\frac{-L_eG(\theta_i)}{\cos\theta_i} \right) \end{align}$$

Canopy-to-total-return ratio:

$$\begin{align} VCC_r &= \frac {\sum_{i=1}^n n_{\text{all}i} \gt 1.25 \text{ m}} {\sum{i=1}^n n_{\text{last}i} + n{\text{single}_i}} \end{align}$$

Canopy-to-total-first-return ratio:

$$\begin{align} VCC_\text{fr} &= \frac {\sum_{i=1}^n n_{\text{all}i} \gt 1.25 \text{ m}} {\sum{i=1}^n n_{\text{first}_i}} \end{align}$$

Intensity-return ratio:

$$\begin{align} VCC_\text{ir} &= \frac {\sum_{i=1}^n I_{\text{ground}i} \gt 1.25 \text{ m}} {\sum{i=1}^n I_{\text{all}_i}} \end{align}$$

Beer's Law modified intensity-return ratio:

$$\begin{align} VCC_\text{bl} &= \frac {\frac {\sum_{i=1}^n I_{\text{ground single}i} } {\sum{i=1}^n I_{\text{all}i}} + \sqrt \frac {\sum{i=1}^n I_{\text{ground last}i}} {\sum{i=1}^n I_{\text{all}i}} } {\frac {\sum{i=1}^n I_{\text{first}i} + I{\text{single}i}} {\sum{i=1}^n I_{\text{all}i}} + \sqrt \frac {\sum{i=1}^n I_{\text{intermediate}i} + I{\text{last}i}} {\sum{i=1}^n I_{\text{all}_i}} } \end{align}$$

Above-height cover index:

$$\begin{align} VCC_\text{aci} &= \frac {\sum_{i=1}^n n_{\text{single}i} + n{\text{all}i \gt 1.25 \text{ m}} + n{\text{intermediate}i} + n{\text{last}i}} {\sum{i=1}^n n_{\text{all}_i}} \end{align}$$

First-echo cover index:

$$\begin{align} VCC_\text{fci} &= \frac {\sum_{i=1}^n n_{\text{single}i \gt 1.25 \text{ m}} + n{\text{first}i \gt 1.25 \text{ m}}} {\sum{i=1}^n n_{\text{single}i} + n{\text{first}_i}} \end{align}$$

Solberg's cover index:

$$\begin{align} VCC_\text{sci} &= \frac {\sum_{i=1}^n n_{\text{single}i \gt 1.25 \text{ m}} + \frac{1}{2} (n{\text{first}i \gt 1.25 \text{ m}} + n{\text{last}i \gt 1.25 \text{ m}})} {\sum{i=1}^n n_{\text{single}i} + \frac{1}{2} (n{\text{first}i} + n{\text{last}_i})} \end{align}$$

Canopy-to-total-pixel ratio:

$$\begin{align} VCC_p &= \frac {\sum_{i=1}^n n_{\text{CHM}i \gt 1.25 \text{ m}}} {\sum{i=1}^n n_{\text{CHM}_i}} \end{align}$$

Cartesian-Voronoi fractional cover:

$$\begin{align} V_a(A_i) &= p|d_a (p,A_i) \le d_a(p,A_j), j \ne i, j=1, \cdots, n \ P_\text{cv} &= 1 - \sum_{i=1}^{n} V_a(A_i) \text{ where } h_i \ge h_\text{threshold} \end{align}$$

Getting started

Data pre-processing

By far, the simplest option for pre-processing data is to follow the lidR package tutorial, which utilizes LAStools under the hood. For LiDAR data without ground point classifications, height-normalized point clouds can be produced either with two LAStools command line functions, lasground and lasheight, or with three functions in USDA Fusion, GroundFilter, GridSurfaceCreate, and CanopyModel. If the ground points are already classified then you only need to use the lasheight function of LAStools, while the process for Fusion still requires the same three functions. Hence, I recommend that you use LAStools, as its ground point classification algorithm is also superior to that of Fusion, producing more accurate height-normalized point clouds. This is because Fusion uses the Kraus and Pfeifer (1998) algorithm^[Kraus and Pfeifer (1998) Determination of terrain models in wooded areas with airborne laser scanner data. http://www.sciencedirect.com/science/article/pii/S0924271698000094], while LAStools implements an optimized version of the Axelsson (1999) algorithm^[Axelsson (1999) Processing of laser scanner data—algorithms and applications. http://www.sciencedirect.com/science/article/pii/S0924271699000088]. For more information, read Maguya, Junttila, and Kauranne (2014)^[Maguya, Junttila, and Kauranne (2014) http://www.mdpi.com/2072-4292/6/7/6524]. An example application of lasground and lasheight, implemented in Command Prompt or Bash, is provided below.

lasground -i lidar.las -o lidar_g.las 
lasheight -i lidar_g.las -o lidar_n.las -replace_z 

In order to run these functions, you need to istall LAStools. For Windows, don't forget to add the LAStools bin directory to your system path^[http://www.computerhope.com/issues/ch000549.htm]. For a single LiDAR plot, this is simple to run without leaving your R session. You can call these functions using the system function included in base R, as shown below.

setwd('C:/lidar')
system(lasground -i lidar.las -o lidar_g.las)
system(lasheight -i lidar_g.las -o lidar_n.las -replace_z)

To loop through LAS files in a folder, the syntax follows this pseudocode:

folder <- 'C:/lidar'
files  <- list.files(folder, pattern="\\.las$", full.names=TRUE)

for (i in 1:length(files)) {
  file   <- files[i]
  basenm <- basename(file)
  filenm <- strsplit(basenm,'.',fixed=TRUE)[[1]][1]
  ground <- paste(folder,filenm,'_ground.las',sep='')
  htnorm <- paste(folder,filenm,'_norm.las',sep='')
  system(paste('lasground -i ',file,' -o ',ground, sep=''))
  while (!file.exists(ground)) { Sys.sleep(1) }
  system(paste('lasheight -i ',ground,' -o ',htnorm,' -replace_z', sep=''))
  while (!file.exists(htnorm)) { Sys.sleep(1) }
}

What this for loop does is read in each LAS file path, extract the name of the file without extension, create the filenames of the ground and height-normalized outputs, execute lasground, wait for the output, execute lasheight using the ground file as the input, wait for the output, then proceed to the next iteration. The code should be simple to parallelize using the foreach package, with each fork running in a new Command Prompt or Bash window.

Again, using the lidR package is much simpler!

Example data

After loading the gapfraction package with the library function, the example data can be loaded by calling data(las). The included data consists of fake $X,Y,Z$ coordinates in UTM 11N and meters, along with 8-bit unsigned interger values for intensity and return number. The data consists of one 100-meter diameter LiDAR plot. Based on previous research, I recommend using plots minimally of this size for comparison to ground data (e.g., hemispherical photography) to capture edge effects.

Let's get started!

knitr::opts_chunk$set(fig.width=6, fig.height=6, echo=FALSE, warning=FALSE, message=FALSE)

require(gapfraction)
data(las)
knitr::kable(head(las,10))

Once the data is loaded, you can proceed to call functions from the gapfraction package.

Example usage

Computing hemispherical Voronoi gap fraction

require(gapfraction)
data(las)
P.hv(las=las, model="equidist", thresh.val=1.25, thresh.var="height", reprojection=NA,
     pol.deg=5, azi.deg=45, col="height", plots=T, plots.each=T, plots.save=F)

Comparison of canopy height model (CHM) algorithms

require(gapfraction)
data(las)
chm(las)
chm.sf(las)

Creating a pit-free CHM and performing individual tree crown (ITC) detection with the standard variable-window and watershed algorithms

require(gapfraction)
data(las)
chm  <- chm.sf(las, silent=TRUE)
mw   <- itc.mw(chm, ht2rad=function(x) 0.15746*x)
wat  <- itc.wat(chm, ht2rad=function(x) 0.15746*x)

Creating a stacked pit-free CHM and performing individual tree crown (ITC) detection with the hierarchical variable-window and watershed algorithms

require(gapfraction)
data(las)
chm  <- chm.sf(las, stacked=TRUE, silent=TRUE)
mw   <- itc.mw.h(chm, ht2rad=function(x) 0.15746*x, silent=TRUE)
wat  <- itc.wat.h(chm, ht2rad=function(x) 0.15746*x, silent=TRUE)

The itc algorithms store values for trees and crown.area in a two-element named numeric vector.

Happy laser scanning :)

adam-erickson/gapfraction documentation built on May 5, 2019, 6:57 p.m.