# source("http://www.climateaudit.org/scripts/mann.2008/utilities.txt") codetablef=function() { #requires details ;transcribes gridproxy.m #this table collates Mann classification of proxies for screening operations codetable=data.frame(as.numeric(as.character(sort(unique(details$code))))) names(codetable)="code" codetable$res="A" codetable$res[codetable$code%%1000==1]="D" codetable$class="O" codetable$class[6:9]="B" #this is the documentary and speleothem treated separately codetable$cortype=factor(paste(codetable$res,codetable$class,sep="")) codetable$whole=codetable$cortype levels(codetable$whole) #[1] "AB" "AO" "DB" "DO" levels(codetable$whole)=c(.136,.106,.425,.337) codetable$EA100=codetable$cortype;levels(codetable$EA100)=c(.164,.128,.521,.418) codetable$full=codetable$cortype;levels(codetable$full)=rep(0,4) for (i in 5:7) codetable[,i]=as.numeric(as.character(codetable[,i])) #codetable transcribes values in gridproxy.m codetable$LA100=codetable$EA100 codetablef=codetable codetablef} (codetable=codetablef()) # code res class cortype whole EA100 full # 1 2000 A O AO 0.106 0.128 0 # 2 3000 A O AO 0.106 0.128 0 # 3 3001 D O DO 0.337 0.418 0 # 4 4000 A O AO 0.106 0.128 0 # 5 4001 D O DO 0.337 0.418 0 # 6 5000 A B AB 0.136 0.164 0 # 7 5001 D B DB 0.425 0.521 0 # 8 6000 A B AB 0.136 0.164 0 # 9 6001 D B DB 0.425 0.521 0 # 10 7000 A O AO 0.106 0.128 0 # 11 7001 D O DO 0.337 0.418 0 # 12 7500 A O AO 0.106 0.128 0 # 13 8000 A O AO 0.106 0.128 0 # 14 8001 D O DO 0.337 0.418 0 # 15 9000 A O AO 0.106 0.128 0 #1. General Information method="online" period_MBH=seq(1800,0,-100);M=length(period_MBH);M #the 19 relevant periods cal=list();cal$earlym=c(1896,1995);cal$latem=c(1850,1949);cal$whole=c(1850,1995) ver=list();ver$earlym=c(1850,1895);ver$latem=c(1950,1995); #these are the calibration periods #2. Make details/information of proxy network download.file("http://www.climateaudit.org/data/mann.2008/details.tab","temp.dat",mode="wb");load("temp.dat") coral=(details$code==7000)|(details$code==7001);sum(coral) #15 tempb= !is.na(match(details$code,c(5000,5001,6000,6001) ) );sum(tempb) #33 #doc|speleo with flips dendro=(details$code==9000)|(details$code==7500);sum(dendro) #1032 luter=(details$code==2000) annual=! (substr(as.character(details$code),4,4)=="1");sum(annual) #1158 NH=(details$lat>0) ;sum(NH) #1036 criteria=list(coral=coral,dendro=dendro,luter=luter,annual=annual,NH=NH,codetable=codetable) #Make information on qualifying "whole", early-miss and late-miss series from Mann SI (stipulating results) #these are 3 logical vectors of length 1209 #the vector distinguishes between annual and "decadal" passingf=function (x,lowf,method="whole",Details=details) { #this is for a list of 1209 proxies target=x decadal=(substr(details$code,4,4)=="1");#sum(decadal) #51 target[decadal]=lowf[decadal] #this selects A and D actuals coral=(details$code==7000)|(Details$code==7001) # #invert coral (7000,7001) need to be inverted target[coral]= - target[coral] #inverts coral sign for test #invert documentary and speleothem with negative correlations tempb= !is.na(match(Details$code,c(5000,5001,6000,6001) ) );sum(tempb) #33 temp1= (x<0) target[tempb&temp1]= - target[tempb&temp1] benchmark= codetable[match(Details$code,codetable$code),method] passingf= (target>=benchmark) passingf } passing=list() passing$whole= passingf(x=details$rtable.r1850_1995,lowf=details$rtable.r1850_1995lf,method="whole") #638 passing$earlym=passingf(x=details$rtable.r1896_1995,lowf=details$rtable.r1896_1995lf,method="EA100") #638 passing$latem=passingf(x=details$rtable.r1850_1949,lowf=details$rtable.r1850_1949lf,method="EA100") #638 #sapply(passing,length) # whole earlym latem # 1209 1209 1209 #make two lists of proxies attached to gridcell: one with Mann fencepost problem and one in a sensible way outerjones=function(x) unique( c(outer (x[1]+1+c(0,0),x[2]+c(0,0),jones) )) outerlist.sensible=apply(details[,c("lat","long")],1,outerjones);names(outerlist.sensible)=1:1209 #this removes the fencepost problem outerjones=function(x) unique( c(outer (x[1]+1+c(-.0001,.0001),x[2]+c(-.0001,.0001),jones) )) outerlist.mbh=apply(details[,c("lat","long")],1,outerjones);names(outerlist.mbh)=1:1209 #this does odd Mann inclusion #3. proxy collection used in various subsets if(method=="online"){ download.file("http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/data/itrdbmatrix","temp.dat") #download.file("http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/data/itrdbmatrix","d:/climate/data/mann.2008/itrdbmatrix") proxy=scan("temp.dat",n=-1)} else { proxy=scan("d:/climate/data/mann.2008/itrdbmatrix",n=-1)} proxy=t( array(proxy,dim=c(1210,length(proxy)/1210) )) ;dim(proxy) #2005 1210 proxy=proxy[4:nrow(proxy),] # from - 5 to 1996 proxy=ts(proxy[,2:ncol(proxy)],start=proxy[1,1]) #dim(proxy) #[1] 2002 1209 tsp(proxy) #- 5 1996 dimnames(proxy)[[2]]=1:ncol(proxy) #make sign adjustments for coral and applicable documentary/speleothem series as in gridproxy.m proxy[,coral]= - proxy[,coral] #invert sign for 15 coral series temp=!is.na(match(details$code,c(5000,6000)))& (details$rtable.r1850_1995<0) ;sum(temp) # 5 doc and speleo proxy[,temp]= -proxy[,temp] temp=!is.na(match(details$code,c(5001,6001))) & (details$rtable.r1850_1995lf<0);sum(temp) #9 doc and speleo proxy[,temp]= -proxy[,temp] #37 MB - you may want to save this #make once-for-all smooth and scale library(signal) proxys=proxy; ipts=10 #ipts set as 10 in Mann lowpass cutfreq=.05; bf=butter(ipts,2*cutfreq,"low");npad=1/(2*cutfreq);npad proxys[,!annual]=apply(proxy[,!annual],2,function(x) mannsmooth(x,M=npad,bwf=bf )) #decadal data cutfreq=.1;ipts=10 #ipts set as 10 in Mann lowpass bf=butter(ipts,2*cutfreq,"low");npad=1/(2*cutfreq);npad proxys[,annual]=apply(proxy[,annual],2,function(x) mannsmooth(x,M=5,bwf=bf ) ) #save(proxys,file="d:/climate/data/mann08/proxys.tab") #this takes a couple of minutes for me temp=(time(proxys)>=1850)&(time(proxys)<=1995) ;sum(temp) #his is period in gridboxcps mean.proxy=apply(proxys[temp,],2,mean,na.rm=T);names(mean.proxy)=dimnames(proxys)[[2]]; mean.proxy["385"] # 7.585251 sd.proxy=apply(proxys[temp,],2,sd,na.rm=T);names(sd.proxy)=dimnames(proxys)[[2]]; sd.proxy["385"] #0.2512297 proxyscaled = ts( scale(proxys,center=mean.proxy,scale=sd.proxy),start=tsp(proxys)[1]) #4. GRidded Instrumental instrumental.infof=function(instfilled=inst) { #requires inst info=1:2592;info=data.frame(info);names(info)="mannid" enough=scan("http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/data/instrument/enough") #1732 toofew=scan("http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/data/instrument/toofew") #860 info$mann=NA; test=match(1:2592, enough);temp=!is.na(test);info$mann[temp]=1 info$lat= c(t( array (rep( seq(-87.5,87.5,5), 72),dim=c(36,72) ) )) info$long= rep( seq(-177.5,177.5,5),36 ) info$jones=jones(lat=info$lat,long=info$long) info$mu=cos(info$lat*pi/180) regrid=data.frame(1:288) lat=seq(82.5,-82.5,-15) long=seq(-172.5,172.5,15) regrid$lat =gl(12,24,labels=lat) regrid$lat=as.numeric(as.character(regrid$lat)) regrid$long= rep(long,12) regrid$mu=cos(regrid$lat*pi/180) fred=t(array( rep(1:24,9),dim=c(24,9)));fred fred=c(t(array(fred,dim=c(3,72)))) fred1=1728+seq(0,72,24) info$regrid=NA info$regrid[info$lat<30]=info$mannid[info$lat<30] info$regrid[info$lat>30]=c(outer(fred,fred1,function(x,y) x+y)) info$regrid.lat=info$lat info$regrid.long=info$long temp=(info$lat>30) info$regrid.lat[temp]= 15* trunc(info$lat[temp]/15)+7.5 info$regrid.long[temp]= 15* (trunc(180+info$long[temp]/15)-180)+7.5 info$mean=apply(instfilled[1:146,],2,mean) info$sd=apply(instfilled[1:146,],2,sd) info$mean_earlym=apply(instfilled[47:146,],2,mean) info$sd_earlym=apply(instfilled[47:146,],2,sd) info$mean_latem=apply(instfilled[1:100,],2,mean) info$sd_latem=apply(instfilled[1:100,],2,sd) instrumental.infof=info instrumental.infof } require("R.matlab");require("Rcompression") download.file("http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/data/instrument/infilled2poles_1850.mat","temp.dat",mode="wb") inst=readMat("temp.dat") inst=inst$anngridinst; #dim(inst) #157 2592} else { load("d:/climate/data/mann08/uc/smoothinst.tab"); inst=smoothinst;rm(smoothinst) instrumental.info=instrumental.infof() #5. Hemisphere Targets g=function(X) ts(X[,2],start=X[1,1]) url="http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/data/instrument" version=c("CRU","iCRU","HAD","iHAD"); hemi=c("NH","SH") version=c(outer(version,hemi, function(x,y) paste(x,y,sep="_"))) instr.target=ts(array(NA,dim=c(157,8)),start=1850) for(i in 1:8) instr.target[,i] = g(read.table(file.path (url, paste(version[i],"reform",sep="_")))) dimnames(instr.target)[[2]]=version #apply(instr.target,2,mean) # CRU_NH iCRU_NH HAD_NH iHAD_NH CRU_SH iCRU_SH HAD_SH iHAD_SH #-0.1258726 -0.1086633 -0.1344140 -0.1078047 -0.1623185 -0.1547042 -0.2276879 -0.1878140 #smooth with 0.1 #http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/code/codeprepdata/doindexxx.m cutfreq=.1;ipts=10 #ipts set as 10 in Mann lowpass bf=butter(ipts,2*cutfreq,"low");npad=1/(2*cutfreq);npad instr.target=ts(instr.target[(1850:1995)- 1849,],start=1850) instr.smooth=ts(apply(instr.target,2,function(x) mannsmooth(x,M=npad,bwf=bf ) ),start=1850) #smooth with 0.05 #http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/code/codeprepdata/doindexxx.m cutfreq=.05;ipts=10 #ipts set as 10 in Mann lowpass bf=butter(ipts,2*cutfreq,"low");npad=1/(2*cutfreq);npad instr.target=ts(instr.target[(1850:1995)- 1849,],start=1850) instr.smooth.05=ts(apply(instr.target,2,function(x) mannsmooth(x,M=npad,bwf=bf ) ),start=1850) outerlist=list(outerlist.sensible,outerlist.mbh);names(outerlist)=c("sensible","mbh") Basics=list(period_MBH=period_MBH,cal=cal,ver=ver,criteria=criteria,passing=passing,outerlist=outerlist,proxy=proxy,proxys=proxys,proxyscaled=proxyscaled, inst=inst,instrumental.info=instrumental.info,instr=list(target=instr.target,smooth=instr.smooth,smooth05=instr.smooth.05) ) save(Basics,file="d:/climate/data/mann08/Basics.tab")