#COLLATION FOR OSBORN aND BRIFFA 2006


#RCS CHRONOLOGY FUNCTION
	#returns time series (chronology), dataframe of 5-year values, delta for each measurement , smooth for each mea\surement, 
	# parameters (complete fit for RCS.spline method), method
		#mthods: linear -= uses straight line fit
		#         "mod.neg.exp" - fits generalzied exponential using start parameters

     RCS.chronology<-function (tree,method,parameters=NULL,interval=5,units=1) {
	dimnames(tree)[[2]][4]<-"x"
	tree<-tree[!is.na(tree$x),]
	tree$agegroup<-median(1:interval)+interval*trunc( (tree$age-1)/interval)
	tree$id<-factor(tree$id)
	tree$agegroupf<-factor(tree$agegroup)
	#year5f<-factor(tree$group)
	year5avg <-tapply(tree$x,tree$agegroupf,mean,na.rm=TRUE)/units
	agegroup<-as.numeric(names(year5avg))
	N<-length(agegroup)

	z<-cbind(agegroup,year5avg)	
	z<-data.frame(z)
	if (method=="mod.neg.exp") {fm <- nls( year5avg ~ A+B*exp(-C*agegroup),data = z,
		start = list( A = mean(z$year5avg), B=10, C= .01 ),
		alg = "default", trace = TRUE,control=nls.control(maxiter=200, tol=1e-05, minFactor=1e-10)); 
		B<-coef(fm); 
		fitted<- B[1]+B[2]*exp(-B[3]*(1:agegroup[N]))  }
	if (method=="nls") {fm <- nls(x ~ A+B*exp(-C*age),data = tree,
		start = list( A = mean(tree$x), B=10, C= .01 ),
		alg = "default", trace = TRUE,control=nls.control(maxiter=200, tol=1e-05, minFactor=1e-10)); 
		B<-coef(fm); 
		fitted<- B[1]+B[2]*exp(-B[3]*(1:max(tree$age)))  }
	if (method=="linear") {fm<-lm(year5avg~agegroup,data=z);
		B<-coef(fm); fitted<-B[1]+B[2]*(1:agegroup[N])}
	if (method=="briffa") { B<-parameters;
 		if (length(B) ==2) fitted<- B[1]+B[2]*(1:agegroup[N]);
		if (length(B) ==3) fitted<- B[1]+B[2]*exp(-B[3]*(1:agegroup[N])) }
	if (method=="RCS.spline") { p<-spline.response(parameters[1]*N);
		fit.model<- smooth.spline( agegroup, year5avg, spar = (1-p));
		fitted<-predict(fit.model,(1:agegroup[N]))$y;
		B<-fitted}
	fitted<-ts( fitted, start=1, end=agegroup[N])
	tree$delta<-NA
	tree$smooth<-NA
	for (k in 1:length(levels(tree$id))) {
		temp<-(tree$id==levels(tree$id)[k])
		y<-tree$year[temp]
		w<-tree$x[temp]
		n<-length(w);
		tree$smooth[temp]<-fitted[tree$age[temp]] 
		tree$delta[temp]<-w / tree$smooth[temp]
		}
	yearf<-factor(tree$year)
	series<-tapply(tree$delta,yearf,mean,na.rm=TRUE)
	series<-ts(series,start=min(tree$year),end=max(tree$year))	
	RCS.chronology<-list(series,z,tree$delta,tree$smooth,B,method)
	names(RCS.chronology)<-c("series","hist.curve","delta","smooth","coefficients","method")
	RCS.chronology
	}


##COLLATE SOURCES
     #14 series
     #from 800 to 1995 #14 series #make array

	osborn<-array(NA,dim=c((1995-799),14))
	id<-c("MBHPC1","Alberta","Foxtail","Quebec","Chesapeake","WGreen","Neth.doc","Tirol","Tornetrask",
		"Yamal","Mangazeja","Taimyr","Mongolia","Chinayang")
	dimnames(osborn)[[2]]<-id

