R version 2.12.1 (2010-12-16)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[Previously saved workspace restored]

> 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
[1] 0.8720133
> # 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)
Fontconfig error: Cannot load default config file
> 
> # 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()
null device 
          1 
> 
> # clean up
> #rm(list=ls(all=TRUE))
> 
> 
> 
> proc.time()
   user  system elapsed 
  0.560   0.020   0.577