Skip to content

Instantly share code, notes, and snippets.

/gist:2813b0035373ddd0dcad6991b92a4d50

Created September 12, 2016 00:42
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save anonymous/2813b0035373ddd0dcad6991b92a4d50 to your computer and use it in GitHub Desktop.
Save anonymous/2813b0035373ddd0dcad6991b92a4d50 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
gistfile1.txt
##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="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"
#loc="http://www.metoffice.gov.uk/hadobs/hadcrut3/data/HadCRUT3v.nc"
download.file(loc,"temp/temp.nc",mode="wb")
v=open.ncdf("temp/temp.nc")
instr <- get.var.ncdf( v, v$var[[1]]) # 1850 2006
save(instr,file="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,"temp/instr.tab",mode="wb")
load("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:28,]))
# Results using the variance-adjusted data set:
id lag_calc rcalc tcalc G16
# 4 X1091_Palmyra_d18O_floating 0 -0.776 -8.834 1
# 30 X1076_Fiji_1F_SrCa 0 -0.700 -7.054 1
# 16 X607_Fiji_AB_d18O 0 -0.640 -6.069 1
# 26 X585_Vostok_d18O -1 0.545 3.418 0
# 8 X897_New_Fenwick_Composite_Signalfree_2 1 0.530 4.291 1
# 22 X937_Stewart_Island_HABI_composite_Signalfree_2 0 0.521 4.156 1
# 32 X600_Fiji_1F_d18O 0 -0.507 -4.422 0
# 9 X885_URW_newz063_Signalfree_2 0 -0.480 -3.758 1
# 21 X850_New_Caledonia_d18O -1 -0.472 -3.163 1
# 13 X667_All.Celery.Top.west.corr_Signalfree_2 -1 0.445 3.783 1
# 3 X1075_Oroko_TT_recon 0 0.427 3.566 1
# 25 X869_Rarotonga_d18O 0 -0.404 -3.136 1
# 12 X879_TAK_newz062_Signalfree_2 0 -0.394 -2.386 0
# 6 X829_Law_Dome_Accumulation 0 0.384 2.306 0
# 45 X1017_Rarotonga.3R_d18O 0 -0.384 -2.888 1
# 24 X857_Rarotonga_SrCa 0 -0.380 -2.973 0
# 29 X913_Savusavu_d18O 0 -0.375 -2.577 1
# 39 X923_Maiana_d18O 0 -0.371 -2.910 0
# 7 X656_CTP_East_BLT_RFR_Signalfree_2 -1 0.370 3.040 1
# 2 X588_Law_Dome_18O_new 0 0.368 2.331 0
# 48 X848_Laing_d18O 0 -0.355 -2.755 1
# 50 X873_Tarawa 0 -0.347 -2.637 0
# 10 X650_BCH_Signalfree_2 1 0.346 2.802 1
# 42 X597_Bunaken_d18O -1 -0.345 -2.713 1
# 1 X669_MtRead_TT_recon 0 0.324 2.599 1
# 44 X976_Rabaul_SrCa 0 -0.304 -2.231 0
# 40 X1093_Moorea_d18O 0 -0.299 -1.927 0
# 38 X334_Avaiki_ASM1 1 0.295 2.278 1
# Results using the data set that is not variance-adjusted:
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
# 3 X1075_Oroko_TT_recon 0 0.432 3.621 1
# 39 X923_Maiana_d18O 0 -0.428 -3.457 0
# 6 X829_Law_Dome_Accumulation 0 0.427 2.647 0
# 50 X873_Tarawa 0 -0.414 -3.256 0
# 25 X869_Rarotonga_d18O 0 -0.410 -3.250 1
# 24 X857_Rarotonga_SrCa 0 -0.392 -3.127 0
# 45 X1017_Rarotonga.3R_d18O 0 -0.387 -2.984 1
# 7 X656_CTP_East_BLT_RFR_Signalfree_2 -1 0.368 3.010 1
# 29 X913_Savusavu_d18O 0 -0.360 -2.521 1
# 42 X597_Bunaken_d18O -1 -0.354 -2.794 1
# 10 X650_BCH_Signalfree_2 1 0.345 2.809 1
# 12 X879_TAK_newz062_Signalfree_2 0 -0.338 -2.250 0
# 1 X669_MtRead_TT_recon 0 0.317 2.546 1
# 33 X599_Bali_d18O -1 0.314 2.518 0
# 48 X848_Laing_d18O 0 -0.303 -2.358 1
# 40 X1093_Moorea_d18O 0 -0.302 -1.997 0
# 38 X334_Avaiki_ASM1 1 0.290 2.273 1
# 14 X839_MANGAWE_Signalfree_2 1 0.286 2.129 1
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment