# class 15 (derived from class12.R)
# 2014-03-27 changes introduced to work
#           with station data in subfolders
#           for each state (NY state currently)

# updated the function that loads the 
# monthly mean anomaly data
source("scripts/loadano2.R")


month<-"Jan"
# Station 1

x1<-loadano2(station="USW00014750",state="NY",month=month,start=1950,end=2012)



# Station 2
x2<-loadano2(station="USW00094728",state="NY",month=month,start=1950,end=2012)




# create a 1x2 panel of figures on one page
# to reset the plot window to a 1-panel plot use:
# par(mfrow(c(1,1)))

par(mfrow=c(1,2))
# left: time series Albany  and NY Central Park January
plot(c(1950,2012),c(-8,8),typ='n',xlab='year',ylab='Tavg [C]',main=month)
lines(x1$time,x1$ano,
     typ='b',pch=1,col=4)
lines(x2$time,x2$ano,
      typ='b',pch=3,col=2)

# right: scatter plot Albany January with NY Central Park 
plot(c(-8,8),c(-8,8),typ='n',xlab='Alb. Tavg [C] ',ylab='NYC Tavg [C]',main=month)
points(x1$ano,x2$ano,pch=3)




#2014-03-27 [NOTE: VALUES ARE NOT CALCULATED ]
#            IF TIME SERIES ARE INCOMPLETE TIME]

# calculate the mean of all years
mx1<-mean(x1$ano)
mx2<-mean(x2$ano)

print(paste("mean of data set 1 :",mx1))
print(paste("mean of data set 2 :",mx2))


# calculate the standard deviation of all years
sdx1<-sd(x1$ano)
sdx2<-sd(x2$ano)
print(paste("standard deviation of data set 1:",sdx1))
print(paste("standard deviation of data set 2:",sdx2))

#covariance:
n<-length(x1$ano)
covx1x2<-1/n*sum((x1$ano-mx1)*(x2$ano-mx2))
print(paste("covariance between the two data sets:",covx1x2))
# correlation:
corx1x2<-covx1x2/(sdx1*sdx2)
print(paste("Pearson product-moment correlation :",corx1x2))



# 2014-03-31 added the calculation of the regression line

bfit<-corx1x2*sdx2/sdx1
afit<-mx2-bfit*mx1

yfit<-afit+bfit*x1$ano

lines(x1$ano,yfit,col=5,lwd=2)

print("regression line: station 2 regressed on station 1")

print(paste("intercept: ",afit,sep=''))
print(paste("slope: ",bfit,sep=''))













###################################################
# 2014-03-31 adding for demonstration purpose only
###################################################
if (FALSE){
  print("regression line 2: station 2 regressed on station 2")
  bfit2<-corx1x2*sdx1/sdx2
  afit2<-mx2-bfit2*mx2
  xfit2<-afit2+bfit2*x2$ano
  lines(xfit2,x2$ano,col=6,lwd=1)
}