##REVIEW OF OSBORN AND BRIFFA [2006]
   #this provides statistics on bristlecones and the MJ04 PC1 used in the review   
   #this is script prepared quickly in the course of reviewing this article
   #at a couple of steps, it uses data that I've collated to make summaries
   # in each case, I've saved the summary, but show the extraction script here for reference	
 
###DIRECTORY AND FUNCTIONS
	mbhdata<-"ftp:/holocene.evsc.virginia.edu/pub/MBH98"
	#mbhdata<-"c:/climate/data/new" #I have copied the MBH98 directory to this location as I'm blocked from Mann's site
	source("http://www.climateaudit.org/scripts/functions.osborn06.txt")

####LOAD PC1 VERSION FROM JONES AND MANN 2004
    #this proxy is identified as series 1 in Table S1 in the O&B SI
    #download version archived at WDCP in connection with Jones and Mann [2004]
 	url<-"ftp://ftp.ngdc.noaa.gov/paleo/contributions_by_author/jones2004/jonesmannrogfig4a.txt"
	g<-read.table(url,skip=9,fill=TRUE,nrow=2000)
	dimnames(g)[[2]]<-c("year","wusa.mann","jacoby","wusa.briffa","jasper")
	g<-ts(g,start=g[1,1],end=g[nrow(g),1])
	temp<-(g== -99.99)
	g[temp]<-NA
	pc.mj04<-ts(g[200:1980,2],start=200) #recovers PC1 as time series

###IDENTIFY SITES IN PC1 STARTING IN 200
	url5<-file.path(mbhdata,"TREE/ITRDB/NOAMER/BACKTO_1000")
	loc<-file.path(url5,"noamer-itrdb-ad1000.txt")
	roster5<-read.table(loc);nrow(roster5)##27 series in AD1000 roster
	roster5<-strip(as.character(roster5[,2]),2) #strips off .txt suffix

	load("c:/climate/data/tree/northamerica.details.tab") 
		#this is my own collation of WDCP data and you'll have to bear with me on this
		#the short output from this calculation listed below is saved and is readable
	temp<-(details$start<201)&!is.na(details$start);sum(temp) #31 series identified
	test<-roster5[!is.na(match(roster5,details$id[temp]))];test #6 series isolated
		#"ca534" "ca535" "nm572" "nv515" "nv516" "ut509"  #Mann and Jones [2003] mention 6 series

	temp<-!is.na(match(details$id,test))
	test1<-details[temp,1:5];test1
	test1$jones<-jones(test1[,4],test1[,5]);test1$jones # 733 733 807 733 734 734
	write.table(details[temp,1:5],file="c:/climate/data/osborn06/bristle.site.txt",sep="\t")
		#this is saved at http://www.climateaudit.org/data/osborn06/bristle.site.txt

	 #       id                     location type   lat    long jones
	#476  ca534               SHEEP MOUNTAIN PILO 37.22 -118.13   733
	#477  ca535              METHUSELAH WALK PILO 37.26 -118.10   733
	#2335 nm572 EL MALPAIS NATIONAL MONUMENT PSME 34.58 -108.06   807
	#2390 nv515                INDIAN GARDEN PILO 39.05 -115.26   733
	#2391 nv516                   HILL 10842 PILO 38.56 -114.14   734
	#2934 ut509                MAMMOTH CREEK PILO 37.39 -112.40   734

	test1<-read.table("http://www.climateaudit.org/data/osborn06/bristle.site.txt",header=TRUE,sep="\t",nrow=6)

