zonal(..., wide = FALSE) produces incorrect result

[jeffreyhanson]

Hi,

I think I might have found a bug when using the zonal() function with wide = FALSE. To illustrate the potential bug, please run the code below with the attached files? I've verified that this behavior is present with the latest version of terra on GitHub (1.7-41) and the latest version on CRAN (1.7-39). I have also included my session information below (with the latest GitHub version of terra).

Please let me know if there's any other information I can provide to help with this?

Reproducible example

# load packages
library(terra)

# load data
x <- terra::rast("pu.tif")
z <- terra::rast("country.tif")

# compute zonal statistics
z1 <- terra::zonal(x, z = z, fun = "sum", na.rm = TRUE, wide = TRUE)
z2 <- terra::zonal(x, z = z, fun = "sum", na.rm = TRUE, wide = FALSE)

# manaully compute correct result for 24th zone
terra::global(x * (z == 24), "sum", na.rm = TRUE)
#>                     sum
#> dem-100m-esri54017 1809

# display result for 24th zone computed with wide = TRUE
z1[z1$layer == 24, ]
#>    layer dem-100m-esri54017
#> 23    24               1809

# display result for 24th zone computed with wide = FALSE
z2[z2$id == 24, ]
#>      value layer id
#> 24.1    73     1 24

Data files

terra-issue.zip

Session information

R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS

Matrix products: default
BLAS:   /opt/R/R-4.3.1/lib/R/lib/libRblas.so 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Pacific/Auckland
tzcode source: system (glibc)

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

other attached packages:
[1] terra_1.7-41    testthat_3.1.10 devtools_2.4.5  usethis_2.2.2  

loaded via a namespace (and not attached):
 [1] miniUI_0.1.1.1    compiler_4.3.1    crayon_1.5.2      brio_1.1.3       
 [5] promises_1.2.0.1  Rcpp_1.0.11       stringr_1.5.0     callr_3.7.3      
 [9] later_1.3.1       fastmap_1.1.1     mime_0.12         R6_2.5.1         
[13] htmlwidgets_1.6.2 profvis_0.3.8     shiny_1.7.4.1     rlang_1.1.1      
[17] cachem_1.0.8      stringi_1.7.12    httpuv_1.6.11     fs_1.6.2         
[21] pkgload_1.3.2.1   memoise_2.0.1     cli_3.6.1         magrittr_2.0.3   
[25] ps_1.7.5          digest_0.6.33     processx_3.8.2    xtable_1.8-4     
[29] remotes_2.4.2     lifecycle_1.0.3   prettyunits_1.1.1 vctrs_0.6.3      
[33] glue_1.6.2        urlchecker_1.0.1  codetools_0.2-19  sessioninfo_1.2.2
[37] pkgbuild_1.4.2    purrr_1.0.1       tools_4.3.1       ellipsis_0.3.2   
[41] htmltools_0.5.5  

[kadyb]

It looks like the ID in wide = FALSE is generated by 1:n, not from x values (there isn't ID = 14).

z1 <- terra::zonal(x, z = z, fun = "sum", na.rm = TRUE, wide = TRUE)
z2 <- terra::zonal(x, z = z, fun = "sum", na.rm = TRUE, wide = FALSE)
z2$id <- unique(z)$layer
z2[z2$id == 24, ]
#>      value layer id
#> 23.1  1809     1 24

[rhijmans]

Looks like this bug only presented itself if the layer name of the zones SpatRaster was "layer"; as in your example.

[jeffreyhanson]

Ah - I see - thank you very much @kadyb and @rhijmans!