jeudi 7 février 2019

Double nested loop with embedded "if/else" statement to iterate over multidimensional array in R

I have 4 dimensional array (ensemble of climate models with dimensions: lon, lat, time, models). I want to iterate over models (4th dim), and apply a function on each decade/10 year chunk (sequence over 3rd dim). The entire length of time series in my real data is 27 years, so I needed to introduce the "if" statement for the last decade which does not offer complete 10 years. However, when I try to manually verify the results - I am getting slightly higher numbers. This does not happen when I remove the if statement or use just perfect 20 years time series. It seems to me that the "if" statement somehow also changes the indices and counters for full decades. Can, please, someone shed the light what is happening? I hope that corrected code maybe useful for many people working on multidimensional arrays in general. I prepared the simplified code, which also generates the data and reproduce this behavior:

array1<- array(1:120, dim=c(3,3,17)) 
array2<- array(500:620, dim=c(3,3,17))
array3<- array(1000:1120, dim=c(3,3,17))

# create a empty array
nlon<-3
nlat<-3
nt<-17
last_decade_length <- 7

play_array<-matrix(0,nlon*nlat*nt*3)
dim(play_array)<-c(nlon,nlat,nt,3)

#Allocate DATA of array1 (model one)
play_array[,,,1]<-array1
#Allocate DATA of array2 (model two)
play_array[,,,2]<-array2
#Allocate DATA of array3 (model three)
play_array[,,,3]<-array3

# create an object to hold final result
mean_decade_object<-matrix(0,2*3)
dim(mean_decade_object)<-c(2,3)

### NESTED DOUBLE LOOP - "a" is index to itirate over decades; "b" index to itirate over arrays/models
a<-1
b<-1

for (array in 1:dim(play_array)[4]) {
  print(paste("array",array))
  for (decade in seq(1,dim(play_array)[3], 10)){ 
    if(length(decade)==10){
      mean_decade_object[a,b] <-mean(play_array[,,decade:(decade+9), array], na.rm = T) 
      print(mean(play_array[,,decade:(decade+9), array], na.rm = T))
      print(paste("decade",decade))
      a <-a+1
    }
    else{
      mean_decade_object[a,b] <-mean(play_array[,,decade:(decade+last_decade_length-1), array], na.rm = T) 
      print(mean(play_array[,,decade:(decade+last_decade_length-1), array], na.rm = T))
      print(paste("decade",decade))
      a <-a+1       
    } 
  }
  a<-1
  b<-b+1
}  
print(mean_decade_object)

##############
#VERIFICATION#
##############

sub1 <- play_array[,,,1]
mean(sub1[,,1:10], na.rm=T)`

Why is the result of verification (mean of 1st decade of 1st model) is not matching the top left cell in object resulting from the loop? Any ideas?

Aucun commentaire:

Enregistrer un commentaire