Skip to content

@nealmcb /replicate_trick.r
Created

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
## Climate analysis R code based on Steve McIntyre's example at
## http://camirror.wordpress.com/2009/11/29/replicating-the-trick-diagram/#more-130
## with some fixes by Neal McBurnett.
## Usage: R --save < replicate_trick.r
## which produces an "Rplots.ps" file
##COMPARE ARCHIVED BRIFFA VERSION TO CLIMATEGATE VERSION#
#1. LOAD BRIFFA (CLIMATEGATE VERSION)
# archive is truncated in 1960: ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/reconstructions/n_hem_temp/briffa2001jgr3.txt”
loc="http://www.eastangliaemails.com/emails.php?eid=146&filename=939154709.txt"
working=readLines(loc,n=1994-1401+109)
working=working[110:length(working)]
x=substr(working,1,14)
writeLines(x,"temp.dat")
gate=read.table("temp.dat")
gate=ts(gate[,2],start=gate[1,1])
#2. J98 has reference 1961-1990
#note that there is another version at ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/jones1998/jonesdata.txt"
loc="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/jones2001/jones2001_fig2.txt"
test=read.table(loc,skip=17,header=TRUE,fill=TRUE,colClasses="numeric",nrow=1001)
test[test== -9.999]=NA
count= apply(!is.na(test),1,sum)
test=ts(test,start=1000,end=2000)
J2001=test[,"Jones"]
#3. MBH : reference 1902-1980
url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/mann1999/recons/nhem-recon.dat"
MBH99<-read.table(url) ;#this goes to 1980
MBH99<-ts(MBH99[,2],start=MBH99[1,1])
#4. CRU instrumental: 1961-1990 reference
# use old version to 1997 in Briffa archive extended
url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/reconstructions/n_hem_temp/briffa2001jgr3.txt"
#readLines(url)[1:50]
Briffa<-read.table(url,skip=24,fill=TRUE)
Briffa[Briffa< -900]=NA
dimnames(Briffa)[[2]]<-c("year","Jones98","MBH99","Briffa01","Briffa00","Overpeck97","Crowley00","CRU99")
Briffa= ts(Briffa,start=1000)
CRU=window(Briffa[,"CRU99"],start=1850)
tsp(CRU) # 1850 1999 #but starts 1871 and ends 1997
delta<-mean(CRU[(1902:1980)-1850])-mean(CRU[(1960:1990)-1850]);
delta # -0.118922
#used to get MBH values with 1961-1990 reference: compare to 0.12 mentioned in Climategate letters
#get updated version of CRU to update 1998 and 1999 values
loc="http://hadobs.metoffice.com/crutem3/diagnostics/hemispheric/northern/annual"
D=read.table(loc) #dim(D) #158 12 #start 1850
names(D)=c("year","anom","u_sample","l_sample","u_coverage","l_coverage","u_bias","l_bias","u_sample_cover","l_sample_cover",
"u_total","l_total")
cru=ts(D[,2],start=1850)
tsp(cru) # 1850 2009
# update 1998-1999 values with 1998 values
CRU[(1998:1999)-1849]= rep(cru[(1998)-1849],2)
#Fig 2.21 Caption
#The horizontal zero line denotes the 1961 to 1990 reference
#period mean temperature. All series were smoothed with a 40-year Hamming-weights lowpass filter, with boundary constraints
# imposed by padding the series with its mean values during the first and last 25 years.
#this is a low-pass filter
source("http://www.climateaudit.info/scripts/utilities.txt") #get filter.combine.pad function
hamming.filter<-function(N) {
i<-0:(N-1)
w<-cos(2*pi*i/(N-1))
hamming.filter<-0.54 - 0.46 *w
hamming.filter<-hamming.filter/sum(hamming.filter)
hamming.filter
}
f=function(x) filter.combine.pad(x,a=hamming.filter(40),M=25)[,2]
Y=X=ts.union(MBH=MBH99+delta,J2001,briffa=gate,CRU=window(CRU,end=1999) ) #collate
temp= time(Y)>1960
Y[temp,"briffa"]=Y[temp,"CRU"]
temp= time(Y)>1980
Y[temp,c("MBH","J2001")]=Y[temp,"CRU"]
smoothb= ts(apply(Y,2,f),start=1000)
xlim0=c(1000,2000) #xlim0=c(1900,2000)
ylim0=c(-.6,.35)
par(mar=c(2.5,4,2,1))
col.ipcc=c("blue","red","green4","black")
par(bg="beige")
plot( c(time(smoothb)),smoothb[,1],col=col.ipcc,lwd=2,bg="lightblue",xlim=xlim0,xaxs="i",ylim=ylim0,yaxs="i",type="n",axes=FALSE,xlab="",ylab="deg C (1961-1990)")
usr <- par("usr")
rect(usr[1],usr[3],usr[2] ,usr[4],col="lightblue") # – the part used to fit the model
for( i in 1:3) lines( c(time(smoothb)),smoothb[,i],col=col.ipcc[i],lwd=2)
axis(side=1)
labels0=labels1=seq(-.6,.4,.1);
labels0[is.na(match(seq(-.6,.4,.1),seq(-.6,.4,.2)))]="";labels0[7]="0"
axis(side=2,at=labels1,labels=labels0,tck=.025,las=1)
axis(side=4,at=labels1,labels=labels0,tck=.025)
box()
abline(h=0)
legend("topleft", fill= col.ipcc[c(2,1,3)],
legend=c( "Jones 1998", "Mann 1999", "Briffa 2000"))
title("WMO 1999 Emulation")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.