pixel_metrics() output values error

[Rpouliot]

Working on a ABA approach, the generation of the suite of ALS metrics on my files collection worked, but it ended up with absurd results, with improbable values and the same value assigned to all output metrics.

Here's a reproductible exemple of the issue using pixel_metrics().

library(lidR)
las_mega <- readLAScatalog("~/win-library/4.3/lidR/extdata/Megaplot.laz")
plot(las_mega)

Change the chunk size to make the process like a lascatalog

opt_chunk_size(las_mega) <- 60
opt_chunk_buffer(las_mega) <- 10
plot(las_mega, chunk = TRUE)

Run pixel_metrics()

w2w_mega_pixel <- pixel_metrics(las_mega, .stdmetrics_z, res = 20)
pixel_metrics_df <- as.data.frame(w2w_mega_pixel)

If I do the same process with grid_metrics() it looks like it works properly:

w2w_mega_grid <- grid_metrics(las_mega, .stdmetrics_z, res = 20)
grid_metrics_df <- as.data.frame(w2w_mega_grid)

[Jean-Romain]

Can you produce a fully reproducible example with R version, terra version, raster version and the list of loaded packages. grid_metrics() simply calls pixel_metrics() internally with pkg = "raster" so I doubt using pixel_metrics() or grid_metrics() can change anything. And everything is fine on my side with pixel metrics().

library(lidR)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las_mega <- readLAScatalog(LASfile)
plot(las_mega)

opt_chunk_size(las_mega) <- 60
#> Be careful, a chunk size smaller than 250 is likely to be irrelevant.
opt_chunk_buffer(las_mega) <- 10
plot(las_mega, chunk = TRUE)

w2w_mega_pixel <- pixel_metrics(las_mega, .stdmetrics_z, res = 20)
plot(w2w_mega_pixel)

pixel_metrics_df <- as.data.frame(w2w_mega_pixel)
pixel_metrics_df
#>      zmax        zmean         zsd       zskew     zkurt   zentropy
#> 1   22.00 12.328000000 7.063748852 -0.48545924  1.756451 0.92847162
#> 2   27.20 20.421703704 5.308555771 -2.06926888  7.573807 0.78415486
#> 3   25.85 17.181541219 6.043577314 -1.25369318  4.176640 0.88843882
#> 4   22.46 15.158794326 4.974598963 -1.60682468  5.468312 0.82511324
#> 5   21.66 14.667542662 4.914715822 -0.90663771  3.801231 0.88021647
#> 6   22.20 14.596972789 5.484525611 -1.28607415  4.480086 0.82246798
#> 7   21.91 14.193561151 5.746722806 -0.74878762  2.926346 0.89451150
#> 8   19.12  9.297339181 5.732670706 -0.24385137  1.886132 0.90953538
#> 9   24.98 14.298700906 6.177650270 -0.63916040  2.799621 0.93554820
#> 10  24.08 15.294304348 6.131553597 -0.97785133  3.241297 0.90512689
[...]

^{Created on 2023-08-15 with reprex v2.0.2}

[flottsam]

I was about to post the exact same issue. For me, grid_metricsworks as expected, but with pixel_metrics, every output metric is exactly the same.

Note that the issue only occurs with multiple files using las catalog, a single file with las catalog works as expected. In the example above, you are reading a single file (on a side note, it would be useful to have a couple of small system files for reproducing las catalog issues/questions/bugs).

Two files I'm using can be accessed here: https://drive.google.com/drive/folders/1ujMK7jF_nICD-VCq88TLF7_GQoqU4Fbl?usp=sharing

Here's an example of my code:

# packages I typically have loaded: 
if (!require("pacman")) install.packages("pacman"); library(pacman)
p_load(lidR, here, tidyverse, sf, terra, magrittr, future, parallel, lmom)

# define function:
myMetrics <- function(x,y,z, rn, i, class){
  first  = rn == 1L
  zfirst = z[first]
  nfirst = length(zfirst)
  firstabove2 = sum(zfirst > 2)
  nz = length(z)
  metrics = list(
    # vertical complexity index (all strata)
    vci = VCI(z, zmax=max(z), by=1),
    # canopy cover (%)
    canopy_cover = ((firstabove2/nfirst)*100),
    # get mean values of returns 
    z_mean = mean(z),
    # standard deviation of returns 
    z_sd = sd(z))
  return(c(metrics))
}

# read files
(ctg <- readLAScatalog(here("path/to/files")))
opt_select(ctg) <- "xyzrnic" 
opt_filter(ctg) <- "-drop_class 6 14" # building & powerlines returns
plan(multisession, workers = detectCores(logical = FALSE)/2) # use 1/2 physical cores

(m <- lidR::grid_metrics (ctg, 
                          ~myMetrics(x=X, y=Y, z=Z, rn=ReturnNumber, i=Intensity, class=Classification),
                          res = 1))

Output using grid_metrics (returns a raster brick):

Output using pixel_metrics (returns a spatrast):

(m <- lidR::pixel_metrics(ctg, 
                          ~myMetrics(x=X, y=Y, z=Z, rn=ReturnNumber, i=Intensity, class=Classification),
                          res = 1))

Output w/ pixel_metrics:

[Jean-Romain]

(on a side note, it would be useful to have a couple of small system files for reproducing las catalog issues/questions/bugs).

You are right but I cannot do anything. I had to reduce drastically the amount of data to fit the package size imposed by the CRAN. Otherwise there would have been much more data available. There are so few data that sometime it is hard to write an example in the documentation :disappointed:

That being said I still cannot reproduce. I used the following code in a fresh session:

myMetrics <- function(x,y,z, rn, i, class){
  first  = rn == 1L
  zfirst = z[first]
  nfirst = length(zfirst)
  firstabove2 = sum(zfirst > 2)
  nz = length(z)
  metrics = list(
    #foliar height diversity
    fhd = (entropy(z, zmax = max(z)) * log(max(z))),  
    # vertical distribution ratio
    vdr = ( (max(z) - median(z)) / max(z)),
    # 'top' rugosity
    top_rug = sd(zfirst),
    # crown relief ratio (i.e. rugosity)
    crr = ((mean(z)-min(z)) / (max(z)-min(z))),
    # vertical complexity index (3m, understory strata)
    uci = VCI(z, zmax=3, by=0.25),
    # vertical complexity index (all strata)
    vci = VCI(z, zmax=max(z), by=1),
    # canopy cover (%)
    canopy_cover = ((firstabove2/nfirst)*100),
    # get mean values of returns 
    z_mean = mean(z),
    # standard deviation of returns 
    z_sd = sd(z))
  return(c(metrics))
}

library(lidR)
library(future)
library(terra)
library(sf)

# read files
ctg <- readLAScatalog("issue 697/")
opt_select(ctg) <- "xyzrnic" 
opt_filter(ctg) <- "-drop_class 6 14" # building & powerlines returns
plan(multisession, workers = 2)

m <- lidR::pixel_metrics(ctg, ~myMetrics(x=X, y=Y, z=Z, rn=ReturnNumber, i=Intensity, class=Classification), res = 1)
m

m is correct.

[flottsam]

Thank you for looking. Just to clarify, you are reading in more than 1 tile with las catalog?

[Jean-Romain]

I'm using the two tiles you shared and I'm getting the correct output you showed

[flottsam]

Hmm ok, I guess I will need to revert back to grid_metrics for now.

[Jean-Romain]

I'd prefer you to help me figuring out what is the problem. Btw grid_metrics does not really exist it is just a call to

res <- pixel_metrics(las, ..., pkg = "raster") 

So it really weird and we should investigate what is going wrong. Please start trying to make the most minimal example possible. Can we simplify the metrics? can we use smaller files? Can we reproduce with no parallelization? No package loaded? Does @Rpouliot's example works for you ? Which packages version are you using? R version? OS?

[flottsam]

It must be related to terra.

I ran pixel_metrics using pkg="raster" and the output was correct (but obviusly returns a raster brick).

m <- 
  lidR::pixel_metrics(ctg, 
                      ~myMetrics(x=X, y=Y, z=Z, rn=ReturnNumber, i=Intensity, class=Classification),
                      res = 1, pkg="raster")

I tried pixel_metrics(pkg='terra') sequentially and still encounter the same issue.

RStudio: 2023.06.1 Build 524 Package versions: Terra: 1.7.39 lidR: 4.0.3 sf: 1.0.14

Session info:

R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8    LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                           LC_TIME=English_United States.utf8    

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] here_1.0.1    sf_1.0-14     terra_1.7-39  future_1.33.0 lidR_4.0.3   

loaded via a namespace (and not attached):
 [1] lwgeom_0.2-13      crayon_1.5.2       dplyr_1.1.2        compiler_4.3.1     tidyselect_1.2.0   Rcpp_1.0.11        parallel_4.3.1    
 [8] progress_1.2.2     globals_0.16.2     lattice_0.21-8     R6_2.5.1           generics_0.1.3     classInt_0.4-9     tibble_3.2.1      
[15] stars_0.6-3        rlas_1.6.3         units_0.8-3        rprojroot_2.0.3    DBI_1.1.3          pillar_1.9.0       rlang_1.1.1       
[22] utf8_1.2.3         sp_2.0-0           lazyeval_0.2.2     cli_3.6.1          magrittr_2.0.3     class_7.3-22       digest_0.6.33     
[29] grid_4.3.1         rstudioapi_0.15.0  hms_1.1.3          lifecycle_1.0.3    prettyunits_1.1.1  vctrs_0.6.3        KernSmooth_2.23-21
[36] proxy_0.4-27       glue_1.6.2         data.table_1.14.8  raster_3.6-23      listenv_0.9.0      codetools_0.2-19   abind_1.4-5       
[43] parallelly_1.36.0  fansi_1.0.4        e1071_1.7-13       tools_4.3.1        pkgconfig_2.0.3   

I will test @Rpouliot data and try and get smaller tiles to work with.

[flottsam]

OK, so I installed terra 1.7.29 and pixel_metricsworks as intended.

[Jean-Romain]

And I installed 1.7.39 and I was able to reproduce. Now I need to figure out what is going wrong in terra to report.

[Jean-Romain]

reported in https://github.com/rspatial/terra/issues/1262