Indexing by other raster returns wrong results

[Nowosad]

The raster package allowed for extracting values of one raster based on the other with the use of the [ operator (see the top example).
However, the same approach does not correctly work for terra (see the bottom example). Is it a bug, or should the users not use this syntax?

set.seed(2017-12-16)
library(raster)
#> Loading required package: sp
elev1 = raster(nrow = 6, ncol = 6, res = 0.5, 
             xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
             vals = 1:36)
clip1 = raster(xmn = 0.9, xmx = 1.8, ymn = -0.45, ymx = 0.45,
             res = 0.3, vals = rep(1, 9))
elev1[clip1]
#> [1] 18 24


set.seed(2017-12-16)
library(terra)
#> terra version 1.4.2
elev2 = rast(nrow = 6, ncol = 6, resolution = 0.5, 
            xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
            vals = 1:36)
clip2 = rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
            resolution = 0.3, vals = rep(1, 9))
elev2[clip2]
#>    lyr.1
#> 1      1
#> 2      2
#> 3      3
#> 4      4
#> 5      5
#> 6      6
#> 7      7
#> 8      8
#> 9      9
#> 10    10
#> 11    11
#> 12    12
#> 13    13
#> 14    14
#> 15    15
#> 16    16
#> 17    17
#> 18    18
#> 19    19
#> 20    20
#> 21    21
#> 22    22
#> 23    23
#> 24    24
#> 25    25
#> 26    26
#> 27    27
#> 28    28
#> 29    29
#> 30    30
#> 31    31
#> 32    32
#> 33    33
#> 34    34
#> 35    35
#> 36    36

^{Created on 2021-09-12 by the reprex package (v2.0.1)}

[Nowosad]

On a similar note, raster allows to keep a RasterLayer structure with drop = FALSE, which does not seem to be working with terra (see my examples below). Is this an expected behavior? I was unable to find anything about it at https://rspatial.github.io/terra/reference/terra-package.html#xxiv-changed-behavior.

set.seed(2017-12-16)
library(raster)
#> Loading required package: sp
elev1 = raster(nrow = 6, ncol = 6, res = 0.5, 
               xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
               vals = 1:36)
r_mask1 = raster(nrow = 6, ncol = 6, res = 0.5, 
              xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
              vals = sample(c(NA, TRUE), 36, replace = TRUE))
elev1[r_mask1, drop = FALSE]
#> class      : RasterLayer 
#> dimensions : 6, 6, 36  (nrow, ncol, ncell)
#> resolution : 0.5, 0.5  (x, y)
#> extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> source     : memory
#> names      : layer 
#> values     : 1, 36  (min, max)


set.seed(2017-12-16)
library(terra)
#> terra version 1.4.2
elev2 = rast(nrow = 6, ncol = 6, resolution = 0.5, 
             xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
             vals = 1:36)
r_mask2 = rast(nrow = 6, ncol = 6, resolution = 0.5, 
               xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
               vals = sample(c(NA, TRUE), 36, replace = TRUE))
elev2[r_mask2, drop = FALSE]
#>    lyr.1
#> 1      1
#> 2      2
#> 3      3
#> 4      4
#> 5      5
#> 6      6
#> 7     13
#> 8     14
#> 9     15
#> 10    16
#> 11    17
#> 12    18

^{Created on 2021-09-12 by the reprex package (v2.0.1)}

[rhijmans]

Thank you. I now get:

set.seed(2017-12-16)
library(terra)
#> terra version 1.4.2
elev2 = rast(nrow = 6, ncol = 6, resolution = 0.5, xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,           vals = 1:36)
clip2 = rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45, resolution = 0.3, vals = rep(1, 9))
elev2[clip2]
#     lyr.1
#[1,]    18
#[2,]    24

elev2[clip2, drop=FALSE]
#class       : SpatRaster 
#dimensions  : 2, 1, 1  (nrow, ncol, nlyr)
#resolution  : 0.5, 0.5  (x, y)
#extent      : 1, 1.5, -0.5, 0.5  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#name        : lyr.1 
#min value   :    18 
#max value   :    24 

[Nowosad]

@rhijmans I found another two cases when drop = FALSE does not work:

set.seed(2017-12-16)
library(terra)
#> terra version 1.4.2
elev = rast(nrow = 6, ncol = 6, resolution = 0.5, 
             xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
             vals = 1:36)

elev[1:2, drop = FALSE]
#>   lyr.1
#> 1     1
#> 2     2

elev[1, 1:2, drop = FALSE]
#>   lyr.1
#> 1     1
#> 2     2

^{Created on 2021-09-13 by the reprex package (v2.0.1)}

[rhijmans]

I could say that it does work:

library(terra)
#terra version 1.4.2
elev = rast(nrow = 6, ncol = 6, resolution = 0.5,   xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,  vals = 1:36)
elev[1:2, drop = TRUE]
#[1] 1 2
elev[1:2, drop = FALSE]
#  lyr.1
#1     1
#2     2

But not as you expected. And what you expected is more interesting. Now I get:

elev[1:2, drop = FALSE]
#class       : SpatRaster 
#dimensions  : 1, 2, 1  (nrow, ncol, nlyr)
#resolution  : 0.5, 0.5  (x, y)
#extent      : -1.5, -0.5, 1, 1.5  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#name        : lyr.1 
#min value   :     1 
#max value   :     2 

And (selecting rows instead of cells)

elev[1:2,  , drop = FALSE]
#class       : SpatRaster 
#dimensions  : 2, 6, 1  (nrow, ncol, nlyr)
#resolution  : 0.5, 0.5  (x, y)
#extent      : -1.5, 1.5, 0.5, 1.5  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#name        : lyr.1 
#min value   :     1 
#max value   :    12 

Selecting columns

