library(Cairo) #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 <- 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", "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,] target2 <- Tann[Tann[,1]<= y1,] # a series for x-axis tyear <- 1880:2010 # some set setup x <- target[,1] y <- target[,2] new <- data.frame(x = seq(y0+1,y1)) # BEST SINE WAVE FOR 20thC s <- cbind(0,0,0,10000) for (A in 5:20) { for (b in 25:29) { for (T in 30:60) { sina <- (A/100) * sin(((x-b)/T)*(2*pi)) res2 <- target[,2] - sina ss <- sum(res2^2) s <- rbind(s,cbind(A,b,T,ss)) }}} s[s[,4]==min(s[,4]),] A <- s[s[,4]==min(s[,4]),1] b <- s[s[,4]==min(s[,4]),2] T <- s[s[,4]==min(s[,4]),3] A b T #A = 16 #b = 28 #T = 50 sina <- (A/100) * sin(((x-b)/T)*(2*pi)) res2 <- target[,2] - sina s2 <- sum(res2^2) s2 # 2.7 cor(target[,2],sina) # 0.578 # file name to write plot CairoPNG(paste("img/tt-b-sine-",src[idx-2],"-",y1,"-",t0,".png",sep=""), bg="white") plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") title(main="20th Century Sine Fit", sub=paste("http://rhinohide.wordpress.com ",dstamp), cex.lab=0.75) lines(target[,1],sina) # residuals for annual global T points(target[,1],target[,2]) # residuals for annual global T dev.off() # BEST LINEAR TREND # get trend; hi and lo for 95% ci trend_lm <- lm(res2 ~ x) lina <- x*trend_lm$coeff[2] + trend_lm$coeff[1] trend_lm$coeff[2] # 0.43 res1 <- res2 - lina s1 <- sum(res1^2) s1 # 1.1 cor(res2,lina) # 0.76 # file name to write plot CairoPNG(paste("img/tt-b-sine-with-line-residuals-",src[idx-2],"-",y1,"-",t0,".png",sep=""), bg="white") # plot # hidden - it just sets the graph axis plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") title(main="20th Century Sine with Linear Residual Trend", sub=paste("http://rhinohide.wordpress.com ",dstamp), cex.lab=0.75) # points in the trend lines(target[,1],lina) # data points for annual global T points(target[,1],res2) # data points for annual global T dev.off() # COMBINE LINEAR AND SINE CairoPNG(paste("img/tt-b-sine-line-",src[idx-2],"-1900-2000.png",sep=""), bg="white") plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") title(main="20th Century Sine Fit with Linear Trend of Residuals\n1900-2000", sub=paste("http://rhinohide.wordpress.com ",dstamp), cex.lab=0.75) lines(target[,1],lina + sina) # residuals for annual global T points(target[,1],target[,2]) # residuals for annual global T dev.off() res3 <- target[,2]-(lina + sina) s3 <- sum(res3^2) s3 # 1.11 cor(target[,2],lina+sina) # 0.858 # add leading and trailing CairoPNG(paste("img/tt-b-sine-line-all-",src[idx-2],"-1880-2000.png",sep=""), bg="white") plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") title(main="20th Century Sine Fit with Linear Trend of Residuals\n1880-2000", sub=paste("http://rhinohide.wordpress.com ",dstamp), cex.lab=0.75) linb <- target2[,1]*trend_lm$coef[2] + trend_lm$coef[1] sinb <- (A/100) * sin(((target2[,1]-b)/T)*(2*pi)) lines(target2[,1],linb+sinb,col="red") # residuals for annual global T points(target2[,1],target2[,2],col="green") # residuals for annual global T lines(target[,1],lina + sina) # residuals for annual global T points(target[,1],target[,2]) # residuals for annual global T res3 <- target2[,2]-(linb + sinb) s3 <- sum(res3^2) s3 # 1.11 cor(target2[,2],linb+sinb) # 0.858 # add leading and trailing CairoPNG(paste("img/tt-b-sine-line-all-",src[idx-2],"-1880-2010.png",sep=""), bg="white") plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Year", ylab="Anomaly (deg C)") title(main="20th Century Sine Fit with Linear Trend of Residuals\n1880-2010", sub=paste("http://rhinohide.wordpress.com ",dstamp), cex.lab=0.75) linb <- Tann[,1]*trend_lm$coef[2] + trend_lm$coef[1] sinb <- (A/100) * sin(((Tann[,1]-b)/T)*(2*pi)) lines(Tann[,1],linb+sinb,col="red") # residuals for annual global T points(Tann[,1],Tann[,2],col="green") # residuals for annual global T lines(target[,1],lina + sina) # residuals for annual global T points(target[,1],target[,2]) # residuals for annual global T res3 <- Tann[,2]-(linb + sinb) s3 <- sum(res3^2) s3 # 1.11 cor(Tann[,2],linb+sinb) # 0.858 # clean up #rm(list=ls(all=TRUE))