##GERGIS CORRELATIONS #used in https://climateaudit.org/2016/08/03/gergis-and-law-dome/ #source("http://www.climateaudit.info/scripts/multiproxy/gergis_2016/gergis_cor_final.txt") #print(tag(B[order(abs(B$rcalc),decreasing=T),][1:10,]) ##LOAD DATA ################ #PROXY loc="http://www.climateaudit.info/data/multiproxy/gergis_2016/proxy_gergis_2016.tab" dest="d:/temp/temp.dat" download.file(loc,dest,mode="wb") load(dest) dim(proxy)#1003 51 id=dimnames(proxy)[[2]] loc="http://www.climateaudit.info/data/multiproxy/gergis_2016/info_gergis_2016.csv" info=read.csv(loc) dim(info) # 51 44 ind= match(id,info$id) info=info[ind,] info$id=factor(info$id) #INSTRUMENTAL #library(ncdf) #loc="http://www.metoffice.gov.uk/hadobs/hadcrut3/data/HadCRUT3.nc" #download.file(loc,"d:/temp/temp.nc",mode="wb") #v=open.ncdf("d:/temp/temp.nc") #instr <- get.var.ncdf( v, v$var[[1]]) # 1850 2006 #save(instr,file="d:/temp/instr.tab") #upload to www.climateaudit.info/data/multiproxy/gergis_2016/instr.tab loc="http://www.climateaudit.info/data/multiproxy/gergis_2016/instr.tab" download.file(loc,"d:/temp/instr.tab",mode="wb") load("d:/temp/instr.tab") had=instr[,36:1,] Long=seq(-177.5,177.5,5) Lat=seq(87.5,-87.5,-5) ##SUMMER had[37,8,][9:14] #-2.3 -3.7 0.8 -0.3 1.7 -0.4 #extract Sep-August work=array(had[,,(1:1956)+8],dim=c(72,36,12,163)) work[37,8,,1][1:6] # -2.3 -3.7 0.8 -0.3 1.7 -0.4 #take summer average where available: other options with missing data available target=work[,,1:6,] summer= apply(target,c(1,2,4),mean,na.rm=T) # dim(summer) #[1] 72 36 163 ##FUNCTIONS require(lmtest) #detrend detrend=function(x) { temp=!is.na(x) y=rep(NA,length(x)) y[temp]=lm(x~ c(time(x)))$residuals return(y) } #cell coordinate latg= function(x) 5*floor(x/5)+2.5 circledist =function(x,R=6372.795) { #fromlat,fromlong, lat,long pi180=pi/180; y= abs(x[4] -x[2]) y[y>180]=360-y[y>180] y[y<= -180]= 360+y[y<= -180] delta= y *pi180 fromlat=x[1]*pi180; tolat=x[3]*pi180; tolong=x[4]*pi180 theta= 2* asin( sqrt( sin( (tolat- fromlat)/2 )^2 + cos(tolat)*cos(fromlat)* (sin(delta/2))^2 )) circledist=R*theta circledist } use0="pairwise.complete.obs" tag=function(A) { B=data.frame(A); row.names(B)=paste("#",row.names(B)); return(B) } ##CALCULATIONS #################### #For predictor selection, both proxy climate and instrumental data were linearly detrended over 1931-90. As detailed in appendix A, only records that were significantly (p <0.05) correlated with temperature variations in at least one grid cell within 500km of the proxy’s location over the 1931-90 period were selected for further analysis. For predictor selection, both proxy climate and instrumental data were linearly detrended over 1931-90…To account for proxies with seasonal definitions other than the target SONDJF season (e.g., calendar year averages), the comparisons were performed using lags of -1, 0, and +1 years for each proxy… #The significance of these correlations was determined using a Student’s t distribution, with the degrees of freedom adjusted for autocorrelation at lag 1 using the method of Bretherton et al. (1999) #Bretherton, C. S., M. Widmann, V. P. Dymnikov, J. M. Wallace, and I. Bladé, 1999: The effective number of spatial degrees of freedom of a time-varying field. #This comparison was only performed for cells containing at least 50 years of data between 1921 and 1990. #DETREND PROXY OVER CALIBRATION PERIOD #set start and end to correspond to Gergis 2016 calibration end0=1990 start0=1931 proxy_det=window(proxy,start0,end0) for ( i in 1:51) proxy_det[,i]=detrend(proxy_det[,i]) #DETRENDED r and t STATISTICS FOR EACH PROXY Cell=NULL for(i in 1:51) { cell= data.frame(lat=rep( c(latg(info$lat[i])+c(-5,0,5)),each=3), long=rep(latg(info$long[i])+c(-5,0,5),3),dist=rep(NA,9)) cell$dist=apply(cell[,1:2],1, function(z) circledist(c(info$lat[i],info$long[i],z[1],z[2])) ) cell$num=1:9 work= cell[cell$dist<=500,] K=nrow(work) work=rbind(work,work,work) work$lag= rep( c(-1,0,1),each=K) work$i= match(work$long,Long) work$j= match(work$lat,Lat) work$icount=NA work$r=NA work$neff=NA work$t=NA work$dw=NA work$dwp=NA work$rain=NA work$rainp=NA for(j in 1: (3*K)) { x= ts(summer[work$i[j],work$j[j],],start=1851) L=work$lag[j] A= window(ts.union(proxy= lag(proxy_det[,i],k=-L), instr=x,detrend=NA),start0+L,end0+L) work$icount[j]=sum(!is.na(A[,"instr"])) if(work$icount[j] >=5) { A[,"detrend"]=detrend(A[,"instr"]) fm=lm(A[,"proxy"]~A[,"detrend"]) N=length(fm$residuals) se.ols=summary(fm)$coef[2,2] beta= fm$coef[2] #bretherton calculation r1=arima(A[,"proxy"],order=c(1,0,0))$coef[1] r2=arima(A[,"detrend"],order=c(1,0,0))$coef[1] work$r[j]=round(cor(A[,1],A[,3],use=use0),3) work$neff[j]= round(N* (1-r1*r2)/(1+r1*r2),1) #Bretherton formula work$t[j] =round(beta/ ( sqrt(( N-2)/(work$neff[j]-2))* se.ols),3); #Sep 2016: was error in this: formerly neff instead of work$neff[j] test=dwtest(fm) work[j,c("dw","dwp")]= round(unlist(test[c("statistic","p.value")]),4) test=raintest(fm) work[j,c("rain","rainp")]=round(unlist(test[c("statistic","p.value")]),4) } #if } #j work$dset=info$id[i] Cell=rbind(Cell,work) }#i Cell$dset=factor(Cell$dset,levels=info$id[1:51]) tapply(!is.na(Cell$num),Cell$dset, sum)[1:3] # X669_MtRead_TT_recon X588_Law_Dome_18O_new X1075_Oroko_TT_recon # 9 12 6 ##EXTRACT CENTRAL CELL AND SHOW DESCENDING Cellz=Cell[Cell$num==5,] info$lag_calc=NA info$tcalc=NA info$rcalc=NA for(i in 1:51){ temp= Cellz$dset==info$id[i] A=Cellz[temp,] if(sum(!is.na(A$t))>0) { info$lag_calc[i]=A$lag[abs(A$t)==max(abs(A$t),na.rm=T)] info$rcalc[i]=A$r[abs(A$t)==max(abs(A$t))] info$tcalc[i]=A$t[abs(A$t)==max(abs(A$t))]} } B=info[!is.na(info$rcalc),c("id","lag_calc","rcalc","tcalc","G16")] print(tag(B[order(abs(B$rcalc),decreasing=T),][1:10,])) # id lag_calc rcalc tcalc G16 # 4 X1091_Palmyra_d18O_floating 0 -0.766 -8.597 1 # 30 X1076_Fiji_1F_SrCa 0 -0.705 -7.201 1 # 16 X607_Fiji_AB_d18O 0 -0.646 -6.201 1 # 21 X850_New_Caledonia_d18O -1 -0.536 -4.006 1 # 2 X588_Law_Dome_18O_new 0 0.529 3.652 0 # 8 X897_New_Fenwick_Composite_Signalfree_2 1 0.525 4.296 1 # 32 X600_Fiji_1F_d18O 0 -0.507 -4.434 0 # 22 X937_Stewart_Island_HABI_composite_Signalfree_2 0 0.502 4.041 1 # 9 X885_URW_newz063_Signalfree_2 0 -0.490 -3.933 1 # 13 X667_All.Celery.Top.west.corr_Signalfree_2 -1 0.440 3.726 1