#REPLICATION CPS #Comments: #workhorse function is manniancps which yields NH, SH composites #manniancps =function(k,criterion=passing$whole,outerlist=outerlist.mbh,smoothmethod="mann", lat.adjustment= -1,info=instrumental.info) { #Stupid Pet Trick patches # 1. outerjones function creates two lists of proxies attached to a gridcell #one is outerlist.mbh and one is outerlist.sensible # 2. patch is available to move latitudes off center in accordance with mannian ##FUNCTIONS source("http://www.climateaudit.org/scripts/utilities.txt") source("http://www.climateaudit.org/scripts/mann.2008/utilities.txt") 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 source("http://www.climateaudit.org/scripts/mann.2008/collation.cps.txt") #0r load("d:/climate/data/mann08/Basics.tab") #Basics=list(period_MBH=period_MBH,cal=cal,ver=ver,criteria=criteria,passing=passing,outerlist=list(outerlist.sensible,outerlist.mbh),proxy,proxys,proxyscaled, # inst=inst,instrumental.info=instrumental.info,instr=list(target=instr.target,smooth=instr.smooth,smooth05=instr.smooth.05) ) use0="pairwise.complete.obs" compare=function(A,B) { #utility function used below compare=array(NA,dim=c(2,ncol(A))) dimnames(compare)[[2]]=dimnames(A)[[2]] for( i in 1:ncol(A)) { test=cbind(A[,i],B[,i]) compare[1,i]=cor(test,use=use0)[1,2] compare[2,i]=max(abs(test[,1]-test[,2]),na.rm=T) } compare } ##EMULATE AD1000 NETWORK WITH MANNIAN OPTIONS #k=9;criterion=passing$whole; outerlist=outerlist.mbh;lat.adjustment= -1; smoothmethod="mann";verbose="verbose" raw=manniancps(k=9,criterion=passing$whole, outerlist=outerlist.mbh,lat.adjustment= -1, smoothmethod="mann",verbose="verbose") names(raw) # "idsort" "Data" "grids" "regts" "rescaled" "recon" "working" ######## ##COMPARE TO MATLAB INTERMEDIATES ####### ##LOAD COLLATION OF INTERIM RESULTS FOR iHAD VERSION FROM UC MATLAB RUN download.file("http://www.climateaudit.org/data/mann.2008/uc/regts9.tab","temp.dat",mode="wb");load("temp.dat") download.file("http://www.climateaudit.org/data/mann.2008/uc/saveregts9.tab","temp.dat",mode="wb");load("temp.dat") download.file("http://www.climateaudit.org/data/mann.2008/uc/grids.mann.tab","temp.dat",mode="wb");load("temp.dat") download.file("http://www.climateaudit.org/data/mann.2008/idsort.mann.tab","temp.dat",mode="wb");load("temp.dat") glat=read.table("http://www.climateaudit.org/data/mann.2008/uc/glat_glon.txt") #2592 #this was collated in collation.recon.txt and is Mann's data dim(grids.mann) #996 21 dim(regts); #1996 1824 #this is the number of regrid cells in info (instrumental) #reshape regts dimnames(regts)[[2]]=1:ncol(regts) count=apply(!is.na(regts),2,sum); sum(!(count==0)) #15 #only count1=apply(!is.na(regts),1,sum); sum(count1==0) #1000 regts.mann=ts(regts[!(count1==0),!(count==0)],start=1000); rm(regts) dim(regts.mann) #996 15 dimnames(saveregts)[[2]]=1:ncol(saveregts) count=apply(!is.na(saveregts),2,sum); sum(!(count==0)) #15 #only count1=apply(!is.na(saveregts),1,sum); sum(count1==0) #1000 saveregts.mann=ts(saveregts[!(count1==0),!(count==0)],start=1000);rm(saveregts) dim(saveregts.mann) #996 15 rbind(idsort.mann,raw$idsort) #these match #raw$regts matches grids.mann (UC Matlab) compare(raw$regts,grids.mann) # 670 714 719 978 1102 1486 1487 1603 #[1,] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.0000e+00 #2,] 5.217661e-08 4.761989e-08 4.944521e-08 4.840176e-08 5.190548e-08 9.342388e-09 4.937191e-08 5.2275e-08 #raw$rescaled is re-grid version: virtually exact match to saveregts.mann t(apply(compare(raw$rescaled,saveregts.mann),1,range)) # [,1] [,2] #[1,] 9.999955e-01 1.000000000 #[2,] 2.978723e-10 0.003153969 ### ##COMPARE SH iHAD RECONSTRUCTION TO ARCHIVED #### #this version is chosen because it is a long unspliced version download.file("http://www.climateaudit.org/data/mann.2008/recon.cps.tab","temp.dat",mode="wb");load("temp.dat") #load("d:/climate/data/mann08/recon.cps.tab") #prior collation from Mann data archived=recon #19 archived series collated from Mann (These are spliced) names(archived) #[1] "NH_all_CRU" "NH_nodendro_CRU" "NH_noluter_CRU" "NH_all_iCRU" "NH_nodendro_iCRU" #[6] "NH_noluter_iCRU" "NH_all_HAD" "NH_nodendro_HAD" "NH_noluter_HAD" "NH_all_iHAD" #[11] "NH_nodendro_iHAD" "NH_noluter_iHAD" "SH_all_CRU" "SH_nodendro_CRU" "SH_noluter_CRU" #[16] "SH_all_iCRU" "SH_nodendro_iCRU" "SH_noluter_iCRU" "SH_all_iHAD" "NH_CRU_composite" #[21] "NH_HAD_composite" "SH_CRU_composite" "SH_HAD_composite" #COMPARE SH iHAD - a long unspliced version test=ts.union(raw$recon[,2], archived[["SH_all_iHAD"]][,1]); cor (test,use=use0) # 0.9996589 plot.ts(test) #scales are different #RESCALE TO SH iHAD smoothed cps = rescale(raw$recon[,"SH"],instr.smooth[,8],index=1850:1995) cor(cps,archived[["SH_all_iHAD"]][,1] ) range(cps- archived[["SH_all_iHAD"]][,1]) # -0.09686743 0.03110847 #range(instr.smooth-Basics$instr$smooth) #0 #cps = rescale(raw$recon[,"SH"],Basics$instr$smooth05[,8],index=1850:1995) #cor(cps,archived[["SH_all_iHAD"]][,1] ) #range(cps- archived[["SH_all_iHAD"]][,1]) # -0.09567279 0.03127168 layout(array(1:2,dim=c(2,1)),heights=c(1.1,1.3)) par(mar=c(0,3,2,1)) plot(c(time(cps)),cps,xlab="",ylab="",axes=FALSE,type="l",col=2) points(1000:1995,archived[["SH_all_iHAD"]][,1],pch=19,cex=.05) axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0,lty=2) legend("bottomright",fill=1:2,legend=c("Emulation","Archived"),cex=.8) title("Emulation of SH iHAD AD1000 CPS") par(mar=c(3,3,0,1)) plot(c(time(cps)),cps- archived[["SH_all_iHAD"]][,1],xlab="",ylab="",axes=FALSE,type="l",ylim=c(-.25,.25)) axis(side=1);axis(side=2,las=1);box();abline(h=0,lty=2)