#SHOW BENCHMARK manniancps #used in http://climateaudit.org/2010/07/30/make-a-stick-make-a-stick/#comment-237164 #re-copied from utilities.v4.txt of Nov 9, 2010; #.org replaced with .info in JUly 2010 # see discussion at #http://www.climateaudit.org/?p=4244 #http://www.climateaudit.org/?p=4274 #http://www.climateaudit.org/?p=4494 #http://www.climateaudit.org/?p=4501 # includes emulation options for various 'stupid pet tricks' ############### # DATA WRAPPER ################### # Basics - wrapper list collecting various Mannian objects (see mann.s008/collation.cps.txt for collation) # the wrapper is unpacked to give the following items used in various places in manniancps # details - dataframe with information on 1209 Mann proxies # period_MBH - vector of length 19 (1800 back to 0) giving start of Mannian step # cal - list of length 3 with start and end of Mann calibration period (earlym, latem, whole) # ver - list of length 3 with start and end of Mann verification period (earlym, latem, whole). "whole" is same as calibration # criteria- list: "coral" "dendro" "luter" "annual" "NH" "codetable" "tempb" # coral - logical vector identifying 15 coral proxies # dendro - logical vector identifying 1032 dendro proxies. These include a few different codes # luter- logical vector identifying 71 Luterbacher instrumental non-proxies. # annual- logical vector identifying 1158 'annual' proxies # NH- logical vector identifying 1036 NH proxies # tempb=logical vector identifying 33 "two-way" proxies - codes 5000,5001,6000,6001 # passing - list of 3 logical vectors (latem, earlym, whole) reporting Mannian 'pass' results # use0="pairwise.complete.obs" #utility for cor function # outerlist.mbh - list of length 1209 attaching Mannian gridcells to proxies. Permits duplicates as in Mann stupid pet trick # outerlist.sensible -vector of length 1209 attaching each proxy to one gridcell # proxy - dimension 2002 1209: original data collated in collation.originalproxy.txt # proxys - dimension 2002 1209: after Mannian smoothing. Avoids re-doing Mannian smoothing thousands of times. # proxyscaled- dimension 2002 1209: after Mannian smoothing and scaling. Avoids re-doing Mannian smoothing thousands of times. # inst - dim 157 2592. Insrumental data as annual time series 1850-2006 for 72x36 . # instrumental.info - 2592x15. Information on gridcells including concordance of different gridding systems # mannid mann lat long jones mu regrid regrid.lat regrid.long mean_whole sd_whole mean_earlym sd_earlym mean_latem #1 1 NA -87.5 -177.5 2521 0.044 1 -87.5 -177.5 -0.196 0.240 -0.172 0.269 -0.296 #2 2 NA -87.5 -172.5 2522 0.044 2 -87.5 -172.5 -0.196 0.240 -0.172 0.269 -0.296 #see instrumental.info.txt # instr.target - 157 8: tsp 1850-2006 NH and SH by CRU, iCRU,HAD, iHAD # instr.smooth - 157 8 tsp 1850-2006. smoothed version of above # alias - names for late- early and whole periods # outerlist.mbh - this is a list of associated gridcells for 1209 proxies. More than one gridcell possible per proxy #this is collated once and saved # outerlist.sensible - this is a list of associated gridcells for 1209 proxies - one gridcell per proxy #this is collated once and saved # smoothmethod (tag): # mbh - Mannian. Done separately for each step, extrapolating values in mannsmooth style. # "sensible" - # lat.adjustment (tag): default is -1. # screenmethod (tag) : Mann #base case returns NH and SH reconstructions, which are rescaled against instrumental targets in function cpsf # verbose: default - results in simplied return; "verbose" - gives intermediates for diagnostics # case: three options- whole, earlym, latem # outerlist: two options- outerlist.mbh; outerlist.sensible # criterion: this permits calculations using dendro, noluter, etc. # smoothmethod: two options- mann (piecemeal smoothing); sensible (smooth entire series) # lat.adjustment: default is quirky Mann displacement in which gridcells are centered 1 deg S of "typical" 2.5 deg centers; # M08 input data Basics is required; also passing; cal; ver; period_MBH # k=9;case="latem";criterion=rep(TRUE,1209);outerlist=outerlist.mbh;smoothmethod="mann";lat.adjustment= -1;verbose="default";Info=Basics # outerlist=outerlist.mbh=Basics$outerlist$mbh;smoothmethod="mann";lat.adjustment= -1;verbose="verbose" # period_MBH=Basics$period_MBH;cal=Basics$cal;ver=Basics$ver; criteria=Basics$criteria;passing=Basics$passing; # coral=criteria$coral;dendro=criteria$dendro;luter=criteria$luter;annual=criteria$annual;NH=criteria$NH;tempb=criteria$tempb # use0="pairwise.complete.obs" # outerlist.mbh=Basics$outerlist$mbh;outerlist.sensible=Basics$outerlist$sensible; # proxy=Basics$proxy;proxys=Basics$proxys;proxyscaled=Basics$proxyscaled; # inst=Basics$inst;instrumental.info=Basics$instrumental.info;instr.target=Basics$instr.target; # instr.smooth=Basics$instr.smooth # details=Basics$details ##WORKING OBJECTS # tempi - logical vector of length 1209 picking out proxies in step under analysis # sd.proxy.long - std deviation for proxy (not proxys) over period under study (not calibration period) # sd.proxy - std deviation for proxys in calibration period # outerlist1000 - this is subset list of associated gridcells for the proxies in use in the step. (22 for AD1000 step). #names are the proxy number # idsort - the associated gridcells in step (20 in AD1000). Obtained from outerlist1000 # Data - data frame of properties of all proxy-gridcell pairs # regid - list of regrid-gridcells used in step (15 in AD1000 step) ########################## #UTILITY FUNCTIONS ################################## ## plot 12 items plotf=function( X,col0=rep(1,12),title0="") { nf=layout( array( 1:12,dim=c(4,3)),heights=c(1.1,1,1,1.3),widths=c(1.2,1,1.1)) lab2= c(rep(TRUE,4),rep(FALSE,8)) lab1= rep(c(rep(FALSE,3),TRUE),3) ylim0=c(-3,5) par0= array(c( c(0,3,2,0), rep(c(0,3,0,0),2),c(3,3,0,0), c(0,0,2,0), rep(c(0,0,0,0),2),c(3,0,0,0), c(0,0,2,1), rep(c(0,0,0,1),2),c(3,0,0,1)), dim=c(4,12) ) for (i in 1:12) { par(mar=par0[,i]) x= window( X[,i] ,start=800,end=1995) plot(c(time(x)), scale( x),type="l",xlim=c(800,2010),ylim=ylim0,axes=FALSE,ylab="",xlab="",col=col0[i]) axis(side=1, labels=lab1[i]); axis(side=2,labels=lab2[i]); box() text(800,4,paste(dimnames(X)[[2]][i]),cex=.8,font=2,pos=4) if(i ==5) mtext(side=3,title0,line=.7,font=2) } } source("http://www.climateaudit.info/scripts/mann.2008/utilities.v4.txt") #nov 7, 20099 version #includes download of wrapper data (Basics) library(GDD) library(signal) ########################################### #####BENCHMARK NO TILJ NO DENDRO STEP ######################################### #SET PARAMETers Info=Basics passing=Basics$passing; names(passing) #[1] "whole" "earlym" "latem" outerlist=Info$outerlist$mbh smoothmethod="mann" lat.adjustment= -1 #weird displacement: see Stupid Pet Tricks screenmethod="mann" verbose="default" proxy=Info$proxy; info=Info$instrumental.info; period_MBH=Info$period_MBH annual=Info$criteria$annual; sum(annual) #1158 alias=c("r1850_1995","r1896_1995","r1850_1949");names(alias)=c("whole","earlym","latem") tempb= !is.na(match(Info$details$code,c(5000,5001,6000,6001) ) );sum(tempb) #33 codetable=Info$codetable Annual=list() #wrapper to collect averages ###################### ##PROXIES ############################## #tempi - selects proxies for step # k=9 case="latem"; calibration=Info$cal[[paste(case)]][1]:Info$cal[[paste(case)]][2] if(! case=="whole") verification=Info$ver[[paste(case)]][1]:Info$ver[[paste(case)]][2] name0=c("id","code","target","benchmark") nodendro= ! ((Info$details$code==9000)|(Info$details$code==3000) ); sum(nodendro) #278 notilj = is.na(match( 1:1209,1061:1064) ) criterion= notilj&nodendro ;sum(criterion) #274 target_id="iHAD" tempi=temp0= !is.na(proxy[period_MBH[k]-tsp(proxy)[1]+1,])&criterion; sum(temp0) ##29 #FLIPPING #collate correlations for case as details$target details=Info$details proxyw=Basics$proxy details$target= round(Info$details[,paste("rtable.",alias[paste(case)],sep="") ],3) details$target[!annual]= round(details[,paste("rtable.",alias[paste(case)],"lf",sep="")][!annual],3) #this collated LF and HF according to their designation details$benchmark= codetable[match(details$code,codetable$code),"EA100"] #this collates Mannian benchmark for each proxy depending on type, period #see collation.cps for construction of this table sum (details$target<0) #604 details$flip=flip= coral| (tempb& details$target<0) #flip proxy and target in details proxyw[,flip]= - proxyw[,flip] details$target[flip]= -details$target[flip] #calculate "passing" details$pass=pass = details$target>details$benchmark; sum(pass) #441 c(sum(flip),sum(flip&tempi)) #32 5 #MAKE PROXY DATA FRAME FOR SUBSET proxydata=details[tempi,c("id","code","lat","flip","target","benchmark","pass")] proxydata$proxy=(1:1209)[tempi] proxydata=proxydata[,c(8,1:7)] proxydata$code=as.numeric(as.character(proxydata$code)) proxydata$annual= proxydata$code%%2 ==0 dim(proxydata) #29 9 tapply(!is.na(proxydata$code),proxydata$code,sum) #calculate distribution #3001 4000 4001 5000 5001 6000 6001 8000 8001 # 1 2 8 1 1 4 5 4 3 temp= proxydata$pass tapply(!is.na(proxydata$code[temp]),proxydata$code[temp],sum) # 4001 5000 5001 6000 6001 8000 8001 # 3 1 1 4 2 4 1 #flipped in network details[flip&tempi,name0] # id code target benchmark #202 burns_2003_socotrad13c 6001 -0.5025322 0.521 #203 burns_2003_socotrad18o 6001 -0.5800865 0.521 #385 dongge 6001 -0.7834777 0.521 #481 lee_thorpe_2001_c13 6000 -0.2022475 0.164 #897 qian_2003_yriver 5000 -0.2049604 0.164 #MAKE STANDARDIZED WORKING (FLIPPED) SUBSET IN STEP tempp=(time(proxyw) >=period_MBH[k])&(time(proxyw)<=1995) working=ts(proxyw[tempp,tempi],start=period_MBH[k]);dim(working) #996 25 sd.proxy.long=apply(working,2,sd,na.rm=T);names(sd.proxy.long)=dimnames(proxy)[[2]][tempi]; mean.long= apply(working,2,mean,na.rm=T) #std deviation over full period of step sd.proxy.long # 192 202 203 204 205 378 #18.490514970 0.522163354 0.463948918 1.921210567 0.631206125 3.956891241 working= scale(working,center=mean.long,scale=sd.proxy.long) ;dim(working) #[1] 996 29 #FIugre 1a and 1b # GDD(file="d:/climate/images/2010/mann08/figure1a.gif",h=540,w=540,type="gif") col1=rep(1,nrow(proxydata)) col1[!proxydata$pass]=3 title1= "AD1000 Notilj Nodendro Unflipped Standardized" plotf(working[,proxydata$lat>0][,1:12],title0=title1,col0 =col1[proxydata$lat>0][1:12] ) # dev.off() # GDD(file="d:/climate/images/2010/mann08/figure1b.gif",h=540,w=540,type="gif") col1=rep(1,nrow(proxydata)) col1[!proxydata$pass]=3 title1= "AD1000 Notilj Nodendro Unflipped Standardized" plotf(working[,proxydata$lat>0][,13:22],title0=title1,col0 =col1[proxydata$lat>0][13:22] ) # dev.off() ## MAKE SUBSET OF PASSING PROXIES proxyw= working[, proxydata$target>proxydata$benchmark]; dim(proxyw) #996 15 proxydata=proxydata[ proxydata$target>proxydata$benchmark,] nh=proxydata$lat>0; sum(nh) #12 # Figure 2: selected standardized proxies - selected and not smoothed col1=rep(1,16); col1[proxydata$flip]=2 # GDD(file="d:/climate/images/2010/mann08/figure2.gif",h=540,w=540,type="gif") title1="AD1000 Notilj Nodendro Flipped - Unsmoothed" plotf (proxyw[,nh],col0=col1[nh] ,title0=title1) # dev.off() tapply(!is.na(proxydata$code),proxydata$code,sum) #calculate distribution #3001 4000 4001 5000 5001 6000 6001 8000 8001 # 1 2 8 1 1 4 5 4 3 ##SMOOTH THEN RE-SCALING proxys=proxyw cutfreq=.05;ipts=10 #ipts set as 10 in Mann lowpass bf=butter(ipts,2*cutfreq,"low");npad=1/(2*cutfreq);npad proxys[,!proxydata$annual]= apply(proxyw[,!proxydata$annual],2,function(x) mannsmooth(x,M=npad,bwf=bf )) cutfreq=.1;ipts=10 #ipts set as 10 in Mann lowpass bf=butter(ipts,2*cutfreq,"low");npad=1/(2*cutfreq);npad proxys[,proxydata$annual]=apply(proxyw[,proxydata$annual],2,function(x) mannsmooth(x,M=5,bwf=bf ) ) #this finishes piecemeal M08 smooth on truncated data dim(proxys) #996 25: butterworth smooth on Mannian #Figure 3 # GDD(file="d:/climate/images/2010/mann08/figure3.gif",h=540,w=540,type="gif") title1="AD1000 Notilj Nodendro Flipped - Smoothed" plotf (proxys[,nh],col0=col1[nh],title0=title1 ) # dev.off() ##ONCE-FOR-ALL SMOOTH ALSO AVAILABLE AS OPTION # if(smoothmethod=="sensible") { # # tempp=(time(proxy) >=period_MBH[k])&(time(proxy)<=1995) # proxys=ts(Basics$proxys[tempp,tempi],start=period_MBH[k]);#dim(working) #997 25 # sd.proxy.long=apply(proxy[tempp,tempi],2,sd,na.rm=T);names(sd.proxy.long)=dimnames(proxy)[[2]][tempi]; # #this is alternative using once-for-all smooth # #flipping done only once # ##RESCALE temp=!is.na(match(time(proxys),calibration)) ;#sum(temp) #his is calibration period mean.proxy=apply(proxys[temp,],2,mean,na.rm=T);names(mean.proxy)=dimnames(proxys)[[2]]; proxydata$sd.proxy=sd.proxy=round(apply(proxys[temp,],2,sd,na.rm=T),3) names(sd.proxy)=dimnames(proxys)[[2]]; range(sd.proxy) # $[1] 0.3343957 1.3389710 proxyrs = ts( scale(proxys,center=mean.proxy,scale=sd.proxy),start=tsp(proxys)[1]) Annual$proxyrs=ts(apply( proxyw[,proxydata$lat>0],1,mean,na.rm=T),start=tsp(proxyw)[1]) #Figure 4 - smoothed and re-scaled # GDD(file="d:/climate/images/2010/mann08/figure4.gif",h=540,w=540,type="gif") title1="Flipped Smoothed Short-scaled" plotf (proxyrs[,nh],title0=title1 ) # dev.off() ########################################################### #make data.frame Data with all proxy-gridcell combinations ########################################################## tempi= tempi& pass; sum(tempi) #16 outerlist1000= sapply(outerlist.mbh[tempi], function(x) instrumental.info$mannid[match(x,instrumental.info$jones)] ) # this yields a list in which proxies are linked to associated gridcells (can be more than one gridcell per proxy) length(outerlist1000) #16 idsort=sort(unique(c(unlist(outerlist1000))));#idsort #this is MAnn id (not jones id) length(idsort) #associated gridcells #[1] 14 #data frame with gridcells and regrid data frame Grid= data.frame(mannid= sort(unique(unlist(outerlist1000)))) Grid$mannid=as.numeric(as.character(Grid$mannid)) Grid$jones=info$jones[match(Grid$mannid,info$mannid)] #looks up Jones gridcell for Mann gridcell Grid$lat= info$lat[Grid$mannid]+lat.adjustment Grid$regrid=info$regrid[match(Grid$jones,info$jones)] Grid$regrid.lat=info$regrid.lat[Grid$mannid]+lat.adjustment dim(Grid) #14 5 #12 NH gridcells #data frame for Regrid gridcells regid=sort(unique(Grid$regrid)); K=length(regid);K #10; Regrid=data.frame(regid=regid) Regrid$regid=as.numeric(as.character(Regrid$regid)) regid # 15 #the set of unique regrid cell; and count Regrid$regrid.lat=Grid$regrid.lat[match(Regrid$regid,Grid$regrid)] dim(Regrid) #10 2 # 8 NH regrid- cells #data.frame for proxy-gridcell combinations - this is used to implement Mannian stupid pet trick of one proxy to more than one gridcell Data=NULL for(i in 1:length(outerlist1000)) Data=rbind(Data,data.frame(rep(names(outerlist1000)[[i]],length(outerlist1000[[i]])),outerlist1000[[i]])) names(Data)=c("proxy","mannid") nrow(Data) #28 #this is a matrix with all proxy-gridcell combinations for step for (i in 1:2) Data[,i]=as.numeric(as.character(Data[,i])) Data$jones=info$jones[match(Data$mannid,info$mannid)] #looks up Jones gridcell for Mann gridcell # mannid mann lat long jones mu regrid regrid.lat regrid.long mean_whole sd_whole mean_earlym sd_earlym mean_latem #1 1 NA -87.5 -177.5 2521 0.044 1 -87.5 -177.5 -0.196 0.240 -0.172 0.269 -0.296 #2 2 NA -87.5 -172.5 2522 0.044 2 -87.5 -172.5 -0.196 0.240 -0.172 0.269 -0.296 Data$gridcount=subcount(Data$proxy,Data$jones) #count the number of proxies in each gridcell Data$sd.proxy=proxydata$sd.proxy[match(Data$proxy,proxydata$proxy)] Data$regrid=info$regrid[match(Data$jones,info$jones)] #this looks up the regrid-cell for each gridcell (Mann merges NH gridcells) Data$lat=info$lat[Data$mannid]+lat.adjustment Data$regrid.lat=info$regrid.lat[Data$mannid]+lat.adjustment #this uses the regrid.lat for later weighting Data[1:2,] # proxy mannid jones gridcount regrid lat regrid.lat #1 192 2123 467 1 1764 56.5 51.5 #2 192 2124 468 1 1764 56.5 51.5 #GRID AVAILABLE PROXIES FOR UTLIZED GRIDCELLS (NOT regridded) grids=array(NA,dim=c(nrow(working),length(idsort) ) ) #length idsort dimnames(grids)[[2]]=idsort average= function(X) if( is.null(ncol(X) )) X else apply(X,1,mean,na.rm=T) for (i in 1:length(idsort) ) { tempr= !is.na(unlist(lapply(outerlist1000,function (A) match(Grid$mannid[i],A)) ) ) # ;sum(tempr) #this picks out the relevant series in gridcell grids[,i]=average( proxyrs[,tempr]) } grids=ts(grids,start=period_MBH[k]) dim(grids) #996 14 Annual$grids=ts(apply(grids[,Grid$lat>0],1,mean,na.rm=T),start=tsp(proxyw)[1]) #FIgure 5 - Gridded # GDD(file="d:/climate/images/2010/mann08/figure5.gif",h=540,w=540,type="gif") title1="Flipped Smoothed Short-scaled Gridded" plotf (X=grids[,Grid$lat>0],title0=title1 ) # dev.off() Data$weight= Data$sd.proxy/Data$gridcount proxydata$weight = round(tapply(Data$weight,Data$proxy,sum),2) #RESCALE GRIDCELL AVERAGE TO GRIDCELL INSTRUMENTAL TARGET mean.obs= info[idsort,paste("mean",case,sep="_")]; #this locates right match in instrumental.info sd.obs=info[idsort,paste("sd",case,sep="_")]; names(sd.obs)=idsort;sd.obs # 670 714 719 978 1102 1486 1487 1603 1714 1787 1858 # 0.3796854 0.5138550 0.2149656 0.3252695 0.5545764 0.1795864 0.2258918 0.2659750 0.2243348 0.3885594 0.3841052 mean.est=apply(grids[calibration-tsp(grids)[[1]]+1,],2,mean);range(mean.est) #-2.758761e-15 2.145853e-15 sd.est=apply(grids[calibration-tsp(grids)[[1]]+1,],2,sd);range(sd.est);sd.est # .6405801 1.0000000 # 670 714 719 978 1102 1486 1487 1603 1714 1787 1858 #0.8556082 1.0000000 1.0000000 0.7496158 0.6405801 1.0000000 1.0000000 0.7468546 1.0000000 1.0000000 1.0000000 regts=array(NA,dim=c(nrow(grids),ncol(grids)) ); dimnames(regts)[[2]]=dimnames(grids)[[2]] for( i in 1:ncol(grids)) { regts[,i]=(grids[,i]-mean.est[i])*sd.obs[i]/sd.est[i]+mean.obs[i] } regts=ts(regts,start=tsp(grids)[1]);dim(regts) #996 14 Grid$regts.weight=regts.weight=round(sd.obs/sd.est,3); regts.weight # 978 1102 1486 1487 1603 1714 1787 1858 1859 1932 2123 2124 2333 2470 #0.418 0.981 0.180 0.226 0.344 0.224 0.389 0.384 0.342 0.421 0.257 0.227 0.777 0.391 Data$weight=Data$weight* sd.obs[paste(Data$mannid)]/sd.est[paste(Data$mannid)] proxydata$weights= round(tapply(Data$weight,Data$proxy,sum),2) Annual$regts=ts(apply(regts[,Grid$lat>0],1,mean,na.rm=T),start=tsp(proxyw)[1]) layout(1); par(mar=c(3,4,2,1)); plot.ts(Annual$regts) #FIgure 6 - Gridded Rescaled # GDD(file="d:/climate/images/2010/mann08/figure6.gif",h=540,w=540,type="gif") title1="Gridded Rescaled" plotf (X=regts[,Grid$lat>0],title0=title1 ) # dev.off() #REGRID N of 30N into 15x15 gridcells weighted.average=function(X,w) if( is.null(ncol(X) )) X else apply(X,1,function(x) weighted.mean(x,w,na.rm=T) ) rescaled=array(NA,dim=c(nrow(regts),length(regid))) #15 dimnames(rescaled)[[2]]=Regrid$regid;K=length(Regrid$regid) for(i in 1:K) { tempr= !is.na(match(info$regrid[Grid$mannid],dimnames(rescaled)[[2]][i])) rescaled[,i]=weighted.average( regts[,tempr],info$mu[idsort[tempr]]) } rescaled=ts(rescaled,start=tsp(regts)[1],end=1995) #ends one year early #996 15 dim(rescaled) #996 10 Annual$rescaled=ts(apply(rescaled[,Regrid$regrid.lat>0],1,mean,na.rm=T),start=tsp(proxyw)[1]) Grid$regrid.count= subcount(Grid$mannid,Grid$regrid,f=function(x) length(unique(x))) #this counts the number of gridcells in each regrid cell and associates count with each line #Figure 7 Rescaled and Regridded # GDD(file="d:/climate/images/2010/mann08/figure7.gif",h=540,w=540,type="gif") title1="Gridded Rescaled Regridded" plotf(rescaled[,Regrid$regrid.lat>0]) # dev.off() #MAKE NH and SH COMPOSITES #weight by cos(latitude) nhtemp=Regrid$regrid.lat>0;sum(nhtemp) Regrid$mu=regrid.weights=cos(Regrid$regrid.lat*pi/180);names(regrid.weights)=regid x= Regrid$mu[Regrid$regrid.lat>0]; x=x/sum(x); x composite= ts( array(NA,dim=c(nrow(rescaled),2)),start=period_MBH[k]) dimnames(composite)[[2]]=c("NH","SH") composite[,"NH"]=apply(rescaled[,nhtemp],1,function(x) weighted.mean(x,regrid.weights[nhtemp],na.rm=T) ) composite[,"SH"]=apply(rescaled[,!nhtemp],1,function(x) weighted.mean(x,regrid.weights[!nhtemp],na.rm=T) ) Annual$composite=composite[,"NH"] #Figure 9 Composites # GDD(file="d:/climate/images/2010/mann08/figure9.gif",h=360,w=540,type="gif") nf=layout(array(1:2,dim=c(2,1)),heights=c(1.1,1.3)) ylim0=c(-1.2,1.2) par(mar=c(0,3,2,1)) plot(c(time(composite)),composite[,"NH"], type="l", axes=FALSE, xlab="",ylab="",ylim=ylim0) axis(side=1,labels=FALSE); axis(side=2,las=1); box(); abline(h=0,lty=3) title("AD1000 Notilj Nodendro CPS in M08 STyle") par(mar=c(3,3,0,1)) plot(c(time(composite)),composite[,"SH"], type="l", axes=FALSE, xlab="",ylab="",ylim=ylim0) axis(side=1); axis(side=2,las=1); box(); abline(h=0,lty=3) # dev.off() ##ADDED TO BLOG SCRIPT # CALCULATE CPS RECONSTRUCTIONS TO TARGET x= Info$cal[[paste(case)]] target_id #iHAD cps0= cpsf(composite,target_id,index=x[1]:x[2] ) cps0[1,] # NH SH GL # 0.1468016 -0.8838390 -0.3685187 (stat=verification.stats(observed= instr.smooth[,paste(target_id,"NH",sep="_")],estimator= cps0[,"NH"],calibration=Info$cal[[paste(case)]], verification=Info$ver[[paste(case)]] )) # RE.cal RE.ver R2.cal R2.ver CE #1 0.5579708 -0.7294076 0.6068183 0.4574853 -6.701885