library(Cairo)

#t <- read.table("http://processtrends.com/Files/RClimate_consol_temp_anom_latest.csv", sep=",", header=T)
t <- read.table("../stats/RClimate_consol_temp_anom_latest.csv", sep=",", header=T)

# trim to 1880
t <- t[t[,1]>=1880,]

# column idx
gis <- 3
cru <- 4
noa <- 5
rss <- 6
uah <- 7
idx <- cru

src = c("GISTEMP", "HADCRU", "NOAA", "RSS", "UAH")

# annualize monthly data
t$year <- as.integer(t[,1])

# get the record of interest, monthly to annual
Tann <- aggregate(t[,idx],list(t$year), FUN=mean, na.rm=T)
target <- Tann[Tann[,1]>=1880 & Tann[,1]<= 2009,]

# some set setup
x <- target[,1]
y <- target[,2]

# monthly
#x <- t[!is.na(t[,idx]),1]
#y <- t[!is.na(t[,idx]),idx]

b = -0.53
m = 0.0059
A = 0.30
T = 60
y1 = (m*(x-1880) + b) + (A * cos(((x-1880)/T)*(2*pi)))
xx = 1880:2100
yy = (m*(xx-1880) + b) + (A * cos(((xx-1880)/T)*(2*pi)))
ym = (m*(xx-1880) + b)

# monthly data has lower correlation
cor(y1,y) # CRU:1880:2009 = 0.872 and CRU:1880:2010 = 0.874
# correlation improves to 0.891 with A = 0.25

# file name to write plot
CairoPNG(paste("img/girma-",src[idx-2],"-1880-2009.png",sep=""), bg="white", width=552, height=336)

# plot
tyear = 1880:2100
# hidden - it just sets the graph axis
plot(tyear,rep(0,length(tyear)), ylim=c(-1,+1), col="white", xlab="Year", ylab="Anomaly (deg C)", axes=FALSE) 
title(main=paste("20th Century",src[idx-2],"Line + Sine\ny1 = (0.0059*YEAR - 0.53) + (0.3 * cos(((YEAR)/60)*(2*pi)))"), sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.sub=0.75)

# points in the trend
axis(1, at=seq(1880,2100,by=10), cex.axis=0.7)
axis(2, at=seq(-1,1,by=0.2))
grid(24,12,col="lightgrey")
lines(xx,yy) # model to 2100
lines(xx,ym) # trend line
zero = rep(0, length(xx))
lines(xx,zero) # zero line
lines(x,y, col="red") # annual anomaly
lines(x,y1,col="blue") # model to 2009
dev.off()

# clean up
#rm(list=ls(all=TRUE))