###COLLATE 6 SITES
     #collate site chronologies: I've saved R-tables for all WDCP chronologies and this is used here
     #this shows the collation; the collation itself is saved at climateaudit
	combine<-NULL
	for (i in 1:6) {
		load(file.path("c:/climate/data/tree/chronologies/r/northamerica",paste(test1$id[i],"crn.tab",sep=".")))
		M<-tsp(chron.crn)[1]
		chron.crn<-ts(chron.crn[(200-M+1):length(chron.crn)],start=200)
		combine<-ts.union(combine,chron.crn)
		}
	combine<-extend.persist(combine)
	combine<-cbind(tsp(combine)[1]:tsp(combine)[2],combine)	
	combine<-combine[1:(1980-199),]
	dimnames(combine)[[2]]<-c("year",as.character(test1$id))
	dim(combine) # 1781    7
	write.table(combine,file="c:/climate/data/osborn06/tree.collation.AD200.txt",sep="\t",row.names=FALSE)
		#this is saved at http://www.climateaudit.org/data/osborn06/tree.collation.AD200.txt

###CALCULATE PC1 USING MBH98 "METHOD"
 	options(digits=4)
	combine<-read.table("http://www.climateaudit.org/data/osborn06/tree.collation.AD200.txt",sep="\t",header=TRUE)
	combine<-ts(combine[,2:7],start=200)
	combine[1,]
		#ca534 ca535 nm572 nv515 nv516 ut509 
		# 1071  1048  1087  1106   791  1257 
	combine.mannomatic<-mannomatic(combine)
	pca.mannomatic<-svd(combine.mannomatic)
	round(pca.mannomatic$v[,1],4)
		## 0.8950 0.1185 0.1297 0.2063 0.2749 0.2235
	pca.mannomatic$v[,1]^2
		## 0.80105 0.01404 0.01682 0.04256 0.07559 0.04994
		#the Sheep Mountain site has over 80% of the eigenvector weighting

##REPLICATION OF ARCHIVED PC1 BY SPLICING
     #this procedure is not documented anywhere but seems to be what Mann did	
     #see comments at www.climateaudit.org 	
	use0<-"pairwise.complete.obs"
	pc1f<-read.table(file.path(mbhdata,"TREE/ITRDB/NOAMER/BACKTO_1000-FIXED/pc1-fixed.dat"))

	pcversions<-ts.union(pc.mj04,ts(pca.mannomatic$u[,1],start=200),ts(pc1f[,2],start=1000))
	dimnames(pcversions)[[2]]<-c("mj04","PC1","pc1f")
	cor(pcversions,use=use0)	
 	#      mj04    PC1   pc1f
	#mj04 1.0000 0.9430 0.8852
	#PC1  0.9430 1.0000 0.7885
	#pc1f 0.8852 0.7885 1.0000

  #show that mj04 is rescaling of mannomatic PC1 up to 1700
	cor(pcversions[1:(1699-199),],use=use0)
   	#       mj04       PC1      pc1f
 	#         mj04       PC1      pc1f
	#mj04 1.0000000 0.9998378 0.8233237
	#PC1  0.9998378 1.0000000 0.8217405
	#pc1f 0.8233237 0.8217405 1.0000000

  #show that mj04 is re-scaling of pc1-fixed from 1700 to 1980
	cor(pcversions[(1700-199):(1980-199),],use=use0)
	#          mj04       PC1      pc1f
	#mj04 1.0000000 0.8087152 0.9999243
	#PC1  0.8087152 1.0000000 0.8091937
	#pc1f 0.9999243 0.8091937 1.0000000


##DO COVARIANCE AND CORRELATION PCs FOR AD200 NETWORK
	pca.covariance<-prcomp(combine)
	pca.correlation<-prcomp(combine,scale=TRUE)

	pcversions<-ts.union(pcversions,ts(cbind(pca.covariance$x[,1],pca.correlation$x[,1]),start=200))
	dimnames(pcversions)[[2]]<-c("mj04","PC200","pc1f","pc.cov","pc.cor")
	cor(pcversions,use=use0)
	#          mj04   PC200    pc1f  pc.cov  pc.cor
	#mj04    1.0000  0.9430  0.8852 -0.7002 -0.8235
	#PC200   0.9430  1.0000  0.7885 -0.7585 -0.8750
	#pc1f    0.8852  0.7885  1.0000 -0.4811 -0.5960
	#pc.cov -0.7002 -0.7585 -0.4811  1.0000  0.9361
	#pc.cor -0.8235 -0.8750 -0.5960  0.9361  1.0000
	
   #calculate "hockeystat"
	pcversions[,4:5]<- -pcversions[,4:5] #flip upside-down hockey sticks; sign of PC is unoriented
	f<-function(test) {N<-length(test)-78; f<-  ( mean(test[N:(N+79)],na.rm=TRUE)- mean(test[1:(N+79)],na.rm=TRUE) )/sd(test[1:(N+79)],na.rm=TRUE); f }
	apply(pcversions[801:(1980-199),],2,f)
	# mj04  PC200   pc1f pc.cov pc.cor 
	#1.0681 1.5903 0.9689 0.6251 0.8206 

