##BRIFFA ET AL 2001 FIGURE 4 ###PREREQUISITES ##DOWNLOAD crutem3 version from cru and install in directory d:/climate/data/gridcell/crutem3 ## install R package ncdf ##LOAD BRIFFA 2001 Plate 2 DATA #1961-1990 mean url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/reconstructions/n_hem_temp/briffa2001jgr2.txt" #readLines(url)[1:50] abd<-read.table(url,skip=22,fill=TRUE,sep="\t",header=TRUE) #1400-1995 temp=(abd< -90);abd[temp]=NA abd=ts(abd[,2:ncol(abd)],start=1400) f=function(x){f=filter.combine.pad(x,truncated.gauss.weights(11))[,2];f} plot(1400:1995,f(abd[,3])[(1400:1995)-1399],lwd=2,xaxs="i",xlim=c(1400,1995),type="l") abline(h=0,lty=2) ##COMPARE BRIFFA SELECTIONS: MAKE REGIONAL TEMP COLLATIONS source("http://www.climateaudit.org/scripts/treeconfig.functions.txt") #utility functions source("http://www.climateaudit.org/scripts/briffa/briffa.region.txt") #briffa.region is list of jones index cells for regions library(ncdf) url="d:/climate/data/gridcell/crutem3/crutem3.nc" #also did with hadcrut3 v<-open.ncdf(url) instr <- get.var.ncdf( v, v$var[[1]]) # 1850 2006 K=dim(instr);K# [1] 72 36 1885 #this is organized in 72 longitudes from -177.5 to 177.5 and 36 latitudes from -87.5 to 87.5 stat=array(NA,dim=c(9,3)) row.names(stat)=names(briffa.region) dimnames(stat)[[2]]=c("all","80_60","reported") stat[,3]=c(.76,.65,.77,.52,.53,.43,.72,.65,.69) region.temperature=NULL for (n in 1:9){ x=briffa.region[[n]];N=length(x) x=array(jonesinv(x),dim=c(N,2)) j<- 19+ floor(x[,1]/5 +.01); i<- 37+ floor(x[,2]/5 +.01); gridcell3<-array(NA,dim=c(K[3],N)) for (k in 1:N) { gridcell3[,k]=instr[i[k],j[k],]} index<-nrow(gridcell3)%%12 gridcell3<-rbind(gridcell3,array(NA,dim=c(12-index,N) ) ) #1884 56 annual<- apply(gridcell3,1,mean,na.rm=TRUE) #1884 annual= t (array(annual,dim=c(12,length(annual)/12) )) annual= ts( apply(annual[,4:9],1,mean,na.rm=TRUE), start=1850) combine=ts.union(abd[,n],annual) dimnames(combine)[[2]][2]=names(briffa.region)[[n]] stat[n,1:2]= c(cor(combine,use=use0)[1,2] , cor(combine[(1880:1960)-1399,],use=use0)[1,2] ) region.temperature=ts.union(region.temperature,annual) } stat #APR-SEP all 80_60 reported SEUR -0.01 0.09 0.76 NEUR -0.04 0.03 0.65 NSIB 0.56 0.74 0.77 ESIB 0.43 0.60 0.52 CAS 0.01 0.26 0.53 TIBP 0.24 0.29 0.43 ECCA 0.18 0.22 0.72 WNA 0.06 0.22 0.65 NWNA 0.02 0.07 0.69 #APR-AUD # all 80_60 reported #SEUR 0.01 0.11 0.76 #NEUR -0.06 0.02 0.65 #NSIB 0.57 0.76 0.77 #ESIB 0.39 0.58 0.52 #CAS 0.01 0.24 0.53 #TIBP 0.28 0.32 0.43 #ECCA 0.18 0.23 0.72 #WNA 0.08 0.25 0.65 #NWNA 0.04 0.09 0.69 par(mar=c(5,3,1,1)) barplot(t(stat),beside=TRUE,las=2) ##figure ## CHECK REGION LOCATIONS library(fields) n=1 x=briffa.region[[n]];N=length(x) x=array(jonesinv(x),dim=c(N,2)) plot(x[,2],x[,1],xlim=c(-177.5,177.5),ylim=c(2.5,87.5) ,pch=19,cex=.8) world(add=TRUE)