#1	MBH PC1 # 6 sites all available. #PC1 at MJ04
#2	Alberta: Jasper/Icefields. Jasper crn at Briffa. No Luckman rwl archive. No archive of newer Icefields version (Wilson)
#3	Boreal/Upperwright foxtails: questions as to provenance - Graumlich has lost rwl data c 1999. Maybe this is other PIBA??
#4	Quebec - no archive for Payette version. OB refer to cana169.rwl which is available (but goes to 1989)
#5	Chesapeake - OK. Both JM04 and WDCP
#6	Fisher's Greenland O18 stack - OK. Both JM04, WDCP (and MBH)
#7	Netherlands documentary. OK version at JM04 a little different than at van Engeln
#8	Austria (TIrol). Esper version NA; rwl version cited here germ21 goes only to 1821 and seems inconsistenbt
#9	Tornetrask. crn at Briffa. Old rwl (probably included) is at WDCP, new rwl NA
#10	Yamal. crn at Briffa. Hantemirov sent me rwl
#11	Mangazeja. No archive of crn. rwl can be compiled, but reconciliation needed
#12	Taimyr. crn at Briffa. No rwl.
#13	Mongolia. a crn version at JM04. partial rwl at WDCP mong003
#14	Yang's China composite - this is compromised by subset issues.

