## 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")