0

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
1

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
  • 2
    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
0

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 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
# [39] 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[[4]])
#          [,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.

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.