I am attempt to write a function to calculate leaf area index using discrete return aerial lidar data using the lidR package in R. My function follows the Beer-Lambert Law with the assumption of randomly oriented leaves so that the Ross-Nilson gap function value = 0.5 in the calculation of the extinction coefficient. With discrete return data, I am calculating the ratio of intercepted to transmitted photosynthetically active radiation as the ratio of canopy returns (defined as above 3 meters) to ground returns (defined as below 2 meters). I have created a "low-level api" function called Calc_LAI()
that does the following:
pixel_metrics()
I am stumped by two issues:
First, I have a helper function called lai_func() that does the arithmetical calculation of leaf area index. Currently I have it nested in the conditional, for if(is(las,"LAS")){}.
However, when I run the whole Calc_LAI()
function, I get an error saying that the lai_func()
isn't found. If I move the helper function out of Calc_LAI()
this error goes away, but I have had helper functions in other "low-level api" calls to catalog_apply()
where they were nested within the if(is(las,"LAS")){}
conditional that worked without error.
Second, the call to pixel_metrics()
results in the following error Error: Duplicated elements found. At least one of the metrics was not a number. Each metric should be a single number.
which perplexes me because I ensure that I only get real values from the calculation of leaf area index by adding 1 to both the numerator and denominator in the ratio of canopy returns to ground returns to ensure that I never divide by zero or take the natural log of 0.
Below is my reproducible example:
##Loading Necessary Package##
library(lidR)
##Reading in Example Lidar Dataset##
lasfile<- system.file("extdata", "MixedConifer.laz", package="lidR")
MixedConifer<-readLAS(lasfile)
##Low-level API function to calculate leaf area index##
Calc_LAI <- function(las)
{
if (is(las, "LAScatalog")) {
options <- list(automerge = TRUE, need_buffer = TRUE)
LAI <- catalog_apply(las, Calc_LAI, .options = options)
return(LAI)
}
else if (is(las, "LAScluster")) {
bbox <- st_bbox(las)
las <- readLAS(las)
if (is.empty(las)) return(NULL)
LAI <- Calc_LAI(las, param1, param2)
LAI <- sf::st_crop(LAI, bbox)
return(LAI)
}
else if (is(las, "LAS")) {
las<-add_attribute(las, paste(las$gpstime, las$ReturnNumber, sep="_"), "UID")
dem<-rasterize_terrain(las, res=1, algorithm = knnidw())
norm<-las-dem
lai_func<-function(ScanAngle, Z, UID){
nground<-length(unique(UID[Z <= 6.6]))
nveg<-length(unique(UID[Z >= 9.9]))
out<-log(nveg+1/(nground+1))*cos((ScanAngle*pi/180)/0.5)
return(list(LAI=out))
}
LAI<-pixel_metrics(las=norm, func=lai_func(ScanAngle=ScanAngleRank, Z=Z, UID=UID), res=9.9)
return(LAI)
}
else {
stop("This type is not supported.")
}
}
##First Attempt with lai_func() nested within Calc_LAI()##
Calc_LAI(MixedConifer)#throws error that lai_func() not found
##Second Attempt with lai_func() outside of Calc_LAI##
lai_func<-function(ScanAngle, Z, UID){
nground<-length(unique(UID[Z <= 6.6]))
nveg<-length(unique(UID[Z >= 9.9]))
out<-log(nveg+1/(nground+1))*cos((ScanAngle*pi/180)/0.5)
return(list(LAI=out))
}
Calc_LAI(MixedConifer)#throws error that Duplicate elements are found and at least one metric was not a number
##Running the function elements directly on example data set
las<-MixedConifer
las<-add_attribute(las, paste(las$gpstime, las$ReturnNumber, sep="_"), "UID")
dem<-rasterize_terrain(las, res=1, algorithm = knnidw())
norm<-las-dem
lai_func<-function(ScanAngle, Z, UID){
nground<-length(unique(UID[Z <= 6.6]))
nveg<-length(unique(UID[Z >= 9.9]))
out<-log(nveg+1/(nground+1))*cos((ScanAngle*pi/180)/0.5)
return(list(LAI=out))
}
LAI<-pixel_metrics(las=norm, func=lai_func(ScanAngle=ScanAngleRank, Z=Z, UID=UID), res=9.9) #throw the duplicate elements error again
throws error that Duplicate elements are found and at least one metric was not a number
out<-log(nveg+1/(nground+1))*cos((ScanAngle*pi/180)/0.5)
return(list(LAI=out))
Here, clearly LAI
is not a metric (i.e a number) but is a vector of size length(ScanAngle)