##MBH DATA INPUT mbhdata="http://www.climateaudit.org/data/mbh98" pcstore<-"http://www.climateaudit.org/data/mbh98/pcs" # Directory were PC will be stored during this collation ##TEMPERATURE PRINCIPAL COMPONENT INFORMATION USED IN PROXY CALIBRATION AND RECONSTRUCTION # temperature principal components download.file( file.path(mbhdata,"pc.tab"),"temp.dat",mode="wb") load("temp.dat") pc.MBH<-pc # temperature eofs used in reconstruction #these are ordered as Jones gridcell order (hour hand N to S; minute hand: -180 to 180) download.file( file.path(mbhdata,"eof.ordermann.tab"),"temp.dat",mode="wb") load("temp.dat") eof.MBH<-eof # temperature eigen-values and choose first 16 as only 16 PCs are used download.file( file.path(mbhdata,"tpca-eigenvals.tab"),"temp.dat",mode="wb") load("temp.dat") lambda<-lambda[1:16,2] lambda.MBH<-lambda # STEPWISE RECONSTRUCTION PARAMETERS #there are 11 steps and different temperature PCs used in each step period.MBH<-c(1981,1820,1800,1780,1760,1750,1730,1700,1600,1500,1450,1400) #sub-interval boundary period<-period.MBH-1399 #eigenvectors selected in each sub-interval- one less than "period" select<-list(c(1:5,7,9,11,14:16),c(1:5,7,9,11,14:16),c(1:5,7,9,11,14:16), c(1:5,7,9,11,15),c(1:3,5,7,9,11,15),c(1,2,5,11,15),c(1,2,5,11,15), c(1,2,11,15),c(1,2),c(1,2),c(1)) select.MBH<-select ##LOAD PROXY #this is a collation of the non-tree ring PC proxies #columns are preserved for 31 non-tree ring proxies from columns 68-99 #for each reconstruction step, the specified tree ring PCs will be pasted in #TO EXTEND SERIES BY PERSISTENCE extend.persist<-function(tree) { extend.persist<-tree for (j in 1:ncol(tree) ) { test<-is.na(tree[,j]) end1<-max ( c(1:nrow(tree)) [!test]) test2<-( c(1:nrow(tree))>end1) & test extend.persist[test2,j]<-tree[end1,j] } extend.persist } download.file( file.path(mbhdata,"proxy.tab"),"temp.dat",mode="wb") load("temp.dat") proxy.MBH<-ts( proxy[1:581,],start=1400,end=1980) proxy.MBH[,c(1:68,100:112)]<-extend.persist(proxy.MBH[,c(1:68,100:112)]) # weights used in proxy regression download.file( file.path(mbhdata,"weights1.tab"),"temp.dat",mode="wb") load("temp.dat") weightsprx=weights1 ##TREE RING PC INFORMATION # there are 6 different tree ring regions with varying number of PCs used and varying calls to subdirectories # Names of tree ring regions region<-c("STAHLE/OK","STAHLE/SWM","ITRDB/NOAMER","ITRDB/SOAMER","ITRDB/AUSTRAL","VAGANOV") regionid<-c("STAHLE.OK","STAHLE.SWM","NOAMER","SOAMER","AUSTRAL","VAGANOV") # for storage # Parameters for tree ring PC collection mbh98.backto<-list(c(1700),c(1400,1450,1500,1600,1700),c(1400,1450,1600,1750),c(1600),c(1600,1750),c(1450,1600,1750)) #Labels on PC files mannpc<-c(3,9,9,3,4,3) # Number of PC pc*.out files to read mannid<-c(69,72,84,93,96,81) # Location of first column of regional PC in proxy collation # Schedules for the number of PCs by calculation step and region; and call subdirectories by calculation step and region series.MBH<-read.table ( file.path(mbhdata,"pc_count.dat"),sep="\t",header=TRUE) #table of retained PCs by period and region directory.MBH<-read.table( file.path(mbhdata,"pc_directory.dat"),sep="\t",header=TRUE) #table of look-up directories for above row.names(series.MBH)=series.MBH[,1] row.names(directory.MBH)=directory.MBH[,1] series.MBH<-series.MBH[2:ncol(series.MBH)] directory.MBH<-directory.MBH[2:ncol(directory.MBH)] ##GRIDCELL INFORMATION # MBH uses 1082 gridcells in their "dense" reconstruction and 219 gridcells in their "sparse" reconstruction specified here # load list of MBH "dense" grid-cells download.file( file.path(mbhdata,"gridpoints.tab"),"temp.dat",mode="wb") load("temp.dat") tempmann<-sort(grid2[,5]) manngrid<-!is.na(match(c(1:2592),tempmann)) ordermann<-order(grid2[,5]) grid2<-grid2[ordermann,] nhtemp<-(grid2[,4]>0) #logical vector for NH gridcells mu<-cos(grid2[,4]*pi/180) #used for areal weighting mu.MBH<-mu names(grid2)=c("i","j","greenwich","lat","jones") grid2$long=grid2$greenwich grid2$long[grid2$long>180]= grid2$long[grid2$long>180]-360 # mask for sparse subset of 219 cells #mask nowhere specified but determiend by careful inspection of diagram download.file( file.path(mbhdata,"mask.edited.tab"),"temp.dat",mode="wb") load("temp.dat") mask.logical<-!is.na(match(grid2[,5],mask.edited)) #sum(mask.logical) #219 #length(mask.logical) #mask<-list(mask.edited,mask.logical) grid2$mask=mask.logical plot(grid2$long[grid2$mask],grid2$lat[grid2$mask],pch=19,xlim=c(-177,177),ylim=c(-70,80));#world(add=TRUE) # allow for gridcells with 0 values #this is no longer used after July 2004 Corrigendum information #load(file.path(mbhdata,"nil.tab"))#load list of 4 grid-cells identified as having no values. See read.jones for discussion. nil<-rep(FALSE,1082) #this uses new SI information (May 2004) nhtemp.nil<-nhtemp&!nil nhtemp.MBH<-nhtemp.nil # load list of grid-cell standard deviations # truncate to exclude nil cells (allowed for but not applicable here) download.file( file.path(mbhdata,"gridcell.stdev.ordermann.tab"),"temp.dat",mode="wb") load("temp.dat") mann.scale<-gridcell.stdev.ordermann #no nil values #nhscale<-1/(mann.scale[nhtemp.nil]*sum(nhtemp.nil)) # #nhscale<-array( nhscale, dim= c(length(nhscale),1) ) #dim ## TEMPERATURE AND RECONSTRUCTION DATA download.file( file.path(mbhdata,"nhmann.tab"),"temp.dat",mode="wb") load("temp.dat") #nhmann #NH.scale.factor<-sd(nhmann[503:581,3]) #NH.MBH<-ts(nhmann[,2],start=1400) # archived "sparse" data #no information at NOAA on zero but dense is zero on 1902-1980 by experiment #url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo/paleocean/by_contributor/mann1998/nhem-sparse.dat" #g=read.table(url,skip=1) #sparse<-ts(g[,2:3],start=g[1,1],end=g[nrow(g),1]) #dimnames(sparse)[[2]]<-c("instr","recon") #sparse[(1981:1993)-1853,2]=NA download.file( file.path(mbhdata,"sparse.tab"),"temp.dat",mode="wb") load("temp.dat") #sparse sparse=sparse[,1] sparse<-ts(scale(sparse,center=mean(sparse[(1902:1980)-1853]),scale=FALSE),start=1854) sparse=sparse[,1] ###PARTICULARS download.file( file.path(mbhdata,"prname.tab"),"temp.dat",mode="wb") load("temp.dat") #prname ###RPCWEIGHT nh=nhtemp sd.rescaling=rep(1,1082) rpcweights.mbh=diag(lambda[1:16])%*% t(eof[nh,1:16]) %*% diag(mann.scale[nh])%*% as.matrix(sd.rescaling[nh]) /sum(mu[nh]) c(rpcweights.mbh) # [1] 6.67309067 -1.45784593 -0.09858368 -2.02105147 1.37839440 sum(mask.logical) #219 rpcweights.mbh_sparse=diag(lambda[1:16])%*% t(eof[nh&mask.logical,1:16]) %*% diag(mann.scale[nh&mask.logical])%*% as.matrix(sd.rescaling[nh&mask.logical]) /sum(mu[nh&mask.logical]) c(rpcweights.mbh_sparse) # [1] 5.516241 -4.501015 2.455701 -2.602087 4.173419 -2.033396 .... rpcweights.mbh=cbind(rpcweights.mbh,rpcweights.mbh_sparse);dimnames(rpcweights.mbh)[[2]]=c("nh","sparse") mannomatic=function(b,M=78,sdmethod="detrended") { N<-nrow(b) m2<-apply(b[(N-M):N,],2,mean) s2<-apply(b[(N-M):N,],2,sd) if (sdmethod=="unscaled") d<-scale(b,center= m2,scale=FALSE) if (sdmethod=="undetrended") d<-scale(b,center= m2, scale= apply(b[(N-M):N,],2,sd) ) if (sdmethod=="detrended") { d<-scale(b,center= m2, scale= apply(b[(N-M):N,],2,sd) ); sdprox<-apply(d[(N-M):N,],2,sd.detrend) d<-scale(d,center= FALSE,scale= sdprox)} mannomatic<-d mannomatic } sd.detrend=function(x) { t<-c(1:length(x)) fm<-lm(x~t) sd.detrend<-sd(fm$residuals) sd.detrend } ##LOADS THE PROXY WITH PCS FOR SPECIFIED STEP #call: proxy<-loadproxy(step0=k,case,proxy.method,mbhdata,series,directory,case.NOAMER,gaspe) ##options are as follows: #case - lookup directory for non-NOAMER proxies #proxy.method: MBH98 - load proxies and PCs; WA - use WA methods; noPC - load individual proxies #mbhdata - is directory address for proxy data. Set at "c:/climate/data/mann" #series.MBH - is lookup table of retained PCs by region and calculation step. Benchmarked off NAture SI #directory.MBH - is lookup table by period and region of period used for PCs in calculation step. Benched off nature SI #case.NOMER - lookup directory for North American tree ring proxies #gaspe.method: default - uses MBH98; archived - NA for first 4 years; updated - uses new version case="MBH98";proxy.method="MBH98";mbhdata0=mbhdata;series=series.MBH;directory=directory.MBH;case.NOAMER="MBH98";gaspe.method="default" loadproxy<-function(step0,case="MBH98",proxy.method="MBH98",mbhdata0=mbhdata,series=series.MBH,directory=directory.MBH,case.NOAMER="MBH98",gaspe.method="default"){ if (proxy.method=="MBH98") {proxy=proxy.MBH proxy<-ts(proxy[1:581,],start=1400,end=1980) k<-step0 case.temp<-c(rep(case,2),case.NOAMER,rep(case,3)) for (m in 1:6) { if (series[m,k]>0) {loc=file.path(pcstore,paste(regionid[m],directory[m,k],"tab",sep=".")) download.file(loc,"temp.dat",mode="wb") load("temp.dat")} #this proides for using CENSORED version if (series[m,k]>0) { proxy[,(mannid[m]-1+(1:series[m,k]))]<-pca[,(1:series[m,k])] } } temp<-!is.na(proxy[period[k+1],]) proxy[,temp]<-extend.persist(proxy[,temp]) if (gaspe.method=="archived") proxy[1:4,53]<-NA if (gaspe.method=="updated") { proxy[,53]<-NA; load("c:/climate/data/erren/gaspe.tab"); proxy[178:581,53]<-gaspe[1:404] } for ( i in 1:ncol(proxy[,temp])){proxy[,temp][,i]<-normalize(proxy[,temp][,i])} loadproxy<-proxy } if (proxy.method=="WA") { data.objects<-format.datamatrices(param.input,12-step0) ;#basic WA input function loadproxy<-data.objects$Ymat.anomaly } if (proxy.method=="noPC") { load(file.path(mbhdata0,"proxy.tab")) proxy<-ts(proxy[1:581,],start=1400,end=1980) name0<-dimnames(proxy)[[2]] temp<-!is.na(proxy[period[step0+1],]) name0<-name0[temp] proxy<-proxy[,temp] for (m in 1:6) { tree<-read.table (file.path("http:/www.climateaudit.org/data/mbh98/UVA",paste(regionid[m],"txt",sep=".")),header=TRUE,sep="\t" ) M0<-max(1400,tree[1,1]) M1<-min(1980,tree[nrow(tree),1]) tree<-tree[(M0-tree[1,1]+1):(M1-tree[1,1]+1),] tree<-ts(tree[,2:ncol(tree)],start=tree[1,1],end=tree[nrow(tree),1]) proxy<-ts.union(proxy,tree)} proxy<-proxy[,!is.na(proxy[period[step0+1],])] proxy<-extend.persist(proxy) dimnames(proxy)[[2]][1:length(name0)]<-name0 if (gaspe.method=="archived") proxy[1:4,"V53"]<-NA if (gaspe.method=="updated") { proxy[,"V53"]<-NA; load("c:/climate/data/erren/gaspe.tab"); proxy[178:581,"proxy.V53"]<-gaspe[1:404] } for ( i in 1:ncol(proxy)){proxy[,i]<-normalize(proxy[,i])} loadproxy<-proxy } loadproxy } normalize=function(x,M1=1902,M2=1980) (x-mean(x[(M1:M2)-tsp(x)[1]+1]))/sd(x[(M1:M2)-tsp(x)[1]+1]) hs=function(x) {N=length(x); (mean(x[(N-78):N]) -mean(x,na.rm=T))/sd(x,na.rm=T)} ssq=function(x) sum(x^2,na.rm=T) use0="pairwise.complete.obs"