####################################################
# 2014-03-27 updated to work with data in
# New York state sub-folder data/NY 
# 2014-02-04 calculate seasonal climatological cycle
#
# data from
# Global Historical Climatology Network -Daily              
# Daily data were retrieved from 
# http://www.ncdc.noaa.gov/oa/climate/ghcn-daily/
# For this course I preprocessed the data 
# into monthly mean data
# ---------------------------------------
# Data table description:
# column 1: year
# column 2: month
# column 3: days in month
# column 4: index counter 
# column 5: monthly mean of the variable
# column 6: number of observations
# column 7: number of missing observations
# column 8: ratio number of observations
#                        / days in month
# If too few observations were made
# I set the monthly mean to 'NA'
# NA is a special value in R indicating
# no value. 
###################################################
# clear memory
rm(list=ls())


###################################################
# User's program control parameters
###################################################
# station
# 2014-03-27 added the subdir (state code)
state<-"NY"

stationname<-"USW00094789"
# variable
varname<-"tavg"
# first and last year for averaging
yr1<-1981
yr2<-2010

###################################################
# if isAno set TRUE the program also writes the 
# monthl mean anomalies to a CSV-file
# (see lines 120-131)
###################################################
isAno<-TRUE


###################################################
# Main program 
###################################################
# 2014-03-27 updated the file name
# (state subdir added)
infile<-paste("data/",state,"/",stationname,"_",varname,"_mon_mean.asc",sep='')
outfile<-paste("data/",state,"/",stationname,"_",varname,"_mon_mean_climc_",yr1,"-",yr2,".csv",sep='')


print ("calculate monthly mean climatology")
print(paste("Station :",stationname))
print(paste("variable :",varname))
buffer<-read.table(infile)

# Creating a time coordinate
inum<-length(buffer[,1])
tstart<-buffer[1,1]
dt<-1.0/12
ihelp<-seq(1,inum)
time<-tstart+(ihelp-1)*dt


#calculate the mean seasonal cycle climatology
#Using the 1981-2010 climate normal period (NOAA definition)

isclim<-(time>=yr1 & time<=(yr2+11/12))
thelp<-time[isclim]
# x contains the monthly mean data
x<-buffer[isclim,5]
# index for the month in the year (1 for Jan, 2 for Feb ..., 12 for Dec)
mhelp<-buffer[isclim,2]
inum2<-length(thelp) 
# xmean is a vector of length 12
# for the caculation of 30-year means
# in each month Jan, Feb, ... ,Dec
xmean<-array(0,dim=c(12))
nsum<-array(0,dim=c(12))
for (i in 1:inum2) { 
  m<-mhelp[i]
  # now we must check at position i 
  # buffer has a  real number or NA
  # is.na() is a builtin funtion 
  # that checks for NA and returns TRUE or FALSE
  if (is.na(x[i]) ) {
    # skip
  } else {
    # real value
    nsum[m]<-nsum[m]+1
    xmean[m]<-xmean[m]+x[i]
  }
}
for (m in 1:12) {
  if (nsum[m]>25) {
    xmean[m]<-xmean[m]/nsum[m]
  } else {
    xmean[m]<-NA
  }
}



##############################################
# write the climatology 
# results to a file (CSV table)
##############################################
label<-c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
output<-list(month=label,mean=xmean,n=nsum)
write.csv(output,file=outfile,row.names=FALSE,quote=FALSE)
print(paste("wrote climatology to file ",outfile))


##############################################################################
# note you can load the result with 
# clim<-read.csv("[REPLACE WITH FILENAME]")
# try then: 
# barplot(clim$mean,names.arg=clim$month,ylab="[REPLACE WITH PROPER LABEL]")
##############################################################################



##############################################################################
# THIS PART WRITES ALSO THE ANOMALIES (deviations from the
# 30-yr seasonal mean cycle)
##############################################################################
if (isAno) {
  source("scripts/myfunctions.R")
  # xmean
  x<-buffer[,5] # monthly mean data of entire period
  res<-anomaly(x,xmean)
  output<-list(time=time,anomaly=res)
  # 2014-03-27 updated to work with state subfolder
  outfile2<-paste("data/",state,"/",stationname,"_",varname,"_mon_mean_ano.csv",sep='')
  write.csv(output,file=outfile2,row.names=FALSE,quote=FALSE)
  print(paste("wrote anomalies to file ",outfile2))
}