global(x, "notNA") yields incorrect results

[jeffreyhanson]

Hi,

Sorry to be a pain, but I think I might have found another bug in the CRAN version of terra. I can also reproduce this with the latest version on GitHub too. It seems that global(x, "notNA") yields an incorrect result, and that na.rm = TRUE must be supplied for correct results. Is this a bug or intended behavior? If the latter, it might be worth mentioning this in the documentation. I've included a reprex below, and also my session information. Please let me know if there's any other information I can provide?

# load package
library(terra)
#> terra 1.7.24

# create SpatRaster with only NA cells,
r <- rast(matrix(NA, ncol = 3, nrow = 3))

# calculate number of not NA cells
## since the raster only has NA cells the correct answer is 0
global(r, "notNA")
#>      notNA
#> lyr.1     9

# calculate number of not NA cells, with na.rm = TRUE
global(r, "notNA", na.rm = TRUE)
#>      notNA
#> lyr.1     0

Session information

R version 4.2.3 (2023-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS

Matrix products: default
BLAS:   /opt/R/R-4.2.3/lib/R/lib/libRblas.so
LAPACK: /opt/R/R-4.2.3/lib/R/lib/libRlapack.so

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       

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

other attached packages:
[1] terra_1.7-24   testthat_3.1.7 devtools_2.4.5 usethis_2.1.6 

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

[rhijmans]

Thank you very much for reporting, and please do not apologize for that. I am sorry for that bug, but is fixed now. I was introduced when I changed global to allow for supplying multiple summarizing functions, which is more efficient than doing them one by one (see below).

This works now:

r <- rast(matrix(NA, ncol = 3, nrow = 3))
global(r, "notNA")
#      notNA
#lyr.1     0
global(r, "notNA", na.rm = TRUE)
#      notNA
#lyr.1     0
global(r, "isNA")
#      isNA
#lyr.1    9
global(r, "isNA", na.rm = TRUE)
#      isNA
#lyr.1    9

As does this

global(r, c("notNA", "isNA", "range", "mean"), na.rm = TRUE)
#      notNA isNA mean min max
#lyr.1     0    9  NaN  NA  NA

global(r, c("notNA", "isNA", "range", "mean"), na.rm = FALSE)
#      notNA isNA mean min max
#lyr.1     0    9   NA  NA  NA

r[1:3] <- 1:3
global(r, c("notNA", "isNA", "range", "mean"), na.rm = TRUE)
#      notNA isNA mean min max
#lyr.1     3    6    2   1   3

global(r, c("notNA", "isNA", "range", "mean"), na.rm =FASLSE)
#      notNA isNA mean min max
#lyr.1     3    6  NaN NaN NaN

[jeffreyhanson]

Awesome - thanks for fixing this so quickly!

[fgoerlich]

Hi,
The issue appears not to be fully solved for big raster files.
The raster does not have NAs.

> library(terra)
terra 1.7.24
> r <- rast("hisdac_es_evol_bufa_v1_100_1900.tif")
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
> ncell(r)
[1] 367452000
> global(r, fun = "isNA")
                                isNA
hisdac_es_evol_bufa_v1_100_1900    0
> global(r, fun = "isNA",  na.rm = TRUE)
                                      isNA
hisdac_es_evol_bufa_v1_100_1900 -387488400
> global(r, fun = "notNA")
                                    notNA
hisdac_es_evol_bufa_v1_100_1900 367452000
> global(r, fun = "notNA", na.rm = TRUE)
                                    notNA
hisdac_es_evol_bufa_v1_100_1900 755365200

The na.rm = TRUE option gives incorrect results.
A negative value for global(r, fun = "isNA", na.rm = TRUE) in nonsense.
Supplying multiple summarizing functions I get the same

> global(r, c("notNA", "isNA", "range", "mean"), na.rm = FALSE)
                                    notNA isNA      mean min  max
hisdac_es_evol_bufa_v1_100_1900 367452000    0 0.6550164   0 9992
> global(r, c("notNA", "isNA", "range", "mean"), na.rm = TRUE)
                                    notNA       isNA      mean min  max
hisdac_es_evol_bufa_v1_100_1900 756427200 -388975200 0.6550164   0 9992

Moreover I don´t always get the same number. Compare the single summarizing functions, with the multiple case.

With version 1.7.23 results were incorrect independently of the na.rm option.

The (zipped) file is attached below.

Session information

R version 4.2.3 (2023-03-15 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 8.1 x64 (build 9600)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252    LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C                   LC_TIME=Spanish_Spain.1252    

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

other attached packages:
[1] terra_1.7-24

loaded via a namespace (and not attached):
[1] compiler_4.2.3   tools_4.2.3      Rcpp_1.0.10      codetools_0.2-19

NB: The warnings after reading, rast, are related to r-spatial/stars#588, I don´t know why I get them and I have been unable to drop.
Thanks a lot!
hisdac_es_evol_bufa_v1_100_1900.zip

[rhijmans]

Thank you! isNA did not work when processing a file in chunks. It works now, I believe:

library(terra)
setwd("../Downloads")
f <- system.file("ex/elev.tif", package="terra")
x <- rast(f)
funs <- c("isNA", "notNA", "min", "mean", "sum", "sd")

terraOptions(steps=1)
global(x, fun=funs)
#          isNA notNA min mean sum  sd
#elevation 3942  4608 NaN  NaN NaN NaN

global(x, fun=funs, na.rm=TRUE)
#          isNA notNA min     mean     sum       sd
#elevation 3942  4608 141 348.3366 1605135 80.21886
 
terraOptions(steps=4)
global(x, fun=funs)
#          isNA notNA min mean sum  sd
#elevation 3942  4608 NaN  NaN NaN NaN

global(x, fun=funs, na.rm=TRUE)
#          isNA notNA min     mean     sum       sd
#elevation 3942  4608 141 348.3366 1605135 80.21886

[fgoerlich]

Thanks for the super-quick response!!!
You mean does it work if I install the devel version from r-universe?

[rhijmans]

Yes, but it can take a couple of hours before R-Universe has rebuilt the package (you can check on their site)