r-lidar / lidR

Airborne LiDAR data manipulation and visualisation for forestry application
https://CRAN.R-project.org/package=lidR
GNU General Public License v3.0
582 stars 130 forks source link

Error in pixel_metrics when called with stdmetrics_i applied on a subset which is empty on some pixels #747

Closed jmmonnet closed 5 months ago

jmmonnet commented 5 months ago

I would like to compute with a single call to pixel_metrics two types of variables :

In case some pixels have no points above the height threshold, pixel_metrics returns an error Below is an example. The test.laz data is available at: https://filesender.renater.fr/?s=download&token=d2fa32eb-2894-4140-b24e-b69e57f85cc9

a <- lidR::readALSLAS("test.laz")

user_func <- function(i, z, hmin = 2) {
  # total point number
  ntot = list(ntot = length(z)) 
  #
  # points above hmin
  phmin <- which(z >= hmin)
  # return ntot and std_metrics of points above hmin
  return(c(lidR::stdmetrics_i(i[phmin]), ntot))
}

# hmin at 0 -> OK
test <- lidR::pixel_metrics(a,
                            ~ user_func(Intensity, Z, hmin = 0),
                            res = 20)

# hmin at 100 -> OK
test <- lidR::pixel_metrics(a,
                            ~ user_func(Intensity, Z, hmin = 100),
                            res = 20)

# hmin at 2 -> fails (some pixels have points, some don't)
test <- lidR::pixel_metrics(a,
                            ~ user_func(Intensity, Z, hmin = 2),
                            res = 20)

I could compute the metrics with two calls to pixel_metrics (one with no filter and one with a filter filter = ~Z >= 2, but I guess it would be more efficient to include all metrics in a single call.

Jean-Romain commented 5 months ago

Here I think it should be your responsibility to check how many points are remaining. user_func() is fed with points (never 0). Then you removed some of them and fed stdmetrics_i() with 0 point. stdmetrics_i() and others could handle this case of course but here I'd rather try to return NULL if there is 0 point.

jmmonnet commented 5 months ago

Thank you for the advice. I will add a test in the beginning of user_func() and return NULL in case no points are fed to stdmetrics_i. It is more straightforward this way, the only drawback I see is that user_func() will not return any value for the other metrics (ntot in the example), even tough it could have been computed.

Jean-Romain commented 5 months ago

I see. In this case it is indeed better if the case is handled internally. I'll make a modification.

Jean-Romain commented 5 months ago

I thought about it and actually to handle the case stdmetrics_* functions must return a list with all the metrics set to NA. However, it is not straightforward to know how many metrics are supposed to be returned without computing them. It is not hard, but not straighforward. I can't simply add a close at the beginning of the functions like

if (lengh(i) == 0)
{
  return (list(... ))
}

This is definitively not a priority to me but I'm keen to accept a PR.