elev[, 1:2,drop = FALSE]
#class       : SpatRaster 
#dimensions  : 6, 2, 1  (nrow, ncol, nlyr)
#resolution  : 0.5, 0.5  (x, y)
#extent      : -1.5, -0.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#name        : lyr.1 
#min value   :     1 
#max value   :    32 

Selecting rows and columns

elev[1:2, 1:2,drop = FALSE]
#class       : SpatRaster 
#dimensions  : 2, 2, 1  (nrow, ncol, nlyr)
#resolution  : 0.5, 0.5  (x, y)
#extent      : -1.5, -0.5, 0.5, 1.5  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#name        : lyr.1 
#min value   :     1 
#max value   :     8 

[Nowosad]

@rhijmans - yes, you are right. Results were different from my expectations. I also have another example of that - raster result is different from the terra one:

# raster ------------------------------------------------------------------
set.seed(2017-12-16)
library(raster)
#> Loading required package: sp
elev1 = raster(nrow = 6, ncol = 6, res = 0.5, 
               xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
               vals = 1:36)
r_mask1 = raster(nrow = 6, ncol = 6, res = 0.5, 
                 xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
                 vals = sample(c(NA, TRUE), 36, replace = TRUE))
r1 = elev1[r_mask1, drop = FALSE]
plot(r1)

# terra -------------------------------------------------------------------
set.seed(2017-12-16)
library(terra)
#> terra version 1.4.2
elev2 = rast(nrow = 6, ncol = 6, resolution = 0.5, 
            xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
            vals = 1:36)
r_mask2 = rast(nrow = 6, ncol = 6, resolution = 0.5, 
              xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
              vals = sample(c(NA, TRUE), 36, replace = TRUE))
r2 = elev2[r_mask2, drop = FALSE]
plot(r2)

^{Created on 2021-09-14 by the reprex package (v2.0.1)}

[rhijmans]

Thank you that is very helpful. terra now should return the same results as raster in this case. That is, the SpatRaster will be cropped and masked by the index SpatRaster. To only use crop, you can do rc = elev2[ext(r_mask2), drop = FALSE]

[rhijmans]

Oh, wait, but what if the rasters do not align? In that case, it should be crop only, I suppose?

[Nowosad]

I would think that the elev2[r_mask2, drop = FALSE] syntax works identically to mask(crop(elev2, r_mask2), r_mask2)).

[ailich]

@rhijmans @Nowosad, this does not seem entirely fixed yet. While you can replace values based on an index raster with a single value, in terra you can't update the values in a way that depends on the values themselves (e.g. for cells that meet this condition multiply them by 2), while this is possible in raster. I've included an example below and you can see that raster provides the expected results while terra outputs a warning and a different result.

library(raster)
library(terra)


#raster
r1<- raster(nrow=3, ncol=3)
values(r1)<- 1:9
r1[r1>5]<- r1[r1>5] * 2
plot(r1)
text(r1)

#terra
r2<- rast(nrow=3, ncol=3)
values(r2)<- 1:9
r2[r2>5]<- r2[r2>5] * 2
# Warning message:
#   In v[as.logical(values(i))] <- value :
#   number of items to replace is not a multiple of replacement length
plot(r2)
text(r2)

[rhijmans]

I made some changes and now both Alex's and Jakub's example have the same results with raster and terra. I would not be surprised if I broke something else though.

[rhijmans]

library(raster)
library(terra)

r1<- raster(nrow=3, ncol=3)
values(r1)<- 1:9
r1[r1>5]<- r1[r1>5] * 2

r2<- rast(nrow=3, ncol=3)
values(r2)<- 1:9
x <- ifel(r2>5, r2*2, r2)
r2[r2>5]<- r2[r2>5] * 2

all(values(r2) == values(r1))
#[1] TRUE
all(values(r2) == values(x))
#[1] TRUE 

Also to point out the alternative route with ifel

[ailich]

Awesome, thanks!

[ailich]

@rhijmans, it is not quite NA proof yet.

library(terra)
r2<- rast(nrow=3, ncol=3)
values(r2)<- 1:9
r2[5]<-NA
r2[r2>5]<- r2[r2>5] * 2
#Error in if (length(value) == sum(i)) { : 
#  missing value where TRUE/FALSE needed

[rhijmans]

Thank you very much Alex. I believe that this last and all the prior examples now work in terra in the same way as in raster. This last one is tricky, it is not clear to me whether it should work as your indexing with NAs, but if the replacement is with a data.frame (that is what r2[r2>5] * 2 evaluates too) I will allow it.

I would reiterate that this type of notation is nice for quick interactive use, but not for production scripts, and certainly not for functions in a package. They are fragile and not they are not memory safe (will fail with large files). I would use ifel or the what is beneath that (cover, classify, mask).

[ailich]

@rhijmans, thanks this works now. I did not realize however that this method is not memory safe. I'll update code in my own functions to avoid this syntax. I have also used base::ifelse and dplyr::case_when within a call to app in some situations. Would that be memory safe? For example:

library(terra)
r2<- rast(nrow=3, ncol=3)
values(r2)<- 1:9
r2[5]<-NA
r2<- app(r2, fun = function(x)(ifelse(x>5, x*2,x)))

r3<- rast(nrow=3, ncol=3)
values(r3)<- 1:9
r3[5]<-NA
#case_when cleaner than ifelse if many different conditions
r3<- app(r3, fun = function(x)(dplyr::case_when(x > 5 ~ x*2,
                                                is.na(x) ~ 0,
                                                x < 3 ~ x*3,
                                                TRUE ~ x)))

[rhijmans]

Yes, that should be memory safe. In the first case, ifel(r2 > 5, r2*2, r2) might be more efficient. In the second case you would need a nested ifel and that would probably make it less efficient.

[ailich]

Thanks!