###SPAGHETTI GRAPH VERSIONS #collate spaghetti graph from original sources #makes spaghetti.txt #spaghetti<-ts.union(MBH99,J98,jones01,esper02,esper.fitted,CL15.emulation,CL2,CL2.Jns11,moberg,mj03,hegerl,juckes) #12 Direct black 1961-90 #1mbh green 1902-80 #2jones red 1961-90 #4Briffa blue 1961-90 #11 Esper 2002 fitted brown 1961-90 #3 Crowley yellow 1902-80 #6 moberg 2005 grey #10 Hegerl 2007 beige #9 Juckes darj grey 1961-90 #8 Darrigo 2006green4 1961-90 #7 Oerlemans 2005 brown 1961-90 #5 huang 2004 indigo 1961-90 #Wilson http://www.climateaudit.org/?p=1586#comment-109690 ##The time-series were smoothed with a 50 year cubic smoothing spline ##In the 2007 version, I scaled (same mean and variance) the smoothed proxy records to the smoothed instrumental data (annual extra-tropical land only temperatures) ## over the 1856-1960 period. scripts="d:/climate/scripts/spaghetti" source(file.path(scripts,"CRU3.txt")) source(file.path(scripts,"CRU3.txt")) source(file.path(scripts,"CRU3.txt")) source(file.path(scripts,"CRU3.txt")) source(file.path(scripts,"CRU3.txt")) source(file.path(scripts,"CRU3.txt")) source(file.path(scripts,"CRU3.txt")) source(file.path(scripts,"CRU3.txt")) source(file.path(scripts,"CRU3.txt")) source(file.path(scripts,"CRU3.txt")) source(file.path(scripts,"CRU3.txt")) source(file.path(scripts,"CRU3.txt")) ##12 CRU NH #potential target 1961-1990 url<-"http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3nh.txt" #1850 2007 h<-scan(url) start0<-h[1] N<-length(h) h<-array(h,dim=c(27,N/27)) h<-t(h) CRU<-ts(h[,14],start=start0,end=2006) #1851 - 2002 # 1 #MBH99 #1902-1980 url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/mann1999/recons/nhem-recon.dat" g<-read.table(url) #this goes to 1980 MBH99<-ts(g[,2],start=g[1,1]) delta.mbh<-mean(CRU[(1960:1990)-1850])-mean(CRU[(1902:1980)-1850]);delta.mbh # -0.1548959 #MBH99=MBH99-delta.mbh #2 JONEs 1998 #1961-1990 #no adjustment url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/jones1998/jonesdata.txt" g<-read.table(url,skip=48,fill=TRUE) dimnames(g)[[2]]<-c("year","NH.norm","NH.anom","NH.count","SH.norm","SH.anom","SH.count") J98<-ts(g[,3],start=1000,end=1991) #1961-1990 centered #The 1961-90 anomalies are related to the normalised series as follows: #nh*0.521 - 0.1134 #sh*0.608 - 0.2394 #The numbers come from the mean and variance of the instrumental data (summer, JJA in NH and DJF in SH) for the 1901-50 period. ##4 BRIFFA JGR 2001 SPAGHETTI DIAGRAM #1961-990 mean per WDCP preamble #this is scaled differently than original sources for non-Briffa 2001 versions url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/briffa1998/briffa2001jgr3.txt" #this has 20N version from Rev Geophys 1999 Apr-Sept #1961-990 mean Briffa<-read.table(url,skip=24,fill=TRUE) temp<-(Briffa< (-999)) Briffa[temp]<-NA dimnames(Briffa)[[2]]<-c("Year","Jones98","MBH99","Briffa01","Briffa00","Overpeck97","Crowley00","CRU99") Briffa<-ts(Briffa,start=Briffa[1,1]) #this gives series 1, 12,3, NA,NA,4,NA. #starts 1000; 1000; 1402; 1000; 1600;1000;1871 #ends 1991; 1980; 1960; 1993; 1990; 1987; 1997 #array from 1000 to 1999 #Briffa00 gores to 1993; Briffa 01 truncated in 1960 (I've emulated extension) briffa01<-ts(Briffa[!is.na(Briffa[,4]),4],start=1402) ##11 ESPER02 RECALIBRATION BY COOK ET AL 2004 #1961-90 by fitting to Briffa #this is a chronology and is not scaled to temperature #Northern Hemisphere Extratropical Regional Curve Standardization (RCS) Chronology #the Cook calibration fits to Briffa CRU basis 1961-1990 url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/reconstructions/n_hem_temp/esper2002_nhem_temp.txt" h<-read.table(url,skip=53,fill=TRUE) esper02<-ts(h[,2],start=h[1,1]) #831 1992 target<-Briffa[,"CRU99"] #already loaded with Briffa z<-ts.union(target,esper02); z<-data.frame(z) # 831 1999 fm<-lm(target~esper02,data=z[(1881:1960)-830,]) esper.fitted<-ts( predict(fm,newdata=z),start=831) #I have not been able to locate an original archive of this #3 CROWLEY 2000 ##CL2 and CL2.Jns11 spliced version #Crowley 2000: Hemisphere climate records and was scaled (12) to temperature using the instrumental record (16) in the overlap interval 1860–1965. url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo/gcmoutput/crowley2000/crowley_fig1_data.txt" crowley<-read.table(url,skip=1,sep="\t") dimnames(crowley)[[2]]<-c("Year","CL2.Jns11","Mn.sm11","CL2","+2s.d.adj","-2s.d.adj") temp<-!is.na(crowley[,2]); crowley<-crowley[temp,] #CL2<-ts(crowley[,"CL2.Jns11"],start=1000) CL2.Jns11<-ts(crowley[,"CL2"],start=1000) #CL2.Jns11 is a commonly used version in which the CRU temperature reconstruction has been spliced after 1880 #CL2 differs from CL00 version which is not archived anywhere #this appears to be centered on 1880-1960 delta.crowley=mean(CRU[(1961:1990)-1850])-mean(CRU[(1850:1960)-1850]);delta.crowley # 0.2600303 mean(CL2.Jns11[(1850:1960)-999]) #mean(CL2.Jns11[(1850:1960)-999]) ts.plot(CL2.Jns11) ##6 MOBERG 1961-1990 #SI at http://www.nature.com/nature/journal/v433/n7026/suppinfo/nature03265.html url<-"http://www.nature.com/nature/journal/v433/n7026/extref/nature03265-s6.doc" #Temperature anomalies in a, b and d are with respect to the 1961–90 average: Figure 2 caption recon<-read.table(url,skip=19,nrow=1998-19) #1979 9 #from 1 to 1979 name<-c("year","recon","lowf","lowf.lb","lowf.ub","ci.lb","ci.ub","ci2.lb","ci2.ub") dimnames(recon)[[2]]<-name recon<-ts(recon[,2:ncol(recon)],start=1,end=1979) temp<-(recon< -9) recon[temp]<-NA moberg<-recon[,"recon"] #mean(moberg[(1902:1980)-0]) ##HEGERL ET AL NATURE 2006 #this is archived in an Excel file ##source is "http://www.nature.com/nature/journal/v440/n7087/extref/nature04679-s2.xls" #I downloaded this and saved into txt file fred<-scan("d:/climate/data/hegerl.2006/nature04679-s2.txt",skip=1) #starts different at 1505 item 1017 dummy<-array(NA,dim=c(1960-1250,7)) dummy[1:(1504-1250),1:4]<-t(array(fred[1:1016],dim=c(4,1016/4))) dummy[(1505-1250):(1960-1250),]<-t(array(fred[1017:length(fred)],dim=c(7,(length(fred)-1016)/7))) hegerl<-ts(dummy[,2:7],start=1251) dimnames(hegerl)[[2]]<-c("long","plus","minus","blend","bplus","bminus") hegerl<-ts(hegerl[,1],start=1251) #mean(hegerl[(1880:1960)-1250]) ts.plot(fg(hegerl)) lines(1850:2006,fg(CRU),col="red") ts.plot(fg(MBH99),lwd=2,ylim=c(-.8,.4)) lines(1000:1991,fg(J98),col="red") lines(1402:1960,fg(briffa01),col="green4",lwd=2) lines(831:1994,fg(esper.fitted),col="cyan",lwd=2) ##JUCKES UNION #download.file("http://home.badc.rl.ac.uk/mjuckes/mitrie_files/data/mitrie_new_reconstructions_v01.nc","d:/climate/data/mitrie_new_reconstructions_v01.nc") #this is collated from Juckes' zipped NCF files http://www.clim-past-discuss.net/2/1001/2006/cpd-2-1001-2006-supplement.zip #url="http://home.badc.rl.ac.uk/mjuckes/mitrie_files/data/mitrie_new_reconstructions_v01.nc" #url="http://home.badc.rl.ac.uk/mjuckes/mitrie_files/data/mitrie_new_reconstructions_v01.csv" ##FROM NC VERSION url<-"d:/climate/data/mitrie" loc<-file.path(url,"newdata", id[6]) #"mitrie_new_reconstructions_v01.nc" ) #[1] fred<-open.ncdf(file.path(url,"newdata","mitrie_new_reconstructions_v01.nc" )) #"http://home.badc.rl.ac.uk/mjuckes/mitrie_files/data/mitrie_new_reconstructions_v01.nc N<-fred$nvars #64 juckes.recon<-array(NA,dim=c(986,N)) id<-rep(NA,N) for (i in 1:N) { v1 <- fred$var[[i]] x<- get.var.ncdf( fred, v1 ) # by default, reads ALL the data juckes.recon[,i] <-x id[i]<-v1$name } temp<-(juckes.recon== -999) juckes.recon[temp]<-NA dim(juckes.recon) #[1] 986 64 juckes.recon<-ts(juckes.recon,start=1000) dimnames(juckes.recon)[[2]]<-id juckes.union=juckes.recon[,1] ##ARCHIVED CA loc<-"http://www.climateaudit.org/data/mitrie" recon<-read.table(file.path(loc,"juckes.recon.txt"),sep="\t",header=TRUE) juckes<-ts(recon[,"union"],start=1000) combine=ts.union(juckes,juckes.union)##the CA archived can be recovered from first series of frm=function(x){frm=filter.combine.pad(x,rep(1/21,21))[,2];frm} ts.plot(frm(juckes),lwd=2);abline(h=0) lines(1000:1980,frm(MBH99),lwd=2,col="green") #Graphs have been smoothed with a 21-year running mean and centred on 1866 to 1970. The #maximum of the “Union” reconstruction in the pre-industrial period (0.25 K, 1091 AD) is shown #by a cyan cross, the maximum of the instrumental record (0.841 K, 1998 AD) is shown as a red #cross. ##DARRGIO url="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/reconstructions/n_hem_temp/nhtemp-darrigo2006.txt" #"YearAD NHLandObs STD recon RCS recon" #[1229] " 1856 -0.55 -0.55 -0.69" #[1230] " 1857 -0.38 -0.39 -0.48" #[1231] " 1858 -0.3 -0.19 -0.3" darrigo=read.fwf(url,skip=85,fill=TRUE,widths=c(6,14,15,13)) names(darrigo)=c("year","instr","STD","RCS") darrigo=ts(darrigo[!is.na(darrigo[,"year"],"RCS"),start=darrigo[1,1]) darrigo=ts(darrigo[!is.na(darrigo)],start=tsp(darrigo)[1]) ##OERLEMANS url="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/oerlemans2005/oerlemans2005.txt" oerl=read.table(url,skip=51,fill=TRUE) oerl=ts(oerl[,2],start=oerl[1,1]) ##HUANG 1961-1990 url="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/huang2004/nhtemp-huang2004.txt" huang=read.table(url,skip=61,fill=TRUE) huang=ts(huang[,2],start=huang[1,1]) #12 Direct black #1mbh green #2jones red #4Briffa blue #11 Esper 2002 brown #3Crowley yellow #6moberg 2005 grey #10 Hegerl 2007 beige #9 Juckes darj grey #8 Darrigo 2006green4 #7 Oerlemans 2005 brown #5 huang 2004 indigo #COLLATE year=ts(900:2000,start=900) spaghetti<-ts.union(year,CRU,MBH99,J98,CL2.Jns11,briffa01,huang,moberg,oerl,darrigo,juckes,hegerl,esper.fitted) spaghetti=ts(spaghetti[900:2000,],start=900) write.table(spaghetti,file="d:/climate/data/spaghetti/newscientist.dat",sep="\t") spaghetti=spaghetti[,2:ncol(spaghetti)] ##EXAMINE CENTERING # year CRU MBH99 J98 CL2.Jns11 briffa01 huang #1920.50000000 -0.23262500 -0.07874625 -0.17957250 0.05700000 -0.16013750 -0.14901791 # moberg oerl darrigo juckes hegerl esper.fitted # -0.15816250 -0.18600000 -0.19675000 0.01041728 0.06887500 -0.16025000 #MBH99 is 1902-1980; as is CL2.Jnns11 ##spline.response spline.response<-function(year) { g<-cos (2*pi/year) spline.response<-1/ ( ( 1-.5)/.5 * (g+2)/ (12*(g-1)^2 ) +1 ) spline.response } p=spline.response(7) spaghetti.smooth=spaghetti M=ncol(spaghetti) for (i in 1:M){ temp=!is.na(spaghetti[,i]) y<-year[temp] w<-spaghetti[temp,i] z<-smooth.spline(y,w,spar=p) ; spaghetti.smooth[temp,i]<-z$y } sd0=apply(spaghetti.smooth[(1856:1960)-899,],2,sd) sd0/sd0[1] # CRU MBH99 J98 CL2.Jns11 briffa01 huang moberg # 0.1723682 0.1486458 0.2144272 0.1041350 0.1013811 0.1746599 0.1046353 # oerl darrigo juckes hegerl esper.fitted # 0.2276833 0.2387728 0.1506280 0.2118783 0.1402965 m0=apply(spaghetti.smooth[(1856:1960)-899,],2,mean) # CRU MBH99 J98 CL2.Jns11 briffa01 huang moberg # -0.25 -0.27 -0.25 0.01 -0.19 -0.20 -0.19 # oerl darrigo juckes hegerl esper.fitted # -0.26 -0.27 -0.03 0.00 -0.21 col0=c("black","green2","red","yellow","blue","cyan","grey70","brown","blue3","grey80","pink","purple3") ts.plot(spaghetti.smooth,xlim=c(700,2000),lwd=2,col=col0,ylim=c(-1,.8)) legend(680,.7,fill=col0,legend=dimnames(spaghetti)[[2]]) [1] "CRU" "MBH99" "J98" "CL2.Jns11" "briffa01" [6] "huang" "moberg" "oerl" "darrigo" "juckes" [11] "hegerl" "esper.fitted" > write.table(spaghetti,file="d:/climate/data/eos/spaghetti.data.txt",quote=FALSE,sep="\t")