## REGIONAL CHRONOLOGY ANALYSIS ## #vaganov data from http://www.cru.uea.ac.uk/cru/people/briffa/yamal2009/data/vaganovetal1996.zip #PDP,PLB,PLL,PLO,PLP,PLR,POM,POZ,PPP #metadata transcribed from FOI source("http://www.climateaudit.info/scripts/utilities.txt") source("http://www.climateaudit.info/scripts/tree/utilities.treering.txt") infob=read.csv("http://www.climateaudit.info/data/tree/rwl/briffa2006/info_yamal_2006.csv") K=nrow(infob);K #chronologies for Briffa 2000, 2008 loc= "http://www.climateaudit.info/data/tree/crn/cru/yamal_2000.crn.tab" download.file(loc, "temp",mode="wb"); load("temp") yamal00=chron.crn loc= "http://www.climateaudit.info/data/tree/crn/cru/yamal_2008.crn.tab" download.file(loc, "temp",mode="wb"); load("temp") yamal08=chron.crn tsp(yamal08) #-202 1996 #measurement data for FOI list Tree=rep(list(NA),K); names(Tree)=infob$id for (i in 1:K) { loc= paste("http://www.climateaudit.info/data/tree/rwl/briffa2006/",infob$id[i],".rwl.tab",sep="") download.file(loc, "temp",mode="wb") load("temp") tree$data=infob$id[i] Tree[[i]]=tree } sapply(Tree, function(A) range(A$year) ) # russ119w russ061w russ064w russ002 russ123w russ180w russ001 pou_la polurula pll yamalAD russ035w plr #[1,] 1603 1674 1630 1696 1588 1828 1570 914 778 1557 -202 1782 1609 #[2,] 1990 1991 1990 1969 1990 1994 1968 1990 1892 1992 1996 1990 1992 # pdp plo russ101w russ084w #[1,] 1505 1684 1740 1767 #[2,] 1992 1992 1990 1990 #make regional dataset work=NULL for(i in 1:K) work=rbind(work, Tree[[i]]) nrow(work) # 144706 count.rcs=countf(work) # barplot(tapply(work$rw,work$data,mean,na.rm=T) ) yamalcount=countf(Tree[["yamalAD"]]) #Figure 1: from IPCC #Figure 2: Compare Yamal and Regional Core Counts #png("d:/climate/images/2012/yamal/counts.png",h=360,w=600) par(mar=c(3,3,2,1)) count=count.rcs;year=c(time(count)); N=length(count) ts.plot(count,col=2,ylim=c(0, 400)); polygon(xy.coords(x=c(year,rev(year) ), y=c(count,rep(0,N) ) ),col="salmon",border="red") polygon(xy.coords(x=c(time(yamalcount),rev(time(yamalcount) )) , y=c(yamalcount,rep(0,length(yamalcount))) ),col="cyan",border="cyan") legend("topleft",fill=c("salmon","cyan"),legend=c("Region","Yamal") ) title("Core Counts") #dev.off() #calculate RCS chronology, adjusting site mean to overall mean after 1600 temp=work$year>=1600 m2=m0=tapply(work$rw[temp],work$data[temp],mean,na.rm=T) ; m0 m1= mean(work$rw[temp],na.rm=T); m1 #49.9 #53.9 overall m0= m0/m1 work$mean=factor(work$data); levels(work$mean)= m0 work$mean=as.numeric(as.character(work$mean)) work$adj=work$rw/work$mean tapply(work$adj,work$data,mean,na.rm=T) #NOT ALL EQUALL #all equal Chron2=RCS.chronology(work[,c(1,2,3,7)],"nls") # 223746128 : 30.68207490 64.73985670 0.01170579 chron2=Chron2$series chron2=window(chron2, end=max(time(chron2)[count>=10]) ) #Figure 3: plot Chronology #png("d:/climate/images/2012/yamal/regional_chron.png",h=400,w=600) par(mar=c(3,3,2,1)) ts.plot(window(chron2,800),type="l",col=1,lwd=1) title("Yamal-Urals Regional Chronology") #dev.off() #Figure 4a: compare regional chronology to Yamal #scale on 800-1900 (for contrast) #png("d:/climate/images/2012/yamal/chron_siteadj2.png",h=400,w=600) f=function(x) filter.combine.pad(x,truncated.gauss.weights(31))[,2] par(mar=c(3,4,2,1)) x= window(yamal08,800); u= mean(window(x,end=1900)); sd0= sd(window(x,end=1900)) year=c(time(x)) plot(year,f(scale(x,center=u,scale=sd0)),type="l",col=1,lwd=2,ylab="SD Units") v= mean(window(chron2,end=1900)); sd1= sd(window(chron2,end=1900)) points(c(time(chron2)),f(scale(chron2,center=v,scale=sd1)),col=3,lwd=2,lty=3,cex=.2) legend("topleft",fill=c(1,3),legend=c("Briffa 2008", "Regional- Site Adj")) title("Yamal and Regional Chronologies") #dev.off() #Figure 4b: scatter plot X=ts.union(yamal=window(yamal08,800,1994),region=window(chron2,800,1994)) X=data.frame(X);X$year=800:1994 #png("d:/climate/images/2012/yamal/scatter.png", h=480,w=480) par(mar=c(4,4,2,1)) plot(yamal~region,X,xlab="Regional",ylab="Yamal") points(yamal~region,X[(1850:1925)-799,],col="yellow",pch=19) points(yamal~region,X[(1926:1975)-799,],col="orange",pch=19) points(yamal~region,X[(1976:1994)-799,],col="red",pch=19) legend("bottomright",fill=c("yellow","orange","red"), legend=c("1850-1925","1926-1975","1976-1994"),cex=1.5 ) title("Yamal vs Regional Chronology") #dev.off() #Compare to counts in Briffa and Melvin 2009 (CRU website) #collate CRU datasets for Yamal id=c("SF","JAH","KHAD","POR","YAD") combo=NULL for (i in 1:5) { loc=paste("http://www.climateaudit.info/data/tree/rwl/briffa2009/",id[i],".rwl.tab",sep="") download.file(loc,"temp",mode="wb") load("temp") tree$data=id[i] combo=rbind(combo,tree) } nrow(combo) # [1] 49787 barplot(tapply(combo$rw,combo$data,mean,na.rm=T) ) combocount=countf(combo) #Plot Core Counts #png("d:/climate/images/2012/yamal/counts_combo.png",h=360,w=600) par(mar=c(3,3,2,1)) count=window(count.rcs,800);year=c(time(count));N=length(year) plot(year,count,type="l",col=2,ylim=c(0, 400),xlim=c(800,2000),xaxs="i"); polygon(xy.coords(x=c(year,rev(year) ), y=c(count,rep(0,N) ) ),col="salmon",border="red") polygon(xy.coords(x=c(time(combocount),rev(time(combocount) )) , y=c(combocount,rep(0,length(combocount))) ),col="violet",border="cyan") polygon(xy.coords(x=c(time(yamalcount),rev(time(yamalcount) )) , y=c(yamalcount,rep(0,length(yamalcount))) ),col="cyan",border="cyan") legend("topleft",fill=c("salmon","violet","cyan"), legend=c("Regional","CRU 2009","Yamal") ) title("Core Counts") #dev.off()