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

Question

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)

Answer 1

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

Answer 2

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