Selecting layers of a raster stack using a sliding windows as a step in R in for loop

I would like to create a loop that every time selects the first ten layers of a raster stack, calculates their mean and saves it to a layer of a new raster stack. Also, I would like to have a step that works like a sliding window. So, in the first loop, it will select layers from 1:10, in the second loop it will select layers from 11:20, and so on. So far, I have written the following code:

file_list=stack(myfile.nc)
mn = stack()
#
for(i in seq(1,length(file_list@layers), by=10)){
mn[[i]] <- mean[[i]] # supposed to calculate the mean of the first ten layers
}
• What does that code do? Does it run? What is mean[[i]]? Did you mean to write mean(somethign)? What's file_list? Does that seq call work? – Spacedman May 3 '18 at 11:53
• @Spacedman, I have edited the question above. The seq call gives layers 1,10,20,30 etc, while I would to get layers 1-10,11-20,21-30 etc instead. – Maria Karypidou May 3 '18 at 12:01
• No, that seq call won't work. Have you tried it? For example, if the length is 100, then you get seq(1:100, by=10) which is an error. You probably meant seq(1,100,by=10), but please test code that you post. – Spacedman May 3 '18 at 12:03
• Yes, you were correct, I 've also changed that in the original post! – Maria Karypidou May 3 '18 at 12:58

I agree with @Spacedman commentaries, I recommend you to read some basic tutorials to know how to handle spatial objects in R. But I have an example of what you want to do.

Reproducible example, a raster stack of 100 layers:

library(raster)

set.seed(123)

r <- raster()

s1 <- list()

for (i in 1:100) {
s1[[i]] <- setValues(r, rnorm(n = ncell(r)))
}

s1 <- stack(s1)

Compute mean by a moving window:

s2 <- list()

for (i in 1:10) {
s2[[i]] <- calc(s1[[((i-1)*10 + 1):((i-1)*10 + 10)]], fun = mean, na.rm = T)
}

s2 <- stack(s2)

Result:

s2

## class       : RasterStack
## dimensions  : 180, 360, 64800, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 1, 1  (x, y)
## extent      : -180, 180, -90, 90  (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,   layer.5,   layer.6,   layer.7,   layer.8,   layer.9,  layer.10
## min values  : -1.353835, -1.314173, -1.336899, -1.330095, -1.572422, -1.276887, -1.413552, -1.292723, -1.293121, -1.366087
## max values  :  1.247697,  1.380478,  1.324504,  1.305732,  1.395083,  1.518274,  1.324168,  1.340525,  1.288073,  1.247176
• Here's a quicker way of making a random stack of 3x4 by 50 layers s = stack(brick(array(runif(3*4*50),dim=c(3,4,50)))) – Spacedman May 3 '18 at 12:05
• @Spacedman thanks for your suggestion! I appreciate it – aldo_tapia May 3 '18 at 12:07
• Dear aldo_tapia, thank you for your answer but what i actually mean by the "moving window", is that I would like to have a moving step in the for loop. So, in the first iteration I will have the mean of the first ten layers of the raster, in the second iteration I will have the following ten layers and so on. – Maria Karypidou May 3 '18 at 12:55
• Yes, actually is exactly what you mean. Test this for (i in 1:10) print(((i-1)*10 + 1):((i-1)*10 + 10)) – aldo_tapia May 3 '18 at 13:08

This can be done using zApply.

Create a test array of 50 3x4 layers:

a = array(runif(3*4*50),dim=c(3,4,50))

Convert to a Raster stack object (have to go via brick):

s = stack(brick(a))

The next step needs the object to have Z values, but they can be anything:

s = setZ(s, 1:nlayers(s))

Now make a vector of group indices:

layers = rep(1:5,rep(10,5))
layers
#  1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4
#  4 4 5 5 5 5 5 5 5 5 5 5

Then apply mean over the layers grouped by the layer indices:

s10mean = zApply(s, by=layers,fun=mean)

This is now a brick with five layers with the means:

s10mean
# class       : RasterBrick
# dimensions  : 3, 4, 12, 5  (nrow, ncol, ncell, nlayers)

To check, first lets see what the mean of the 31st to 40th parts of the original array are:

apply(a[,,31:40],c(1,2),mean)
#          [,1]      [,2]      [,3]      [,4]
#[1,] 0.5128400 0.4061487 0.5870473 0.2036074
#[2,] 0.4610763 0.4732387 0.4928084 0.4972154
#[3,] 0.6670825 0.5265644 0.5172894 0.4435808

and that should be the same as the values in the fourth layer of the output:

as.matrix(s10mean[])
#          [,1]      [,2]      [,3]      [,4]
#[1,] 0.5128400 0.4061487 0.5870473 0.2036074
#[2,] 0.4610763 0.4732387 0.4928084 0.4972154
#[3,] 0.6670825 0.5265644 0.5172894 0.4435808

Looks good.