## Rasmus Benestad ## R-script for analysis. ## R is freely available from http://cran.r-project.org ## It's a 'free lunch' :-) ## The esd-library can be installed by un-commenting the two following lines: ##library(devtools) ##install_github('metno/esd') ## Alternatively, it can be retieved from: https://github.com/metno/esd ## The news article addressed with this R-script ##http://www.telegraph.co.uk/news/earth/environment/globalwarming/11395516/The-fiddling-with-temperature-data-is-the-biggest-science-scandal-ever.html ## ## There is also the comment from the Guardian: ##http://www.theguardian.com/commentisfree/2008/sep/23/controversiesinscience.health ## Function to read the GISS data stored as ASCII and convert into suitable structure: readgiss <- function(fname,path='Dropbox/R/') { data <- read.table(paste(path,fname,sep=''),header=TRUE) y <- as.station.data.frame(data,loc=substr(fname,6,nchar(fname)-4), param='t2m',unit='deg C',lon=NA,lat=NA,alt=NA, cntr=NA,longname=NA,stid=NA,quality=NA,src=NA,url=NA, reference=NA,info=NA, method= NA) x <- coredata(y) x[abs(x) > 99] <- NA coredata(y) <- c(x) invisible(y) } ## Activate the esd-package: library(esd) ## Retrieve the Nordli et al. (2014) Svalbard temperature (air port/lufthavn) data(Svalbard) ## The GISS data are retrieved from http://data.giss.nasa.gov/gistemp/station_data/ ## Date of retrieval is 2015-02-10 giss <- list.files(path='Dropbox/R',pattern='GISS-') ## Get the GISS equivalent of the Svalbard temperature ## (here the index is 18, but this may differ according to the number of GISS sites ## retrieved) y.giss <- readgiss(giss[18]) ## Only use matching times with valid data in both records. y <- matchdate(Svalbard,y.giss) miss <- !is.finite(coredata(y)) & !is.finite(coredata(y.giss)) y[miss] <- NA; y.giss[miss] <- NA ## Make the first figure plot(annual(y),errorbar=FALSE,lwd=5) grid() lines(annual(y.giss),lwd=3) lines(trend(annual(y)),col='red',lty=2) lines(trend(annual(y.giss)),lty=2) legend(2000,-7,c('Nordli et a. (2014)','GISS'), col=c('red','black'),lwd=c(5,3),bty='n',cex=0.75) dev2bitmap('Svalbard-Nordli+giss.png') ## Get the NACD data - stored in theesd-package nacd <- station(src='nacd') ## Loop - compare annual trends between corresponding sites: X <- rep(NA,length(giss)); Y <- X; cex=X; col=rgb(rep(0,length(giss)),rep(0,length(giss)),rep(1,length(giss)),0.2) for (i in 1:length(giss)) { x <- readgiss(giss[i]) j <- grep(tolower(loc(x)),tolower(loc(nacd))) if (length(j)>0) { y <- subset(nacd,is=j) y <- matchdate(y,x) x <- matchdate(x,y) miss <- !is.finite(coredata(y)) & !is.finite(coredata(x)) y[miss] <- NA; x[miss] <- NA x <- annual(x); y <- annual(y) plot.zoo(x,y) X[i] <- trend.coef(zoo(x)) Y[i] <- trend.coef(zoo(y)) cex[i] <- 0.05*length(zoo(x)) col[i] <- rgb(0,0,1,0.3) } } ## Figure 2: jpeg(filename = "trends-Nordli+giss.png",width=600,height=600,units="px",bg = "white") par(bty='n') plot(c(-0.5,1.2),c(-0.5,1.2), xlab=expression(paste('GISS ',(degree*C/decade))) ,ylab=expression(paste('NACD ',(degree*C/decade))), type='l',lwd=5,col="grey",main='Comparison of annual trends') points(X,Y,pch=19,cex=cex,col=col) grid() dev.off()