##COLLATE RWL DATA AS NEWLY ARCHIVED #1. MAKE R TABLES FROM RECENTLY ARCHIVED DATA source("d:/climate/scripts/utilities.treering.txt") id= paste("cana",223:242,sep="") N<-length(id) url1= "ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/updates" for (i in 1:length(id) ){ tree<-make.rwl(id[i],url0=url1) #MALKE RWL out<-file.path("d:/climate/data/tree/measurements/r",country.out(id[i]),paste(id[i],"rwl.tab",sep=".") );out if(class(tree)=="data.frame") save(tree,file=out) } id= paste("cana",c(225,226,230,231,234,237,238),"x",sep="") N<-length(id) for (i in 1:length(id) ){ tree<-make.rwl(id[i],url0=url1) #MALKE RWL out<-file.path("d:/climate/data/tree/measurements/r",country.out(id[i]),paste(id[i],"rwl.tab",sep=".") );out if(class(tree)=="data.frame") save(tree,file=out) } ##2. MAKE AND SAVE JACOBY-STYLE CHREONOLOGIES source("d:/climate/scripts/treeconfig.functions.txt") id= paste("cana",223:242,sep="") N<-length(id) for (i in 2:length(id) ){ inbox=file.path("d:/climate/data/tree/measurements/r/northamerica",paste(id[i],"rwl.tab",sep=".") ) load(inbox); chron=jacoby.chronology(tree) out<-file.path("d:/climate/data/tree/chronologies/r/northamerica",paste(id[i],"crn.tab",sep=".") );out chron.crn=chron$series save(chron.crn,file=out) } ##3. LOAD ARCHIVED WILSON 2007 #readme:ftp://ftp.agu.org/apend/jd/2006jd008318 url="ftp://ftp.agu.org/apend/jd/2006jd008318/2006jd008318-ts01.txt" wilson=read.table(url,header=TRUE,fill=TRUE)#255 16 wilson=ts(wilson[,2:16],start=wilson[1,1]) #BC SERIES is IBC ##TEMPERATURE RECONSTRUCTION The TR proxy series must correlate at >0.40 against an optimal seasonal parameter of ‘‘local’’ gridded mean temperature data from the CRU3 [Brohan et al., 2006] land only data set. Even if a series correlates at <0.40, but the inferred association is significant at the 95% confidence limit, it was still rejected from further analysis. source("d:/climate/scripts/gridcell/collation.functions.txt") library(ncdf) url= "d:/climate/data/gridcell/crutem3/crutem3.nc" v<-open.ncdf(url) instr <- get.var.ncdf( v, v$var[[1]]) # 1850 2006 #dim(instr)# [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 lat=52; long= -117.5 j<- 19+ floor(lat/5 +.01); i<- 37+ floor(long/5 +.01); gridcell3<-instr[i,j,] index<-length(gridcell3)%%12 gridcell3<-c(gridcell3,rep(NA,12-index) ) gridcell3<-ts(gridcell3,start=c(1850,1),freq=12) X=array(gridcell3,dim=c(12,length(gridcell3)/12)) test1=ts(apply(X[5:8,],2,mean,na.rm=TRUE),start=tsp(gridcell3)[1]) ##4. COLLATE JACOBY-STYLE CHRONOLOGIES #Ring-width series were detrended using negative exponential or negative linear regression functions # and averaged to form site standard chronologies. id= paste("cana",223:242,sep="") chron=NULL for (i in 1:length(id) ){ inbox=file.path("d:/climate/data/tree/chronologies/r/northamerica",paste(id[i],"crn.tab",sep=".") ) load(inbox); chron=ts.union(chron,chron.crn) } dimnames(chron)[[2]]=id tsp(chron) # 1477 1998 annual=ts(apply(chron,1,mean,na.rm=TRUE),start=tsp(chron)[1]) count=apply(!is.na(chron),1,sum); temp=(count>2) annual=ts(annual[temp],start=min(time(annual)[temp])) plot(time(annual)[temp],annual[temp],type="l") chron.rw=chron;annual.rw=annual id= paste("cana",c(225,226,230,231,234,237,238),"x",sep="") chron=NULL for (i in 1:length(id) ){ inbox=file.path("d:/climate/data/tree/chronologies/r/northamerica",paste(id[i],"crn.tab",sep=".") ) load(inbox); chron=ts.union(chron,chron.crn) } dimnames(chron)[[2]]=id tsp(chron) # 1512 1997 annual=ts(apply(chron,1,mean,na.rm=TRUE),start=tsp(chron)[1]) count=apply(!is.na(chron),1,sum); temp=(count>2) annual=ts(annual[temp],start=min(time(annual)[temp])) chron.mxd=chron;annual.mxd=annual combine=ts.union(annual.rw,annual.mxd,wilson[,"IBC"],test1) combine=data.frame(combine);names(combine)=c("rw","mxd","wilson","temp") fm=lm(wilson~rw+mxd,data=combine) cor(combine,"kendall",use=use0) # annual.rw annual.mxd wilson[, "IBC"] test1 #annual.rw 1.0000000 0.1557734 0.3424418 0.1340600 #annual.mxd 0.1557734 1.0000000 0.7125308 0.5205223 #wilson[, "IBC"] 0.3424418 0.7125308 1.0000000 0.7638708 #test1 0.1340600 0.5205223 0.7638708 1.0000000 ##5. PLOT ANNUAL AVERAGES annual=ts.union(annual,wilson[,"IBC"]) f=function(x){f=filter.combine.pad(x,truncated.gauss.weights(21))[,2];f} index=time(annual) nf=layout(array(1:3,dim=c(3,1)),heights=c(1.1,1,1.3)) par(mar=c(0,3,2,1)) plot(c(index),annual[,1],type="l",axes=FALSE,xlab="",ylab="RW",col="grey70") lines(c(index),f(annual[,1]),lwd=2) axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=1,lty=2) text(1560,1.33,"Average of 20 RW ",font=2,pos=4) title(main="Wilson Data Archived July 2007: Jacoby-type Chronologies") par(mar=c(0,3,0,1)) plot(c(index),annual[,2],type="l",axes=FALSE,xlab="",ylab="RW",col="grey70") lines(c(index),f(annual[,2]),lwd=2) axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=1,lty=2) text(1560,1.1,"Average of 7 MXD ",font=2,pos=4) par(mar=c(3,3,0,1)) plot(c(time(annual)),annual[,3],type="l",axes=FALSE,xlab="",ylab="RW",col="grey70") lines(c(time(annual)),f(annual[,3]),lwd=2) axis(side=1);axis(side=2,las=1);box();abline(h=1,lty=2) text(1560,2,"Wilson et al 2007 IBC",font=2,pos=4) ##6. RCS-type CHRONOLOGIES id= paste("cana",223:242,sep="") N<-length(id) chron=NULL;fit=NULL for (i in 1:length(id) ){ inbox=file.path("d:/climate/data/tree/measurements/r/northamerica",paste(id[i],"rwl.tab",sep=".") ) load(inbox); test=RCS.chronology(tree,method="RCS.spline") ##neg exp fits don't work chron=ts.union(chron,test$series) fit=ts.union(fit,ts(test$coefficients,start=1)) } annual=ts(apply(chron,1,mean,na.rm=TRUE),start=tsp(chron)[1]) count=apply(!is.na(chron),1,sum); temp=(count>2) annual=ts(annual[temp],start=min(time(annual)[temp])) chron.rw.RCS=chron;annual.rw.RCS=annual;fit.RCS=fit id= paste("cana",c(225,226,230,231,234,237,238),"x",sep="") N<-length(id) chron=NULL;fit=NULL for (i in 1:length(id) ){ inbox=file.path("d:/climate/data/tree/measurements/r/northamerica",paste(id[i],"rwl.tab",sep=".") ) load(inbox); test=RCS.chronology(tree,method="linear") chron=ts.union(chron,test$series) fit=ts.union(fit,ts(test$coefficients,start=1)) } annual=ts(apply(chron,1,mean,na.rm=TRUE),start=tsp(chron)[1]) count=apply(!is.na(chron),1,sum); temp=(count>2) annual=ts(annual[temp],start=min(time(annual)[temp])) chron.mxd.RCS=chron;annual.mxd.RCS=annual;fit.mxd=fit ##7. PLOT RCS ANNUAL AVERAGES annual=ts.union(annual.rw.RCS,annual.mxd.RCS,wilson[,"IBC"]) f=function(x){f=filter.combine.pad(x,truncated.gauss.weights(21))[,2];f} index=time(annual) nf=layout(array(1:3,dim=c(3,1)),heights=c(1.1,1,1.3)) par(mar=c(0,3,2,1)) plot(c(index),annual[,1],type="l",axes=FALSE,xlab="",ylab="RW",col="grey70") lines(c(index),f(annual[,1]),lwd=2) axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=1,lty=2) text(1560,1.33,"Average of 20 RW ",font=2,pos=4) title(main="Wilson Data Archived July 2007: RCS-type Chronologies") par(mar=c(0,3,0,1)) plot(c(index),annual[,2],type="l",axes=FALSE,xlab="",ylab="RW",col="grey70") lines(c(index),f(annual[,2]),lwd=2) axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=1,lty=2) text(1560,1.1,"Average of 7 MXD ",font=2,pos=4) par(mar=c(3,3,0,1)) plot(c(time(annual)),annual[,3],type="l",axes=FALSE,xlab="",ylab="RW",col="grey70") lines(c(time(annual)),f(annual[,3]),lwd=2) axis(side=1);axis(side=2,las=1);box();abline(h=1,lty=2) text(1560,2,"Wilson et al 2007 IBC",font=2,pos=4) ##SEMIPAR FITS library(SemiPar);i=222 load("d:/climate/data/tree/northamerica.details.tab") test=details[!is.na(match(details$id,paste("cana",223:242,sep=""))),1:8];name0=as.character(test$location) rm(test);name0 # [1] "Glacier" "Revelstoke" "Meadow Mountain" "Kokanee" "Idaho Lookout" "Ymir" "Old Glory" # [8] "Kootenay" "Park Mountain" "Silverstar" "Galloping Mountain" "Big White" "Mt. Baldy" "Grey Creek" #[15] "Pavilion Mountain" "Cornwall Hills" "Mt. Chuwells" "Promontory" "Zum Peak" "Lone Man" index= paste("cana",223:242,sep="") N<-length(index) nf=layout(array(1:20,dim=c(5,4)),heights=c(1.1,1,1,1,1.3)); par(mar=c(2,2,1,1)) for (i in 223:242){ load(file.path("d:/climate/data/tree/measurements/r/northamerica",paste(index[i],"rwl.tab",sep="."))) temp=!is.na(tree$rw) tree=tree[temp,] attach(tree) mean(tree$rw) if(mean(tree$rw)>500) detach(tree) tree$rw=tree$rw/10 attach(tree) fit <- spm(rw~f(age),random=~1,group=id) plot(fit,main=paste(name0[i-222]),ylim=c(40,170),xlab=""); detach(tree) } index= paste("cana",c(225,226,230,231,234,237,238),"x",sep="") test=details[!is.na(match(details$id,index)),1:8];name0=as.character(test$location) N<-length(index) nf=layout(array(1:20,dim=c(4,2)),heights=c(1.1,1,1,1.3)); par(mar=c(2,2,1,1)) for (i in 1:7){ load(file.path("d:/climate/data/tree/measurements/r/northamerica",paste(index[i],"rwl.tab",sep="."))) temp=!is.na(tree$rw) tree=tree[temp,] attach(tree) fit <- spm(rw~f(age),random=~1,group=id) plot(fit,main=paste(name0[i]),xlab=""); }