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,] target3 <- Tann # 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)) # SOLVE BEST LINE AND SINE WAVE FOR 20thC TOGETHER #s <- cbind(0,0,0,0,0,10000) #for (A in seq(11,16,1)) { #for (b in seq(25,35,1)) { #for (T in seq(60,65,1)) { #for (m in seq(40,50,1)) { #for (a in seq(94,101,1)) { # #sina <- (A/100) * sin(((x-b)/T)*(2*pi)) #lina <- x*m/10000 - a/10 #res1 <- target[,2] - (lina + sina) #ss <- sum(res1^2) #s <- rbind(s,cbind(A,b,T,m,a,ss)) #}}}}} #s[s[,6]==min(s[,6],na.rm=T),] #A <- s[s[,6]==min(s[,6]),1] #b <- s[s[,6]==min(s[,6]),2] #T <- s[s[,6]==min(s[,6]),3] #m <- s[s[,6]==min(s[,6]),4] #a <- s[s[,6]==min(s[,6]),5] #A #b #T #m #a A = 9 # GIS b = 33 # GIS T = 61 # GIS m = 46 # GIS a = 90 # GIS A = 12 # CRU b = 32 # CRU T = 61 # CRU m = 49 # CRU a = 97 # CRU sina <- (A/100) * sin(((x-b)/T)*(2*pi)) lina <- x*m/10000 - a/10 # file name to write plot CairoPNG(paste("img/tt-d-both-line-and-sine-",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 Simultaneous Line and Sine Fit\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() res1 <- target[,2] - (lina + sina) s2 <- sum(res1^2) s2 # 0.889 cor(target[,2], lina+sina) # 0.886 # add leading and trailing CairoPNG(paste("img/tt-d-both-line-sine-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 Simultaneous Line and Sine Fit\n1880-2000", sub=paste("http://rhinohide.wordpress.com ",dstamp), cex.lab=0.75) linb <- target2[,1]*m/10000 - a/10 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 res1 <- target2[,2] - (linb + sinb) s2 <- sum(res1^2) s2 # 0.889 cor(target2[,2], linb+sinb) # 0.886 # add leading and trailing CairoPNG(paste("img/tt-d-both-line-sine-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 Simultaneous Line and Sine Fit\n1880-2010", sub=paste("http://rhinohide.wordpress.com ",dstamp), cex.lab=0.75) linb <- Tann[,1]*m/10000 - a/10 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 res1 <- Tann[,2] - (linb + sinb) s2 <- sum(res1^2) s2 # 0.889 cor(Tann[,2], linb+sinb) # 0.886 s <- NULL #for (A1 in seq(12,16,1)) { #for (b1 in seq(1,5,1)) { #for (T1 in seq(65,75,1)) { #for (A2 in seq(20,25,1)) { #for (b2 in seq(120,130,1)) { #for (T2 in seq(135,145,1)) { for (A1 in seq(90,100,2)) { for (b1 in seq(0,5,1)) { for (T1 in seq(1,11,2)) { for (A2 in seq(65,85,2)) { for (b2 in seq(110,130,2)) { for (T2 in seq(120,135,2)) { sin1 <- (A1/100) * sin(((target[,2]-b1)/T1)*(2*pi)) sin2 <- (A2/100) * sin(((target[,2]-b2)/T2)*(2*pi)) res2 <- target[,2] - (sin1 + sin2) ss <- sum(res2^2) s <- rbind(s,cbind(A1,b1,T1,A2,b2,T2,ss)) }}}}}} s[s[,7]==min(s[,7],na.rm=T),] #A1 <- s[s[,7]==min(s[,7]),1] #b1 <- s[s[,7]==min(s[,7]),2] #T1 <- s[s[,7]==min(s[,7]),3] #A2 <- s[s[,7]==min(s[,7]),4] #b2 <- s[s[,7]==min(s[,7]),5] #T2 <- s[s[,7]==min(s[,7]),6] quit() A1=15# GIS b1=2# GIS T1=71# GIS A2=23# GIS b2=123# GIS T2=141# GIS A1=15# CRU b1=2# CRU T1=71# CRU A2=23# CRU b2=123# CRU T2=141# CRU sin1 <- (A1/100) * sin(((x-b1)/T1)*(2*pi)) sin2 <- (A2/100) * sin(((x-b2)/T2)*(2*pi)) # file name to write plot CairoPNG(paste("img/tt-d-both-sine-",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 Simultaneous Two Sine Fit\n1900-2000", sub=paste("http://rhinohide.wordpress.com ",dstamp), cex.lab=0.75) lines(target[,1],sin1+sin2) # residuals for annual global T points(target[,1],target[,2]) # residuals for annual global T dev.off() res1 <- target[,2] - (sin1 + sin2) s2 <- sum(res1^2) s2 # 0.920 cor(target[,2], sin1+sin2) # 0.885 # add leading and trailing CairoPNG(paste("img/tt-d-both-sine-",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 Simultaneous Two Sine Fit\n1880-2000", sub=paste("http://rhinohide.wordpress.com ",dstamp), cex.lab=0.75) sin1a <- (A1/100) * sin(((target2[,1]-b1)/T1)*(2*pi)) sin2a <- (A2/100) * sin(((target2[,1]-b2)/T2)*(2*pi)) lines(target2[,1],sin1a+sin2a,col="red") # residuals for annual global T points(target2[,1],target2[,2],col="green") # residuals for annual global T lines(target[,1],sin1+sin2) # residuals for annual global T points(target[,1],target[,2]) # residuals for annual global T res1 <- target2[,2] - (sin1a + sin2a) s2 <- sum(res1^2) s2 # 0.920 cor(target2[,2], sin1a+sin2a) # 0.885 # add leading and trailing CairoPNG(paste("img/tt-d-both-sine-",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 Simultaneous Two Sine Fit\n1880-2010", sub=paste("http://rhinohide.wordpress.com ",dstamp), cex.lab=0.75) sin1a <- (A1/100) * sin(((Tann[,1]-b1)/T1)*(2*pi)) sin2a <- (A2/100) * sin(((Tann[,1]-b2)/T2)*(2*pi)) lines(Tann[,1],sin1a+sin2a,col="red") # residuals for annual global T points(Tann[,1],Tann[,2],col="green") # residuals for annual global T lines(target[,1],sin1+sin2) # residuals for annual global T points(target[,1],target[,2]) # residuals for annual global T res1 <- Tann[,2] - (sin1a + sin2a) s2 <- sum(res1^2) s2 # 0.920 cor(Tann[,2], sin1a+sin2a) # 0.885