# This is a script to detrend and deseasonalize mm SPEI 
#===================================================================================================================
library(zoo)
library(reshape2)
#===================================================================================================================

spei.ind = c('01','02','03')
years=c(2000:2013)

for (n in 1:length(spei.ind)){
  print(spei.ind[n])
  
  #---Load gridded SPEI
  load(paste('./data/LOCA.obs/spei',spei.ind[n],'_obs_195001-201312.Rdata',sep=''))
  ss=res
  itime = which(ss$time >= "1999-12-01" & ss$time <= "2013-11-01")
  lon=ss$lon; lat=ss$lat
  df = ss$data[,,itime]
  
  #--- Remove linear trend
  res1 = array(NA, dim=dim(df))
  time=c(1:dim(df)[3])
  
  for (i in 1:length(lon)){
    for (j in 1:length(lat)){
      x = df[i,j,]
      x.cmpl = x[complete.cases(x)]
      if (length(x.cmpl) >= 2){ 
        fit = lm(df[i,j,] ~ time, na.action=na.exclude)
        res1[i,j,] = resid(fit) # extract residuals 
      }
    }}
  
  #--- Remove seasonality
  temp = array(res1, dim=c(length(lon), length(lat), 12, dim(df)[3]/12)) # lon x lat x month x year
  mm = apply(temp, c(1,2,3), FUN=mean, na.rm=TRUE)
  
  new.res1 = array(res1, dim=dim(temp))
  res2 = array(NA, dim=dim(temp))
  
  for (i in 1:length(lon)){
    for (j in 1:length(lat)){
      for(k in 1:length(month.abb)){
        res2[i,j,k,] = new.res1[i,j,k,] - mm[i,j,k]
      }}}
  
  iav = array(res2, dim=dim(res1))
  ss= list('iav'=iav, 'lon'=lon, 'lat'=lat, 'time'=ss$time[itime])
  save(ss, file=paste('./data/LOCA.obs/dtds/dtds.mm-spei',spei.ind[n],'.2000-2013.RData',sep=''))
}

#===================================================================================================================
#===================================================================================================================