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

# note: this does not yet handle NAs in the trend well
# be careful with the bounderies - 1979 for sats, pre-1900 for ground
# need boundary tests

y1 <- 2007 # end of trend
t0 <- 25 # years of trend
y0 <- y1 - t0 # start of trend
t1 <- 3 * t0 # number of years to display 
y2 <- y0 + t1 # max of x-axis
ymin <- -0.5 # min of y-axis
ymax <- 2.5 # max of y-axis

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

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]> y0 & Tann[,1]<= y1,]

# a series for x-axis
tyear <- y0:y2

# some set setup
x <- target[,1]
y <- target[,2]
new <- data.frame(x = seq(y0+1,y1))

# get trend; hi and lo for 95% ci
trend_lm <- lm(y ~ x)
trend_lo = cbind(x,predict(lm(y ~ x), interval="prediction")[,2])
trend_hi = cbind(x,predict(lm(y ~ x), interval="prediction")[,3])
lm_lo <- lm(trend_lo[,2] ~ trend_lo[,1])
lm_hi <- lm(trend_hi[,2] ~ trend_hi[,1])

# extend trend +/- ci into future years
trend_lm_lo <- tyear*lm_lo$coef[2] + lm_lo$coef[1]
trend_lm_hi <- tyear*lm_hi$coef[2] + lm_hi$coef[1]
trend_lm_mid <- tyear*trend_lm$coef[2] + trend_lm$coef[1]

# set y axis
ymax = max(trend_lm_hi) + 0.5
ymin = min(trend_lm_lo) - 0.5

# file name to write plot
png(paste("trendtest-",src[idx-2],"-",y1,"-",t0,".png",sep=""),height=360, width=405, bg="white")

# plot
dstamp=Sys.Date()
plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") # hidden - it just sets the graph axis
title(main=paste("Did Global Warming Stop in the Year",y1,"?\n",t0,"Year", src[idx-2],"Trend Tester"),
	sub=paste("http://rhinohide.wordpress.com       ",dstamp), cex.lab=0.75)

# warm, cold, no trend boxes
lines(c(y1,y1),c(ymin,ymax), type="l", col="gray") # beginning year of trend test
rect(y1,trend_lm_hi[t0+1],y2,ymax, col="pink")
rect(y1,trend_lm_lo[t0+1],y2,ymin, col="light blue")
vtx <- min(min(tyear[trend_lm_lo > trend_lm_hi[t0+1]]), y2)
polygon(rbind(c(y1,trend_lm_lo[t0+1]),c(vtx,trend_lm_lo[vtx-y0+1]),c(y2,trend_lm_hi[t0+1]),c(y2,trend_lm_lo[t0+1])), col="light gray", lty="solid")

# label trend boxes
text(y2,(trend_lm_lo[t0+1]+trend_lm_hi[t0+1])/2, "No Trend", pos=2, cex=0.8) 
text(y1+1,(ymin+trend_lm_lo[t0+1])/2, "Cooling", pos=4, cex=0.8) 
text(y1+1,(ymax+trend_lm_hi[t0+1])/2, "Warming", pos=4, cex=0.8) 
text(y0+1,(ymax+trend_lm_hi[t0+1])/3, paste(t0,"year trend\n95% confidence interval"), pos=4, cex=0.8) 

# trend slope data
a1 <- format(trend_lm$coeff[2]*10,digits=2)
a2 <- format(trend_lm_hi[t0+1] - trend_lm_mid[t0+1],digits=2)
a3 <- format(-trend_lm_lo[t0+1] + trend_lm_mid[t0+1],digits=2)
text(y0+1,ymin+0.1, paste(a1,"C per decade\n+",a2,"/ -",a3), pos=4, cex=0.8) 

# points in the trend
points(target[,1],target[,2]) # data points for annual global T

# trend lines and +/- confidence interval
abline(trend_lm,col="black")  # trendline for annual global T
lines(trend_lm_hi ~ tyear,col="red")  # 2 SD upper limit of monthly global T
lines(trend_lm_lo ~ tyear,col="blue") # 2 SD lower limit of monthly global T

# add the trend test points
t_test <- Tann[Tann[,1]> y1,]
points(t_test[,1],t_test[,2], col="green") # data points beyond the trend setters
dev.off()

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