##PLOT DIFFERENT VERSIONS
	nf<-layout(array(1:3,dim=c(3,1)),heights=c(1.1,1,1.4))
	par(mar=c(0,3,1,1))
	plot(200:1980,pcversions[,1],type="l",axes=FALSE,xlab="",ylab="")
	axis(side=2,las=1);axis(side=2,labels=FALSE);box()
	text(250,3,pos=4,font=2,"Jones and Mann PC1")
	par(mar=c(0,3,0,1))
	plot(200:1980,scale(pcversions[,4]),type="l",axes=FALSE,xlab="",ylab="")
	axis(side=2,las=1);axis(side=2,labels=FALSE);box()
	text(250,2.7,pos=4,font=2,"Covariance PC1")
	par(mar=c(3,3,0,1))
	plot(200:1980,scale(pcversions[,5]),type="l",axes=FALSE,xlab="",ylab="")
	axis(side=2,las=1);axis(side=2);box()
	text(250,2.7,pos=4,font=2,"Correlation PC1")


##MAKE TEMPERATURE COLLATIONS
    #this uses collation of HadCRU2 into R table which was made using script posted up at climateaudit.org
    #the collation is shown here, but the output is separately saved
	load("c:/climate/data/jones/hadcrut2.tab") #object v.tab is monthly starting (1856,1); nrow 1791
	vbristle<-v[,c(733,734,807)] #these are 3 cells used in MJ04 AD200 network
	rm(v)

	annualb<-apply(vbristle[1:1788,],2,annavg)
	annualb<-cbind(1856:2004,annualb)
	annualb[is.na(annualb)]<-NA
	annualb<-data.frame(annualb);names(annualb)<-c("year","733","734","807")

	f<-function(x) { f<-sum(x*weights0,na.rm=TRUE)/sum(weights0[!is.na(x)]);f}
	weights0<-c(1/3,1/3,1/3)
	annualb$unweighted<-apply(annualb[,2:4],1,f)

	weights0<-tapply(pca.mannomatic$v[,1],test1$jones,sum)
	annualb$weighted1<-apply(annualb[,2:4],1,f)
		#this is weighting in same proportion as PC1

	weights0<-tapply(pca.mannomatic$v[,1]^2,test1$jones,sum)
	annualb$weighted2<-apply(annualb[,2:4],1,f)
	annualb[,2:7]<-round(annualb[,2:7],1)
	write.table(annualb,file="c:/climate/data/osborn06/bristle.gridcells.txt",sep="\t")
	
	round(cor(annualb,use=use0),2)
	#            X733 X734 X807 unweighted weighted1 weighted2
	#X733       1.00 0.70 0.54       0.86      0.97      1.00
	#X734       0.70 1.00 0.75       0.92      0.86      0.79
	#X807       0.54 0.75 1.00       0.90      0.76      0.69
	#unweighted 0.86 0.92 0.90       1.00      0.96      0.93
	#weighted1  0.97 0.86 0.76       0.96      1.00      0.99
	#weighted2  1.00 0.79 0.69       0.93      0.99      1.00

