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
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