terra
terra copied to clipboard
rspatial
terra::aggregate - results of `'modal'` are inconsistent with `terra::modal`
It seems that 'modal' returns different results than terra::modal/raster::modal when used as the function to apply within terra::aggregate. This is particularly apparent when the raster contains NAs, but also occurs in the absence of NAs. I'm not sure whether in the absence of NAs the discrepancy may be a result of the treatment of ties, but I suspect not.
library(raster)
#> Loading required package: sp
library(terra)
#> terra 1.7.25
r1 <- rast(system.file("external/rlogo.grd", package="raster"), lyrs='red')
set.seed(1)
r1a <- terra::aggregate(r1, fact=5, 'modal', na.rm=TRUE, ties='first')
r1b <- terra::aggregate(r1, fact=5, terra::modal, na.rm=TRUE, ties='first')
plot(c(r1a, r1b, r1b - r1a), main=c('r1a', 'r1b', 'r1b minus r1a'))
When the raster has NAs, I expect the result to ignore those NAs when na.rm is TRUE, returning the mode of the non-NA cells if there are any, or returning NA otherwise. This is the behaviour we see with terra::modal (though I have a suspicion that my terra::modal actually deploys raster::modal below, since terra::modal is unhappy unless we first load the raster package; https://github.com/rspatial/terra/issues/1121).
However, when using 'modal', the NA values returned by terra::modal are now floats that are close to zero (see head below).
r2 <- r1
r2[r2 > 240] <- NA
r2a <- terra::aggregate(r2, fact=5, 'modal', na.rm=TRUE, ties='first')
r2b <- terra::aggregate(r2, fact=5, terra::modal, na.rm=TRUE, ties='first')
plot(c(is.na(r2a), is.na(r2b), r2b - r2a),
main=c('is.na(r2a)', 'is.na(r2b)', 'r2b minus r2a'))
head(c(r2a, r2b))
#> red red
#> 1 2.470328e-323 NA
#> 2 4.940656e-323 NA
#> 3 7.410985e-323 NA
#> 4 9.881313e-323 NA
#> 5 1.235164e-322 NA
#> 6 2.280000e+02 237
I'd like to use 'modal' since it is considerably faster, but I'd first like to understand why it's behaving the way it does.
Created on 2023-04-19 with reprex v2.0.2
Session info
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.2.1 (2022-06-23)
#> os macOS Big Sur ... 10.16
#> system x86_64, darwin17.0
#> ui X11
#> language (EN)
#> collate en_US.UTF-8
#> ctype en_US.UTF-8
#> tz Australia/Melbourne
#> date 2023-04-19
#> pandoc 2.19.2 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date (UTC) lib source
#> cli 3.6.0 2023-01-09 [1] CRAN (R 4.2.1)
#> codetools 0.2-18 2020-11-04 [1] CRAN (R 4.2.1)
#> curl 4.3.3 2022-10-06 [1] CRAN (R 4.2.0)
#> digest 0.6.31 2022-12-11 [1] CRAN (R 4.2.0)
#> evaluate 0.19 2022-12-13 [1] CRAN (R 4.2.0)
#> fansi 1.0.3 2022-03-24 [1] CRAN (R 4.2.0)
#> fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.2.0)
#> fs 1.5.2 2021-12-08 [1] CRAN (R 4.2.0)
#> glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.0)
#> highr 0.10 2022-12-22 [1] CRAN (R 4.2.0)
#> htmltools 0.5.4 2022-12-07 [1] CRAN (R 4.2.0)
#> httr 1.4.4 2022-08-17 [1] CRAN (R 4.2.1)
#> knitr 1.41 2022-11-18 [1] CRAN (R 4.2.0)
#> lattice 0.20-45 2021-09-22 [1] CRAN (R 4.2.1)
#> lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.2.0)
#> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.0)
#> mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)
#> pillar 1.8.1 2022-08-19 [1] CRAN (R 4.2.0)
#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.0)
#> purrr 1.0.0 2022-12-20 [1] CRAN (R 4.2.0)
#> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.2.0)
#> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.2.0)
#> R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.2.0)
#> R.utils 2.12.2 2022-11-11 [1] CRAN (R 4.2.0)
#> R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.0)
#> raster * 3.6-13 2023-01-07 [1] CRAN (R 4.2.0)
#> Rcpp 1.0.10 2023-01-22 [1] CRAN (R 4.2.0)
#> reprex 2.0.2 2022-08-17 [1] CRAN (R 4.2.0)
#> rlang 1.0.6 2022-09-24 [1] CRAN (R 4.2.0)
#> rmarkdown 2.19 2022-12-15 [1] CRAN (R 4.2.0)
#> rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.0)
#> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.0)
#> sp * 1.5-1 2022-11-07 [1] CRAN (R 4.2.0)
#> stringi 1.7.8 2022-07-11 [1] CRAN (R 4.2.0)
#> stringr 1.5.0 2022-12-02 [1] CRAN (R 4.2.0)
#> styler 1.7.0 2022-03-13 [1] CRAN (R 4.2.0)
#> terra * 1.7-25 2023-04-19 [1] Github (rspatial/terra@fa8195b)
#> tibble 3.1.8 2022-07-22 [1] CRAN (R 4.2.0)
#> utf8 1.2.2 2021-07-24 [1] CRAN (R 4.2.0)
#> vctrs 0.5.1 2022-11-16 [1] CRAN (R 4.2.0)
#> withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.0)
#> xfun 0.36 2022-12-21 [1] CRAN (R 4.2.0)
#> xml2 1.3.3 2021-11-30 [1] CRAN (R 4.2.0)
#> yaml 2.3.6 2022-10-18 [1] CRAN (R 4.2.0)
#>
#> [1] /Library/Frameworks/R.framework/Versions/4.2/Resources/library
#>
#> ──────────────────────────────────────────────────────────────────────────────
I do think the difference you see is due to ties. I do not see a difference when using ties='lowest'
library(raster)
library(terra)
r1 <- rast(system.file("external/rlogo.grd", package="raster"), lyrs='red')
set.seed(1)
r1a <- terra::aggregate(r1, fact=5, 'modal', na.rm=TRUE)
r1b <- terra::aggregate(r1, fact=5, raster::modal, na.rm=TRUE, ties='lowest')
all.equal(r1a, r1b, maxcell=Inf)
#[1] TRUE
minmax(r1a - r1b)
minmax(r1a - r1b)
# red
#min 0
#max 0
What needs to happen, I think is to allow providing a "ties" argument to aggregate(<SpatRaster>, fun='modal')
Thanks Robert. That won't explain the differing treatment of NAs, though, will it (the second example above)?
Thanks, I overlooked that one, the terra method now also returns NAs (not zeros) for cells where there are no values.