# This is a script to calculate detrended and deseasonalized monthly mean fine dust anomalies 
#===================================================================================================================
library(zoo)
library(reshape2)
#===================================================================================================================

#---Load site info
load('./data/IMPROVE/IMPROVE_Southwest_siteinfo.RData') # df
site.info=df

#---Load fine dust data and subset 2000-2013
load('/groups/anenberggrp/PA/02.dust/IMPROVE/raw/IMPROVE_2000-2016-mm_noncd-Fe.RData')
head(ss)
years=c(2000:2013)
smon=month.abb[1:12]

df = ss[ss$site %in% site.info$SiteCode & ss$year %in% years, c(1,2,3,5)]
df$ym = (paste(df$mon, df$year, sep =" ")) 
df$date = as.Date(as.yearmon(df$ym))
df2 = df[,c(1,6,4)]
df3 = df2[order(df2$site, df2$date),]

#--- Reshape into time x sites array, and screen out sites with < 50% complete data
temp = acast(df3, date~site, value.var="FD") 
ts=rownames(temp)
min.na = dim(temp)[1]/2
na.count <- sapply(temp, function(x) sum(length(which(is.na(x)))))
na.count = array(na.count, dim=dim(temp))
na.sum = apply(na.count, 2, FUN=sum)
ind=which(na.sum > min.na)
ss = temp[,-ind] # 35 sites

#--- Calculate 2000-2013 annual and regional mean
x = array(ss, dim=c(12,14,35))
x2 = apply(x, c(2,3), FUN=mean, na.rm=TRUE) # annual-mean
x3 = apply(x2, 1, FUN=mean, na.rm=TRUE) # regional-mean
mean(x3) # 1.03 ug m-3

### FOR EACH SITE:
#--- (1) Remove linear trend
sites = colnames(ss)
time=c(1:dim(ss)[1])
res1 = array(NA, dim=dim(ss))

for (i in 1:length(sites)){
    fit = lm(ss[,i] ~ time, na.action=na.exclude)
    res1[,i] = resid(fit) # extract residuals 
}

#--- (2) Remove seasonality (subtract long-term monthly means and divide by s.d. of res1)
temp = array(res1, dim=c(12, dim(ss)[1]/12, dim(ss)[2])) # month x year x site
mm = apply(temp, c(1,3), FUN=mean, na.rm=TRUE)
sd = apply(temp, c(1,3), FUN=sd, na.rm=TRUE)

new.res1 = array(res1, dim=dim(temp))
res2 = array(NA, dim=dim(temp))
res3 = array(NA, dim=dim(temp))

for (i in 1:length(sites)){
  for(j in 1:length(month.abb)){
  res2[j,,i] = (new.res1[j,,i] - mm[j,i])/sd[j,i]
  res3[j,,i] = (new.res1[j,,i] - mm[j,i]) # non-standardized
}}

temp = array(res2, dim=dim(ss))
iav = rbind(rep(NA, length=length(sites)), temp[-dim(temp)[1],]) # Dec 1999 - Nov 2013
colnames(iav)=sites
save(iav, file='./data/IMPROVE/EOF/dtds+st.mm-fd.2000-2013.Rdata') # metyear

temp = array(res3, dim=dim(ss))
iav = rbind(rep(NA, length=length(sites)), temp[-dim(temp)[1],]) # Dec 1999 - Nov 2013
colnames(iav)=sites
save(iav, file='./data/IMPROVE/EOF/dtds.mm-fd.2000-2013.Rdata') # metyear

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

#--- Save out info of final 35 sites
ss = site.info[site.info$SiteCode %in% sites, c(1:6)]
write.csv(ss, file='../outputs/IMPROVE.sites.csv')

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