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
}
4
  • 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, 2018 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. May 3, 2018 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, 2018 at 12:03
  • Yes, you were correct, I 've also changed that in the original post! May 3, 2018 at 12:58

2 Answers 2

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
4
  • 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, 2018 at 12:05
  • @Spacedman thanks for your suggestion! I appreciate it
    – aldo_tapia
    May 3, 2018 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. May 3, 2018 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, 2018 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 and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.