Creating a 3D voxel plot with LAD values in lidR

Question

I am interested in the vertical distribution of the leaf area density (LAD) in my plot. That is why I am seeking to create voxels that represent LAD. I've been trying to achieve this with the voxel_metrics() function. According to the documentation of this function the inputs necessary are a LAS object, a function to be applied to the voxels, and a resolution of the voxels. The functions that can be applied to the voxels work the same way as the grid_metrics() function. One of the pre-made functions is the LAD function which computes a LAD value per voxel.

My limited experience in scripting got me this far:

library(lidR)
library(lidRplugins)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las = readLAS(LASfile)
voxelsLAD = voxel_metrics(las, ~LAD(z=las@data$Z, dz = 1, k= 0.5, z0= 1), 5)
plot(voxelsLAD, color = "lad", colorPalette = heat.colors(5), bg = "black", legend = TRUE)

The resulting LAD values range from 0.00009805‬ to 0.247 which are in line with the values retrieved by the original research of which the function is based on. So far so good you would say, however plotting these values turned out to be difficult. The colorpalette does not apply itself to the voxels as can be seen in this image (all voxels retrieve the same color):

Right now I'm not sure whether this problem can be fixed by computing the LAD voxels otherwise (I might be misusing voxel_metrics()), or by plotting the voxels otherwise.
FYI; I'm using Windows 10 with R 3.5.0

Answer 1

Your question is closely related to this one: Apply a function to a catalog. It is the exact same issue.

There are several confusion in your example. You tried to use the "metrics syntax" in a wrong way. When defining a "metric" function the input should be some attributes of a LAS object and the output must be a number or a list of numbers.

LAD(z=las@data$Z, dz = 1, k= 0.5, z0= 1)

Is a valid statement (you can write las$Z) but the output is not a number. The output is a data.frame is thus it is not a metrics and it cannot be used in *metrics functions.

LAD has been designed from Bouvier et al. (2015). If you read the paper the metric used was the coefficient of variation of the LAD above 2m (if I remember well).

cv(LAD(z, dz = 1, k= 0.5, z0= 2)$lad)

might be a valid metric. The following works (notice the use of uppercase Z. This is an attribute of the LAS object):

library(lidR)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las = readLAS(LASfile)
pixelsLAD = grid_metrics(las, ~as.numeric(cv(LAD(Z, dz = 1, k= 0.5)$lad)), 20)
plot(pixelsLAD)

However it makes sense at the pixel level but it does not really make sense at the voxel level

voxelsLAD = voxel_metrics(las, ~as.numeric(cv(LAD(Z, dz = 1, k= 0.5)$lad)), 5)

won't crash because the computation is valid and works but the numbers won't really have a valid meaning in my opnion (edit: it returns NA everywhere actually).

In your case you wrote:

voxelsLAD = voxel_metrics(las, ~LAD(z=las@data$Z, dz = 1, k= 0.5, z0= 1), 5)

This did not failed because R is very permissive but what it actually does is the computation of the same data.frame from the whole point cloud as many times as the number of voxel. All the data.frame were rbinded.

---

Bouvier, M., Durrieu, S., Fournier, R. a, & Renaud, J. (2015). Generalizing predictive models of forest inventory attributes using an area-based approach with airborne las data. Remote Sensing of Environment, 156, 322-334. http://doi.org/10.1016/j.rse.2014.10.004