0

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)
2
  • Why don't you just use focal? It looks like you are emulating a focal function anyway. Nov 29 '21 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.
    – user44796
    Nov 29 '21 at 17:05
2

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) )
2
  • 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!).
    – user44796
    Nov 29 '21 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. Nov 29 '21 at 21:53
1

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
2
  • 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. Nov 30 '21 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!) Nov 30 '21 at 1:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.