# This is a script to derive linear sensitivity of fine dust to SPEI02 anomalies (regional means)
#===================================================================================================================
library(boot)
source('./my.functions.R')
#===================================================================================================================

years=c(2000:2013)
seas=c('DJF','MAM','JJA','SON')

#--- Load fine dust anomalies and calculate regional means

load('./data/IMPROVE/EOF/dtds.mm-fd.2000-2013.Rdata') # iav
sites=colnames(iav)
temp=array(iav, dim=c(12,dim(iav)[1]/12,length(sites)))
temp2=array(temp, dim=c(3,4,length(years),length(sites))) # month x seas x year x site
df=apply(temp2, c(2,3,4), FUN=mean, na.rm=TRUE) # seasmean
fd = apply(df, c(1,2), FUN=mean, na.rm=TRUE) # regional mean
colnames(fd)=years
rownames(fd)=seas

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

#---Load SPEI anomalies and calculate regional means for different seasons
spei = array(NA, dim=c(length(seas),length(years)))

for (k in 1:length(seas)){
  print(seas[k])
  files=c('spei02','spei02','spei02','spei02')
  load(paste('./data/LOCA.obs/dtds/dtds.mm-',files[k],'.2000-2013.RData',sep=''))

  if(k==1){xlim=c(-117,-105); ylim=c(30,38)}
  if(k==2){xlim=c(-115,-103); ylim=c(28,38)}
  if(k==3){xlim=c(-115,-103); ylim=c(29,40)}
  if(k==4){xlim=c(-113,-103); ylim=c(30,39)}

  ilon=which(ss$lon >= xlim[1] & ss$lon <= xlim[2])
  ilat=which(ss$lat >= ylim[1] & ss$lat <= ylim[2])

  df = ss$iav[ilon,ilat,]
  df2 = array(df, dim=c(length(ilon), length(ilat), 12, dim(df)[3]/12))
  df3 = df2[,,(3*k-2):(3*k),]
  df4 = apply(df3, c(1,2,4), FUN=mean, na.rm=TRUE) # seasmean
  spei[k,] = apply(df4, 3, FUN=function(data) calc_awa(data, lon=ss$lon[ilon], lat=ss$lat[ilat]) ) # regmean
}

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

#--- Perform SLR

bs <- function(formula, data, indices) {
  d <- data[indices,] # allows boot to select sample 
  fit <- lm(formula, data=d)
  return(coef(fit)) 
} 

#layout(matrix(c(1:4), 2, 2, byrow = TRUE)) # layout.show(nf)

for(k in 1:4){ # seas
  print(seas[k])

  # (1) Check XY scatterplot  
  y = fd[k,]
  x = spei[k,]
  #plot(x,y)  
  #print(cor.test(x,y))

# (2) Regression analysis
  df=cbind(y, x)
  df=as.data.frame(df)
  colnames(df)=c('fd','spei')
  final.lm = lm(df$fd~df$spei)
  print(summary(final.lm)) # 
  y.hat = fitted(final.lm)
  
  # bootstrapping with 1000 replications 
  results <- boot(data=df, statistic=bs, R=10000, formula=fd~spei)
  temp=(boot.ci(results, type="bca", index=2)) # 95% CI slope
    
  #---Save results
  ss = list('data'=df, 'SLR'=summary(final.lm), 'yhat'=y.hat,  
            'beta'=final.lm$coefficients[2], 'beta.CI'=c(temp$bca[4], temp$bca[5]))
  save(ss, file=paste('./data/slr/FD/spei02/slr_2000-2013-',seas[k],'.RData',sep='')) 
  }

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