#CALCULATE ANNUAL CORRELATIONS TO GRIDCELL
	annualb<-read.table("http://www.climateaudit.org/data/osborn06/bristle.gridcells.txt",sep="\t",header=TRUE)
	annualb<-ts(annualb[,2:7],start=1856)
	loc1<-c(1,1,3,1,2,2) #manually look up gridcell for each series

     #calculate correlations for all 6 tree ring sites from 1870-1980; all cells available from 1870 on
	stat<-rep(NA,6)
	for (i in 1:6) {
		z<-ts.union(combine[,i],annualb[,loc1[i] ])
		stat[i]<-cor(z[(1870-199):(1980-199),],use=use0)[1,2]
		}
	cor.bristle<-stat;cor.bristle
		#0.02861566 -0.14388843 -0.27971763 -0.12918283 -0.35482184 -0.27280318
		#only one site has positive correlation; 5 have negative correlation to temperature

    #calculate correlations for mannomatic PC and "fixed" mannomatic PC to weighted gridcell temperature
	cor(annualb[(1870-1855):(1980-1855),5],pcversions[(1870-199):(1980-199),2])
		#   -0.1223255
	cor(annualb[(1870-1855):(1980-1855),5],pcversions[(1870-199):(1980-199),"pc1f"])
		#   0.0365

#DECADALLY SMOOTHED
	annual.smooth<-ts(annualb[,1:3],start=1856)
	f<-function(x){f<-filter.combine.pad(x,gaussian.filter.weights(13,6))[,2];f}
	for (i in 1:3) {annual.smooth[,i]<-f(annual.smooth[,i])}

	pc.smooth<-pcversions
	for (i in 1:ncol(pcversions)){ pc.smooth[,i]<-f(pc.smooth[,i])}

	bristle.smooth<- combine
	for (i in 1:ncol(bristle.smooth)) {bristle.smooth[,i]<-f(bristle.smooth[,i])}

	stat<-rep(NA,6)
	for (i in 1:6) {
		z<-ts.union(bristle.smooth[,i],annual.smooth[,loc1[i] ])
		stat[i]<-cor(z[(1870-199):(1980-199),],use=use0)[1,2]
		}
	cor.smooth<-stat;cor.smooth
		#0.05942632 -0.15839882 -0.18256974 -0.11032780 -0.32755484 -0.22864735
		#only one site slightly positive in decadally smoothed correlation

    #check against weighted average gridcell temperature
	f<-function(x) { f<-sum(x*weights0,na.rm=TRUE)/sum(weights0[!is.na(x)]);f}
	weights0<-tapply(pca.mannomatic$v[,1],test1$jones,sum)
	annual.smooth<-data.frame(annual.smooth)
	annual.smooth$weighted1<-apply(annual.smooth[,1:3],1,f)
	
	cor(pc.smooth[(1856-199):(1980-199),2],annual.smooth$weighted1[1:(1980-1855)],use=use0)
		# -0.3540
	cor(pc.smooth[(1856-199):(1980-199),"pc1f"],annual.smooth$weighted1[1:(1980-1855)],use=use0)
		# -0.2628
	

	plot(1856:1980,scale(test.smooth[(1856-199):(1980-199),"pc1f"]),type="l",xlab="",ylab="")
	lines(1856:1980,scale(annual.smooth$weighted1[1:(1980-1855)]),col="red")

#PLOT PC AGAINST GRIDCELL
    nf<-layout(array(1:2,dim=c(2,1)),heights=c(1.1,1.4))
	par(mar=c(0,3,1,1))
	plot(1870:1980,scale(annualb[(1870-1855):(1980-1855),5]),type="l",xlab="",ylab="",axes=FALSE)
	lines(1870:1980,scale(pcversions[(1870-199):(1980-199),"pc1f"]),col="red")
	axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0)
	text(1870,2.5,pos=4,font=2,"Annual")

	par(mar=c(3,3,0,1))
	plot(1870:1980,scale(annual.smooth$weighted1[(1870-1855):(1980-1855)]),type="l",xlab="",ylab="",axes=FALSE,ylim=c(-3,3))
	lines(1870:1980,scale(pc.smooth[(1870-199):(1980-199),"pc1f"]),col="red")
	axis(side=1);axis(side=2,las=1);box();abline(h=0)
	text(1870,2.5,pos=4,font=2,"Decadal Smooth")