Improbability
Improbability

Reputation: 21

R: Replacing several layers in a raster stack with layers from a second raster stack

I have a raster stack and need to replace some of the layers it contains with layers from another stack.

I guess the solution is pretty simple and straightforward but so far I was unable to find it.

To give you a reproducible example, here is what I have tried to do:

library(raster)

# create some raster
r1 <- raster(nrows=10, ncols=10, xmn=0, xmx=10, ymn = 0, ymx = 10)
r1[] <- sample(1:100, 100, replace = TRUE)

r2 <- raster(nrows=10, ncols=10, xmn=0, xmx=10, ymn = 0, ymx = 10)
r2[] <- sample(1:100, 100, replace = TRUE)

r3 <- raster(nrows=10, ncols=10, xmn=0, xmx=10, ymn = 0, ymx = 10)
r3[] <- sample(1:100, 100, replace = TRUE)

r4<- raster(nrows=10, ncols=10, xmn=0, xmx=10, ymn = 0, ymx = 10)
r4[] <- sample(1:100, 100, replace = TRUE)

# put the raster into a stack
r_stack <- stack(r1, r2, r3, r4)

#calculate the mean of the raster
r_mean <- mean(r_stack)

# What I would like to do is to subtract the mean from some of the
# raster in the stack (layers_to_replace), but not from all, and to             
# replace the raster in the stack with the new difference.
# For example I would like to replace the second and fourth layer with     
# the difference.
l_replace <- c(2, 4)

# Note: I place the difference into a second stack for the sake
# of the example as my original data comes in a second stack
rep_stack <- r_stack[[l_replace]] - r_mean

r_stack[[l_replace]] <- rep_stack

Unfortunately this approach does not work and throws the following error:

Error in v[] <- value : incompatible types (from S4 to logical) in subassignment type fix
In addition: Warning messages:
1: In if (i < 1) { : the condition has length > 1 and only the first element will be used
2: In if (i > nl + 1) { : the condition has length > 1 and only the first element will be used
3: In if (i > nl) { : the condition has length > 1 and only the first element will be used

Any ideas on how to solve this issue will be more than welcome.

Upvotes: 1

Views: 599

Answers (1)

lbusett
lbusett

Reputation: 5932

This should work:

for(ind in seq(along = l_replace)) r_stack[[l_replace[ind]]][] <- rep_stack[[ind]][]

r_stack
class       : RasterStack 
dimensions  : 10, 10, 100, 4  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 0, 10, 0, 10  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
names       : layer.1, layer.2, layer.3, layer.4 
min values  :    1.00,  -54.25,    2.00,  -54.25 
max values  :   100.0,    51.5,   100.0,    51.5

Alternatively, for somehow easier readability:

for(ind in seq(along = l_replace)){
  setValues(r_stack, getValues(rep_stack[[ind]]), layer = l_replace[ind])
} 

Upvotes: 0

Related Questions