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,]

# 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 <- 2000 # end of trend
t0 <- 100 # years of trend
y0 <- y1 - t0 # start of trend
t1 <- t0 # number of years to display 
y2 <- y0 + t1 # max of x-axis
ymin <- -0.5 # min of y-axis
ymax <- 0.8 # max of y-axis
dstamp=Sys.Date()

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

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

#--------------------------------
# GISS 1880 - 2010 annual
#--------------------------------

# 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)
t_trend <- Tann[Tann[,1]>=1880 & Tann[,2]<=2009,]

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

# fit a simple exponential
# initial value from my hand fits earlier
nls(y ~ -beta1/100 + beta2/10 * exp((beta3/10000)*(x-1880)), start=list(beta1=44,beta2=1,beta3=210), trace=T)

# GIS: 1.458176 :   32.7136099   0.5939041 209.0974014
# CRU: 2.008889 :   44.3962016   0.6017519 210.9072272
# NOA: 1.602463 :   26.8185131   0.4356877 231.8945787 

if (idx == gis) {
    a <- 32.71
    b <- 0.59
    k <- 209.10
}
if (idx == cru) {
    a <- 44.40
    b <- 0.61
    k <- 210.91
}
if (idx == noa) {
    a <- 26.82
    b <- 0.44
    k <- 231.89
}
a
b
k
y1 = -a/100 + b/10 * exp((k/10000)*(x-1880))

# get correlation
cor(y1,y) # GIS = 0.910 # HADCRU = 0.889 # NCDC = 0.901

# plot fit
CairoPNG(paste("img/nls-",src[idx-2],"-1880-2010.png", sep=""))
plot(x,y1,type="l", ylab="", xlab="")
points(x,y)
title(paste(src[idx-2],"-1880-2010.png", sep=""))
dev.off()

# plot projection
xx <- 1880:2100
yy = -a/100 + b/10 * exp((k/10000)*(xx-1880))

CairoPNG(paste("img/nls-",src[idx-2],"-1880-2100.png", sep=""))
plot(xx,yy,type="l")
points(x,y)
title(paste(src[idx-2],"-1880-2100.png", sep=""))
dev.off()