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:

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):

all voxels end up with 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

  • 1
    This question will help you. It is basically the same user mistake gis.stackexchange.com/questions/342578/…
    – JRR
    Commented Dec 10, 2019 at 15:07
  • 1
    Do not refer to the documentation of lidR 1.6.1. It is fully outdated! If you prefer to read documentation on third party website this is your choice but please pay attention at reading an up-to-date version of the documentation.
    – JRR
    Commented Dec 10, 2019 at 16:56

1 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):

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)

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

  • Thanks for your detailed explanation JRR. I am starting to understand why my code is wrong. One thing that is not completely clear to me is if it is possible to modify the LAD function in a way that it is applicable to voxels, or is that complete out of the question when using the Bouvier et al. approach?
    – Lout
    Commented Dec 11, 2019 at 13:35
  • You are free to do whatever you want and take inspiration from the current LAD function and the paper. But as it is in the package the function is not intended to be used in voxels.
    – JRR
    Commented Dec 11, 2019 at 14:02
  • Alright, thanks for clarifying! In the documentation it stated that the metrics for the grid_metrics3d can be applied in the same way as for the grid_metrics, that is why I figured applying LAD to the voxels would work right off the bat.
    – Lout
    Commented Dec 12, 2019 at 8:44
  • 1
    I will clarify the doc. It means that you can apply the same computations the same way but you are the only one responsible of ensuring your computation has a meaning. For example grid_metrics( ~mean(Z)) and grid_metrics3d(~mean(Z)) are both valid but the second one roughly gives you the position of the voxel and is useless.
    – JRR
    Commented Dec 12, 2019 at 9:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.