###############
##A. FROM JONES AND JONES 2004
   #sites 1,2,5,6,7, 13
	#url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/mann2003b/mann2003b.txt"
	#All series are anomalies based on 1961-1990 instrumental reference period. 
	#Data are presented as both decadally-resolved series (added to this file 12/2004)
	#and 40 year smoothed versions of the decadally-resolved reconstrucitons 
	#(as shown in Figure 2 of Mann and Jones).

   ##NORTH AMERICA - 2
	url<-  "ftp://ftp.ngdc.noaa.gov/paleo/contributions_by_author/jones2004"
	loc<-file.path(url,"jonesmannrogfig4a.txt")
	g<-read.table(loc,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
	tsp(g)#1 2000
	apply(g[(1750:1950),2:5],2,mean,na.rm=TRUE) #all zeroed on 1750-1950
		#    wusa.mann        jacoby   wusa.briffa        jasper 
		#-3.866448e-18  1.840796e-03 -5.970149e-04  4.029851e-03 
	osborn[1:(1995-799),c(1,2)]<-g[800:1995,c(2,5)]

   ##NORTH ATLANTIC
	url<-"ftp://ftp.ngdc.noaa.gov/paleo/contributions_by_author/jones2004/jonesmannrogfig4b.txt"
	g<-read.table(url,skip=8,fill=TRUE,nrow=2000)
	dimnames(g)[[2]]<-c("year","wgreen.fisher","chesapeake","vinther")
	g<-ts(g,start=g[1,1],end=g[nrow(g),1])
	temp<-(g== -99.99)
	g[temp]<-NA
	osborn[1:(1995-799),c(6,5)]<-g[800:1995,c(2,3)]

    #EUROPE
        #for Tornetrask, Briffa version is used
	url<-"ftp://ftp.ngdc.noaa.gov/paleo/contributions_by_author/jones2004/jonesmannrogfig4c.txt"
	g<-read.table(url,skip=8,fill=TRUE,nrow=2000)
	name1<-c("tornetrask","polarurals","vanengeln")
	dimnames(g)[[2]]<-c("year",name1)
	g<-ts(g,start=g[1,1],end=g[nrow(g),1])
	temp<-(g== -99.99)
	g[temp]<-NA
	osborn[1:(1995-799),7]<-g[800:1995,4]

    ###EAST ASIA
	url<-"ftp://ftp.ngdc.noaa.gov/paleo/contributions_by_author/jones2004/jonesmannrogfig4d.txt"
	g<-read.table(url,skip=6,fill=TRUE,nrow=2000)
	name1<-c("mongolia")
	dimnames(g)[[2]]<-c("year",name1)
	g<-ts(g,start=g[1,1],end=g[nrow(g),1])
	temp<-(g== -99.99)
	g[temp]<-NA
	osborn[1:(1995-799),13]<-g[800:1995,2]

 
##B. WDCP  SOURCE DATA FOR YANG 2002
   #series 14
   #not used in JM04; no MJ03 archive
   #this is only in decadal version in original; many components are smoothed even more (50 years)
  	url<-"ftp://ftp.ngdc.noaa.gov/paleo/contributions_by_author/yang2002/china_temp.txt"
	g<-read.table(url,skip=74,fill=TRUE) #goes from 0 to 1990
	dimnames(g)[[2]]<-c("Year","Complete","H-res","Sample Number","Weighted")
	f<-approxfun(g[,1],g[,2])
	h<-ts(f(1:1990),start=1,end=1990)
	osborn[1:(1990-799),14]<-h[800:1990]

##C. FROM BRIFFA WEBSITE
    #series 9,10,12 
      #Tornetrask e1r-N-14       Yamal -f1r T-20  Taimyr  gr-Z -26   
	url<-"http://www.cru.uea.ac.uk/cru/people/briffa/qsr1999/qsrfig1.csv"
	g<-read.csv(url,fill=TRUE,skip=11)#dim 1999x30
	new<-g[,c(1,  c(14,20,26))]
 	dim(new)#[1] 2000    4
 	new[1,1] #[1] 1
	osborn[1:(1995-799),c(9,10,12)]<-as.matrix(new[800:1995,2:4])

##D.  ESPER - NO ARCHIVED VERSIONS

  #4. Quebec - Briffa says Esper (S1) then S12 - cana169 (presumable cana169w)
 		#www.sciencemag.org/cgi/content/full/295/5563/2250/DC1

	load("c:/climate/data/tree/chronologies/r/northamerica/cana169w_crns.crn.tab")
	tsp(chron.crn) # 1352 1989

	load("c:/climate/data/tree/measurements/r/northamerica/cana169w.rwl.tab")
	range(tree$year) # 1352 1989
		#Black Spruce located at +50  5520N -7750 (W)
	RCS.cana169<-RCS.chronology(tree,"mod.neg.exp") #fails after 200 iterations
	RCS.cana169<-RCS.chronology(tree,"nls") #converges
	RCS.cana169$coefficients
		#            A             B             C 
		#-54.000520938  94.363071679   0.001192954 
	tsp(RCS.cana169$series) # 1352 1989   
	osborn[(1352-799):(1947-799),4]<-RCS.cana169$series[1:(1947-1351)]
		#this is inexplicably terminated at 1947: Briffa SI
		
		#however Esperr cited data from S. Payette and L. Filion (Quebec) provided Quebec data. 
		#and Esper graphic ends about 1820 as in PAyette articles
		#did not supply data to Esper directly and not archived

  #8. Tirol

	#Esper credits V. Siebenlist-Kerner ( Tirol) 
	#Osborn-Briffa cite germ21 by Eckstein from 46N,11E
	#the end date is inconsistent

	load("c:/climate/data/tree/measurements/r/europe/germ21.rwl.tab")
	range(tree$year) # 1466 1837
	RCS.germ21<-RCS.chronology(tree,"nls") #converges
 	RCS.germ21$coefficients
 		#          A            B            C 
 		#29.75575367 108.77715724   0.01487110 
	tsp(RCS.germ21$series) # 1466 1837
	osborn[(1466-799):(1837-799),8]<-RCS.germ21$series
	
  #11. Mangazeja

	#Esper credits S. Shiyatov (Polar Urals, Mangazeja),
	#OB say that it is russ067, russ068
	load("c:/climate/data/tree/measurements/r/asia/russ068w.rwl.tab")
	range(tree$year)
	tree0<-tree

      #there is no russ067w; therefore go to Schweingruber archive at WDCP

     #function to read from SChweingruber archive
      	WSL<-function(id) {	
		loc<-file.path("ftp://ftp.ngdc.noaa.gov/paleo/treering/updates/wsl/raw-data/singles",id)
		test<- class(try(readLines(loc)))
		if(test=="character") fred<-readLines(loc) else {
			loc<-file.path(url,id);
			fred<-readLines(loc) }
		tree.WSL<-data.frame( substr(fred,1,8),as.numeric(substr(fred,13,16)),as.numeric(substr(fred,29,32)),as.numeric(substr(fred,53,56)))
		temp<-(tree.WSL[,2]==9999)
		tree.WSL<-tree.WSL[!temp,]
		temp<- (tree.WSL[,3] < 0);
		tree.WSL[temp,3]<-NA
		temp<- (tree.WSL[,4] < 0);
		tree.WSL[temp,4]<-NA
		names(tree.WSL)<-c("id","year","mxd","rw")
		tree<-tree.WSL
		tree$age<-NA
		tree$start<-tree$id
		levels(tree$start)<-tapply(tree$year,tree$id,min)
		tree$start<-as.numeric(levels(tree$start))[tree$start]
		tree$age<-tree$year-tree$start+1
		WSL<-tree[,c(1:2,5,3:4)]
		WSL
		}

      #load Mangazeja data
	tree.mangazpc<-WSL("mangazpc.her")
	dim(tree.mangazpc) #9371 5
	range(tree.mangazpc$year)#[1] 1246 1969
	tree.mangazla<-WSL("mangazla.hbr")
	dim(tree.mangazla) #5673 5
	range(tree.mangazla$year)#[1] 1307 1698
	tree<-rbind(tree.mangazla,tree.mangazpc)
	RCS.mangaz<-RCS.chronology(tree,"nls") 
		#doesn't converge
	RCS.mangaz<-RCS.chronology(tree,"mod.neg.exp",parameters=c(25,50,.01))
		#doesn't converge
	RCS.mangaz<-RCS.chronology(tree,"RCS.spline",parameters=c(.125,5)) #experiment with curve
		#does converge - not a mod neg exp series	
	
	tsp(RCS.mangaz$series) #1246 1969
	osborn[(1246-799):(1969-799),11]<-RCS.mangaz$series


    #3. BOREAL/UPPERWRIGHT FOXTAILS
      #illustration in article starts in 800s, not at beginning
      #this is a real puzzle: 
	#not archived at WDCP
	#sites not discussed in cited article Lloyd and Graumlich
	#sites are discussed in Bunn et al [2005]
	#inquired for data from Andrew Bunn on Feb 8, 2006
	#Bunn said that he had never seen rwl files and, to his knowledge, they had been lost
		#when Graumlich moved to Montana (which seems to be 1999)		

	#as a substitute, I tried to use a couple of other archived PIBA sites

	#cirque Peak, Flower Lake, Timber Gap Upper, Timber Gp Lower
	#470          ca528                            FLOWER LAKE PIBA 36.46 -118.22     3291   898 1987  DONALD A. GRAYBILL
	#471          ca529                       TIMBER GAP UPPER PIBA 36.27 -118.36     3261   699 1987  DONALD A. GRAYBILL
	#472          ca530                            CIRQUE PEAK PIBA 36.27 -118.13     3505   917 1987  DONALD A. GRAYBILL
	#473          ca531                           ONION VALLEY PIBA 36.46 -118.21     2865  1027 1987  DONALD A. GRAYBILL
	#474          ca532                       TIMBER GAP LOWER PIBA 36.27 -118.37     3017  1050 1987  DONALD A. GRAYBILL
	#497   

	#north are ca528; ca531.  south are ca529,ca531, ca532
	#try Cirque Peak for Boreal Plateau and Flower Lake for Upperwright

	load("c:/climate/data/tree/measurements/r/northamerica/ca530.rwl.tab")
	range(tree$year) #917 1987
	RCS.ca530<-RCS.chronology(tree,"nls") #converges
 	RCS.ca530$coefficients
		#      A        B        C 
		#42.31875 22.96911  0.01654 

	load("c:/climate/data/tree/measurements/r/northamerica/ca528.rwl.tab")
	range(tree$year) #917 1987
	RCS.ca528<-RCS.chronology(tree,"nls") #converges
 	RCS.ca528$coefficients
		#       A        B        C 
		# 34.39954 13.67097  0.00269 

	z<-ts.union(RCS.ca528$series,RCS.ca530$series)
	dim(z)## 1090    2
	h<-apply(z,1,mean,na.rm=TRUE)  #898 1987
	
     #work into osborn
	osborn[(898-799):(1987-799),3]<-h


###SAVE OSBORN COLLATION
	osborn<-cbind(800:1995,osborn)
	dimnames(osborn)[[2]]<-c("year",id)
	save(osborn,file="c:/climate/data/osborn06/osborn.tab")
	write.table(osborn,file="c:/climate/data/osborn06/osborn.txt",sep="\t",row.names=FALSE)
	

##COMPARE GRAPHS