Multiply raster brick by vector based on codition from another vector in R

I am trying to multiply a raster brick by a vector based on the conditions of a second raster. The question is similar to one I posted here, but is a little more complicated. Here is an example.

Suppose I have a brick `x` and I want to multiply it by a vector `trans` that has values in three classes: `1`, `2`, and `3`. I have a second matrix, `y`, that has values for the class for each cell. I know that I could do this by coding the `trans` vector in the function and using `overlay`, but I have reasons for wanting to keep the trans vector outside the function. Here is what I have so far, but realize that the f_trans function isn't correct.

``````trans = data.frame(class = c("1","1","1","2","2","2","3","3","3"),value = c(.33,.66,.99,.25,.5,0.75,.1,.4,.9))

mat = matrix(data = rep(1,9),nrow=3,ncol=3)
x = brick(raster(mat),raster(mat),raster(mat))

class_matt = matrix(data = c(1,1,1,2,2,2,3,3,3),nrow=3,ncol=3)
class_rast = raster(class_matt)

f_trans = function(x,y,trans){
z = rep(NA, length(x))
i <- which(y == 1)
z[i] <- x[i] * trans\$value[trans\$class==1]
i <- which(y == 2)
z[i] <- x[i] * trans\$value[trans\$class==2]
i <- which(y >= 3)
z[i] <- x[i] * trans\$value[trans\$class==3]

out = raster(data = z, ncol=ncol(x),nrow = nrox(x))
return(out)
}

results = f_trans(a,class_rast,trans)
``````
• Why don't you just use focal? It looks like you are emulating a focal function anyway. Commented Nov 29, 2021 at 16:54
• I don't think it is a focal type operation, unless I am mistaken. I must not have explained the question well. What I am trying to do it multiply the layers at cell (x,y) by the values in trans, conditional on the values in class_rast. Commented Nov 29, 2021 at 17:05

First, lets create your example data. Note that in your example the "a" object is not created so, I am assuming that it is your raster stack. I am stacking the conditional raster with the "a" raster stack to simplify the function and the object passed to `raster::overlay`.

``````library(raster)
s = data.frame(class = c("1","1","1","2","2","2","3","3","3"),
value = c(.33,.66,.99,.25,.5,0.75,.1,.4,.9))
cr = raster(matrix(c(1,1,1,2,2,2,3,3,3),nrow=3,ncol=3))
a <- do.call(stack, replicate(3,
raster(matrix(sample(1:9, 9),
nrow=3, ncol=3))))
a <- stack(cr, a)
``````

The first thing to note is that, in the "trans" (now "s") data.frame object, you have replicated values for each conditional "class". As such, in your function `trans\$value[trans\$class==1]` will return three scalars and not one, which is what is currently indexed. If these represent weights for each class this is not a bad thing as long as they are ordered ie., in class 1, c(0.33, 0.66, 0.99) represents the scalars for the associated raster layers 1:3.

In a function like `raster::overlay` values are recycled so, if your function writes out 3 values then they will be assigned to three separate rasters in a stack object. We can rely on the fact that in the function x will always have the same number of values as layers in the stack, even if they are NA's.

If this is the case then your function can be written as such. Please note that due to limitations in passing additional function arguments to `raster::overlay`, I am calling your scalar data.frame object from the global environment and not as an argument. Your function will need to reflect this explicit object naming condition.

``````f_trans = function(y){
if(y[[1]] == 1) {
x = y[-1] * s[s\$class==1,]\$value
} else if(y[[1]] == 2) {
x = y[-1] * s[s\$class==2,]\$value
} else if(y[[1]] >= 3) {
x = y[-1] * s[s\$class==3,]\$value
} else {
x = rep(NA,length(y[-1])) }
return(x)
}
``````

Now, lets test the function for each condition. We can see that it returns the expected scalar values. Now we can be confident that the scalars match the vector that they are being applied to.

``````  f_trans(c(1, rep(1,3)))
f_trans(c(2, rep(1,3)))
f_trans(c(3, rep(1,3)))

``````

The function is now appropriate to pass to a function like `raster::overlay`.

``````( r <- overlay(a, fun=f_trans) )

# Check first value
( x <- as.numeric(as.matrix(a)[1,]) )
f_trans(x)
as.matrix(r)[1,]
``````

Note that you can also do this in terra, which is the replacement for the raster package, and it will be much faster in practice. Here is the same analysis in terra. Note that we can now have a function argument for our scalar object that is then specified in `terra::app`.

