# class 11/12  examples for correlation
# loader function for Albany, NY Central Park
# station data (monthly mean anomalies)
source("scripts/loadano.R")



# Albany station data: USW00014735

month<-"Ann"
x1<-loadano(station="USW00014735",month=month,start=1950,end=2012)



# New York Central Park station data: USW00094728
x2<-loadano(station="USW00094728",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)




# in-class development:
# calculate and tabulate for each month:
# mean, standard deviations, covariance, and correlation


# 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))



# last step: adjust the title option in the plot-functions