Recategorise a terra raster based on 2 and 3 other layers with conditions

Question

I have 2 rasters (x1 and x2) that I use to recategorize a new layer 'x' in 4 different categories :

• if values in x1 is 0 AND values in x2 is 0, then category 1

• if values in x1 is 1 AND values in x2 is 0, then category 2

• if values in x1 is 0 AND values in x2 is 1, then category 3

• if values in x1 is 1 AND values in x2 is 1, then category 4

  set.seed(1234)
  
  # Set 2 rasters 
  x1 <- rast(xmin=-110, xmax=-80, ymin=40, ymax=70, ncols=30, nrows=30)
  x2 <- rast(xmin=-110, xmax=-80, ymin=40, ymax=70, ncols=30, nrows=30)
  
  # With different values 
  values(x1) <- sample(c(0,1), size = ncell(x1), replace = T)
  values(x2) <- sample(c(0,1), size = ncell(x2), replace = T)
  
  # Show them 
  par(mfrow = c(1,2))
  plot(x1)
  plot(x2)
  
  # Make a new 'base' raster that will be filled with 4 different categories
  values(x) <- 0
  par(mfrow = c(1,1))
  plot(x)
  
  # Set 4 different categories
  x[x1 == 0 & x2 == 0] <- 1
  x[x1 == 1 & x2 == 0] <- 2
  x[x1 == 0 & x2 == 1] <- 3
  x[x1 == 1 & x2 == 1] <- 4
  plot(x)
  

This works, but as stated here, this is not advised for larger raster. In addition, I'd like to generalize and have more categories (i.e., if I have a 3rd layer 'x3', then I'd have to generate all tree-way categories (expand.grid(c(0,1), c(0,1), c(0,1))).

The only way I found with 2 layers to make it easy is to use 0-1 values for the first layer and 0-50 values for the second such that :

x2.2 <- rast(xmin=-110, xmax=-80, ymin=40, ymax=70, ncols=30, nrows=30)

values(x2.2) <- sample(c(0,50), size = ncell(x2.2), replace = T)
plot(x2.2)

new_x = sum(x1, x2.2)
plot(new_x)
new_cat = rbind(c(0, 1),
                c(1, 2),
                c(50, 3),
                c(51, 4))
z <- classify(new_x, rcl = new_cat)
plot(z)

For 3 layers, this would be an example :

x3 <- rast(xmin=-110, xmax=-80, ymin=40, ymax=70, ncols=30, nrows=30)
values(x3) <- sample(c(0,150), size = ncell(x3), replace = T)

new_x_all = sum(x1, x2.2, x3)
plot(new_x_all)

Is there a better way of doing this?

eg = expand.grid(c(0,1), c(0,50), c(0,200))
cbind(eg, apply(eg, 1, sum))
  Var1 Var2 Var3 apply(eg, 1, sum)
1    0    0    0                 0
2    1    0    0                 1
3    0   50    0                50
4    1   50    0                51
5    0    0  200               200
6    1    0  200               201
7    0   50  200               250
8    1   50  200               251

Answer 1

You can use terra::unique with as.raster=TRUE

u <- unique(c(x1, x2), as.raster=TRUE)

u
#class       : SpatRaster 
#size        : 30, 30, 1  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : -110, -80, 40, 70  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source(s)   : memory
#categories  : label, lyr.1, lyr.1.1 
#name        : label 
#min value   :   0_0 
#max value   :   1_1 

levels(u)
#[[1]]
#  ID label
#1  0   0_0
#2  1   0_1
#3  2   1_0
#4  3   1_1

This also works for more layers, and if numbers are not restricted to 0 and 1.

For the 0/1 case, it is possible to use one's base2 approach (see below) but it is probably less efficient. Here is how I would implement that approach using (memory safe) terra idiom:

v <- app(c(x1, x2), \(i) strtoi(paste(i, collapse=""), 2))

Answer 2

If we only works with 0-1 values, we could collapse all rasters' value and treat it as base2 (binary) number. Then we can convert it to base10.

eg = expand.grid(c(0,1), c(0,1), c(0,1))
apply(eg,1,\(x) strtoi(paste(x,collapse = ""),base=2))
[1] 0 4 2 6 1 5 3 7

Note that it is not necessary to generate all n-way possibilities since there is 1-1 relationship between base2 and base10.

One way to use this approach is:

library(terra)

set.seed(1234)

# Set rasters 
x <- rast(xmin=-110, xmax=-80, ymin=40, ymax=70, ncols=30, nrows=30)
x1 <- rast(xmin=-110, xmax=-80, ymin=40, ymax=70, ncols=30, nrows=30)
x2 <- rast(xmin=-110, xmax=-80, ymin=40, ymax=70, ncols=30, nrows=30)
x3 <- rast(xmin=-110, xmax=-80, ymin=40, ymax=70, ncols=30, nrows=30)


values(x1) <- sample(c(0,1), size = ncell(x1), replace = T)
values(x2) <- sample(c(0,1), size = ncell(x2), replace = T)
values(x3) <- sample(c(0,1), size = ncell(x2), replace = T)


# use base conversion
values(x)<-strtoi(paste(values(x1),values(x2),values(x3),sep=""),base=2)
head(values(x))

     lyr.1
[1,]     4
[2,]     4
[3,]     5
[4,]     4
[5,]     2
[6,]     6