``````library(terra)
sl = data.frame(class = c("1","1","1","2","2","2","3","3","3"),
value = c(.33,.66,.99,.25,.5,0.75,.1,.4,.9))
a <- c(rast(matrix(c(1,1,1,2,2,2,3,3,3),nrow=3,ncol=3)),
do.call(c, replicate(3, rast(matrix(sample(1:9, 9),
nrow=3, ncol=3)))))

f_trans = function(y, s){
if(y[[1]] == 1) {
x = y[-1] * s[s\$class==1,]\$value
} else if(y[[1]] == 2) {
x = y[-1] * s[s\$class==2,]\$value
} else if(y[[1]] >= 3) {
x = y[-1] * s[s\$class==3,]\$value
} else {
x = rep(NA,length(y[-1])) }
return(x)
}
( r <- app(a, fun=f_trans, s=sl) )
``````
• That works, and let's me keep the scalar parameters as global. I wasn't even thinking about globals in R (which are useful and violate so many rules!). Commented Nov 29, 2021 at 21:33
• Oh so many rules. The problem with `raster::overlay` is that to pass additional arguments you have to write a function that wraps your function and, you still cant do it directly in `raster::overlay`. However, I would point out that `terra::app` has a ellipse `...` argument that does allow for passing additional arguments to a custom function, it just cannot be another raster object. Commented Nov 29, 2021 at 21:53

Here are three alternative approaches to `terra::app` that Jeffrey Evans proposed; instead using `ifel`, `selectRange` and `lapp`. I am not saying these are better in any way, but they may be easier to use.

Jeffrey's example data with x and y separate

``````library(terra)
set.seed(1)
x <- rast(nrow=3, ncol=3, nlyr=3, xmin=0)
values(x) <- rep(1:3, each=9)
y <- rast(nrow=3, ncol=3, nlyr=1, xmin=0, vals=rep(1:3, 3))
z <- data.frame(class = c(1,1,1,2,2,2,3,3,3),
value = c(1,2,3,4,5,6,7,8,9))
``````

The `ifel` approach

``````r1 <- ifel(y==1, x*z[z\$class==1,]\$value,
ifel(y==2, x*z[z\$class==2,]\$value,
ifel(y==3, x*z[z\$class==3,]\$value, NA)))
``````

The `selectRange` approach

``````# all products (9 layers, a bit inefficient)
xx <- x * z\$value
#The groups start at layers 1, 4, and 7
yy <- y+(y-1) * 2
r2 <- selectRange(xx, yy, 3)
``````

`lapp` is good for this, but the tricky thing for this problem is that you need to transpose the data in the function to get the right multiplications.

``````ftr <- function(x, y, s){
x <- t(x)
for (i in 1:3) {
j <- y==i
x[,j] = x[,j] * s\$value[s\$class==i]
}
t(x)
}

r3 <- lapp(sds(x, y), ftr, z)
``````

For comparison, Jeffrey's approach with `app`

``````f_trans = function(y, s){
if(y[[1]] == 1) {
x = y[-1] * s[s\$class==1,]\$value
} else if(y[[1]] == 2) {
x = y[-1] * s[s\$class==2,]\$value
} else if(y[[1]] >= 3) {
x = y[-1] * s[s\$class==3,]\$value
} else {
x = rep(NA,length(y[-1])) }
return(x)
}
r0 <- app(c(y,x), fun=f_trans, s=z)

all(values(r0) == values(r1))
#[1] TRUE

all(values(r0) == values(r2))
#[1] TRUE

all(values(r0) == values(r3))
[1] TRUE
``````
• So in this case `terra::ifel` produces multiple values? Using `base::ifelse` will only return the first element in the scaled values because the input condition is a single value. You can observe this behavior here: `y=1; x=c(.33,.66)` this syntax will return a single value `ifelse(y == 1, 1 * x)` but this one returns the correct length `ifelse(rep(y,2) == 1, 1 * x)`. I had assumed that `terra::ifel` would behave in a similar way so, good to know that it does not. BTW, I have been very happy with the availability of `ifel` and use it frequently. Commented Nov 30, 2021 at 0:47
• I suppose a `SpatRaster` can be considered a single value, irrespective of the number of layers. If you use a vector with numbers in `ifel`, only the first one is used. This computation `x*z[z\$class==1,]\$value` returns a single `SpatRaster`. (And glad you like it!) Commented Nov 30, 2021 at 1:00