Reputation: 21
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
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