############################################################################################################## ## ## This is R code to compare Hadley with CMIP5 mean ## ############################################################################################################## rm(list=ls(all=TRUE)) # remove all objects in memory require(doBy) # will use to collapse data loc="http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.2.0.0.monthly_ns_avg.txt" ht1 = read.table(file = loc, dec=".") ht=as.matrix(ht1) h = ht[,2] h = as.numeric(h) z = strsplit(ht[,1],"/") # split date string into year, month zm = matrix(unlist(z), ncol=2, byrow=TRUE) # convert list into 2-col matrix year = as.numeric(zm[,1]) # year as number endyear=max(year) month = as.numeric(zm[,2]) # month as number rm(ht1, ht, z, zm) h = h[year>=1861] y = year[year>=1861] loc="data/knmi-model-mean.txt" # Downloaded from http://climexp.knmi.nl/ ht1 = read.table(file = loc, dec=".", skip=4) g = as.matrix(ht1[,2]) dat = data.frame(y,h,g) z = summaryBy(g+h ~ y, data=dat ) # collapse by year, month (mean) y = as.numeric(z[,1]) # year h = as.numeric(z[,3]) # HadCRUT4 h = h - mean(h[100:130]) # Centred on 1961-1990 g = as.numeric(z[,2]) # GCM Mean g = g - mean(g[100:130]) # Centred on 1961-1990 dat = data.frame(y,h,g) write.table(dat, file = "comparison_data.csv", sep=",") # Write data to CSV file # Plot the 2 series together par(mar=c(6,6,2,1) ) plot(h ~ y, ylab = "deg C anomaly", xlab = "Year", ylim = c(-1, 1), xlim = c(1850,2020), type="l", lwd=2, col="black" ) abline(h=0, lty=2) par(new=T) plot(g ~ y, axes=F, type="l", lty=1, xlab="", ylab="", ylim = c(-1, 1), xlim = c(1850,2020), lwd=6, col="gray")