R/canopyporosityfilledvolume.R

In akamoske/LiDARforestR: R library for estimating leaf area density (LAD) and leaf area index (LAI) from airborne LiDAR point clouds

#' Convert LAD estimates into two rasters - volume of filled canopy and volume of canopy porosity
#'
#' This function reads in a the LAD estimates that were previously calculated, 
#' finds the volume of voxels in a given column that contain a LAD estimate and the volume of 
#' voxels in a given column that are empty (i.e. no LAD estimates). The output is a list 
#' containing two rasters, one for each calcuation.
#' 
#' These forest structure attributes are based off calculations from: 
#' 
#' Hardiman, B., Bohrer, G., Gough, C., & Curtis, P. (2013). 
#' Canopy structural changes following widespread mortality of canopy dominant trees. 
#' Forests, 4, 537-552. https://doi.org/10.3390/f4030537
#' 
#' Lefsky, M.A., Cohen, W.B., Acker, S.A., Parker, G.G., Spies, T.A., and Harding, D. (1999).
#' Lidar Remote Sensing of the Canopy Structure and Biophysical Properties of Douglas-Fir Western Hemlock Forests.
#' Remote Sensing of the Environment, 70, 339-361. https://doi.org/10.1016/S0034-4257(99)00052-8
#' 
#' @param lad.array LAD estimate array that was generated using the machorn.lad function. 
#' @param laz.array Voxelized LiDAR array that was generated using the laz.to.array function. This contains
#' spatial information for all arrays.
#' @param ht.cut Height that calculations will exclude. This is to remove understory LAD estimates from
#' further calculations. If 5 is entered then all voxels 5 meters and above will be included. Enter 0 if
#' you want to include all calculations
#' @param xy.res Resolution of xy coordinates - if it is a 10x10 meter pixel then enter 10 here
#' @param z.res Vertical resolution of voxel - if it is 1 meter tall then enter 1
#' @param epsg.code EPSG code so that the rasters can be projected into the appropriate projection
#' @return A list containing max LAD and height of max LAD rasters.
#' @export

canopy.porosity.filled.volume <- function(lad.array, laz.array, ht.cut, xy.res, z.res, epsg.code) {
  
  #Lets create an empty matrix that corresponses with this row of data
  porosity.mat <- matrix(data = NA, nrow = dim(lad.array$rLAD)[2], ncol = dim(lad.array$rLAD)[3])
  filled.mat <- matrix(data = NA, nrow = dim(lad.array$rLAD)[2], ncol = dim(lad.array$rLAD)[3])
  
  #loop through the array and calculate the volume of filled and empty voxels
  for (r in 1:dim(lad.array$rLAD)[2]) {
    for (c in 1:dim(lad.array$rLAD)[3]) {
      
      canopy.column <- lad.array$rLAD[(ht.cut + 1):dim(lad.array$rLAD)[1],r,c]
      
      filled <- sum(canopy.column > 0, na.rm = TRUE)
      porosity <- sum(canopy.column == 0, na.rm = TRUE)
      
      filled.volume <- round(filled * xy.res * xy.res * z.res, digits = 4)
      porosity.volume <- round(porosity * xy.res * xy.res * z.res, digits = 4)
      
      filled.mat[r,c] <- filled.volume
      porosity.mat[r,c] <- porosity.volume
    }
  }
  
  #this is the projection code for your site, change it as needed.
  crs.proj <- base::paste0("+init=epsg:", epsg.code)
  
  #now we can create a raster for each level of the canopy using the original x,y data from the laz data
  filled.raster <- raster::raster(filled.mat,
                                  xmn = laz.array$x.bin[1],
                                  xmx = laz.array$x.bin[length(laz.array$x.bin)],
                                  ymn = laz.array$y.bin[1],
                                  ymx = laz.array$y.bin[length(laz.array$y.bin)],
                                  crs = crs.proj)
  
  empty.raster <- raster::raster(porosity.mat,
                                 xmn = laz.array$x.bin[1],
                                 xmx = laz.array$x.bin[length(laz.array$x.bin)],
                                 ymn = laz.array$y.bin[1],
                                 ymx = laz.array$y.bin[length(laz.array$y.bin)],
                                 crs = crs.proj)
  
  #we have to flip these rasters so that they are orientated the correct direction
  #this is standard when converting from an array to a raster
  filled.raster.flip <- flip(filled.raster, direction = "y")
  empty.raster.flip <- flip(empty.raster, direction = "y")
  
  #return the final rasters
  final.data <- list("filled.raster" = filled.raster.flip, 
                     "porosity.raster" = empty.raster.flip)
  return(final.data)
}