1
votes

I want to modify a subset of layers in a raster brick by multiplying those layers by another raster.

For example, if we have a raster brick called 'r.brick' and we try to multiply its layers 2:4 by a raster, 'r.mult', with same row&column dimensions:

  • r.brick[[2:4]] returns layers 2:4, as expected

  • r.brick[[2:4]] * r.mult successfully multiplies those layers, as expected

but, if I try to assign the result back into the subset layers I get an error

r.brick[[2:4]] = r.brick[[2:4]] * r.mult
# Error in value[] <- val : 
#   incompatible types (from S4 to double) in subassignment type fix

The error message suggests that the assignment is trying to assign the raster values rather than the raster itself. But if I try the assignment with getValues, I get a different error:

r.brick[[2:4]] = getValues(r.brick[[2:4]] * r.mult)
# Error in .local(x, values, ...) : length(values) is not equal to ncell(x)

What is the correct way to do this?

Some reproducible data:

library(raster)
r.list = vector("list", 20)
set.seed(123)
for (i in 1:20) {
  r.list[[i]] = raster(vals=runif(100), nrows=10, ncols=10, ext=extent(c(0,25,0,25)))
}
r.brick = brick(r.list)
r.mult = raster(vals=sample(2,100,T), nrows=10, ncols=10, ext=extent(c(0,25,0,25)))
3
Did you try using setValues ? (Don't know if it allows to assign on multiple layers, though...) - lbusett
@LorenzoBusetto setValues only seems to work on single layers, so far as I can tell. - dww

3 Answers

3
votes

Then I only got a workaround (using a loop) for you:

layers <- 2:4
for(i in layers) {
  r.brick[[i]]  <- r.brick[[i]] * r.mult
}

Note: Apparently the assignment with the [] subsetting only works with single layers.

2
votes

I think that is a missing feature. Thanks for pointing that out. I would go with the loop following maRtin (perhaps after creating a RasterStack, that may be more efficient). If the dataset is not too large, you can do

# example data
library(raster)
b <- brick(nrows=2, ncols=2, nl=6)
values(b) <- rep(1:4, 6) 
r.mult <- raster(vals=10, nrows=2, ncols=2)


values(b)[,3:4] <- values(b[[3:4]] * r.mult)
# values(b)
1
votes

The following code with multiple layers 2-4 of r.brick raster values with r.mult values and assign the results to r.brick layers 2-4.

attr(attr(r.brick, 'data'), 'values')[,2:4] = attr(attr(r.brick, 'data'), 'values')[,2:4] * attr(attr(r.mult, 'data'), 'values')

or

attr(attr(r.brick, 'data'), 'values')[,2:4] = getValues(r.brick[[2:4]] * r.mult)