############################################################
#3. AVERAGE TEMPERATURE TIME SERIES REPLICATION
####################################
# the new SI contains a listing of the grid cell standard deviations used in MBH98
# these calculations were emulated and compared to new SI values

	mbhscripts<-"http://www.climate2003.com/toolbox/scripts"
	source(file.path(mbhscripts,"MBH98config3.txt"))
	source(file.path(mbhscripts,"treeconfig.functions.txt"))

###########
#data sets for sparse and dense at NOAA
# summary data for NH temperature at NOAA and the new SI

####DENSE

  ##LOAD NOAA VERSIONS
	#no information at NOAA on zero
	#only from 1902-1998
	url<-"ftp://ftp.ngdc.noaa.gov/paleo/paleocean/by_contributor/mann1998/nhem-dense.dat"
	g<-read.table(url,skip=1,nrow=96)
	g<-rbind(g,c(1998,0,0.85)) #manual to pick up *
	g<-ts(g[,2:3],start=1902,end=1998)
	dimnames(g)[[2]]<-c("recon","instr")
	dense<-g
	tsp(dense) #1854 1998  

  	#this is the same as http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/nhmean.txt from 1902-1995
	#Nature just goes to 1995
	#this is same as col 2 of 
	url<-"ftp://ftp.ngdc.noaa.gov/paleo/paleocean/by_contributor/mann1998/mannnhem.dat"
	g<-read.table(url,skip=1)
	temp<-(g[,3]==0)
	g[temp,3]<-NA
	NH.NOAA<-ts(g[!temp,3],start=g[!temp,1][1],end=g[!temp,1][length(g[!temp,1])])

	apply(dense[1:79,],2,mean) # 3.797468e-09 -1.265823e-09 
	#the dense data is zeroed to 1902-1980
	apply(dense[1:92,],2,mean) #3.260870e-09 4.342307e-02  #values of "recon" are zero from 1981-1993
	#dense is zeroed on 1902:1980

	#except in 1998

 ##LOAD INSTRUMENTAL DATA (UNREGULARIZED)
	#what is anomaly: Jones and Briffa 1992. 
	load(file.path(dirname(mbhdata),"MBH98.2004/v.combined.tab"))
	dim(v) # 1680 2592
	temp<-!is.na(match(c(1:2592),tempmann[nhtemp])) #172 cells #NH cells are in first 1296 cells
	v<-v[,temp] #1680 707
	v<-v[((1901-1853)*12+1):nrow(v),]/100
	dim(v) # 1104  707

  #CHECK ANOMALY PERIOD FOR 1961-1990
	test<-v[,200]
	test<-array(test,dim=c(12,1104/12))
	test<-t(test)
	apply(test[(1961-1901):(1990-1901),],2,mean) # [1] -0.12033333 -0.17066667  0.05800000  0.04566667 -0.06200000 -0.03400000 -0.18666667 -0.09866667  0.03633333
	#[10]  0.00900000 -0.21333333 -0.29333333
	apply(test[(1902-1901):(1980-1901),],2,mean) # [1] -0.297974684 -0.477974684 -0.001139241 -0.166075949 -0.041772152 -0.157215190  0.066075949 -0.017974684
 	#[9] -0.034050633  0.058860759 -0.297088608 -0.453417722
	apply(test[(1950-1901):(1979-1901),],2,mean)
	# [1]  0.017000000 -0.042333333 -0.087666667  0.006333333  0.008333333  0.001333333  0.064000000  0.074000000
 	#[9]  0.074333333 -0.036666667 -0.025000000  0.009333333
	apply(test[(1902-1901):(1993-1901),],2,mean)
	#[1] -0.236086957 -0.519347826 -0.033152174 -0.169782609 -0.115217391 -0.134239130 -0.001304348  0.006847826
 	#[9] -0.010434783  0.094565217 -0.397826087 -0.477717391


  #CALCULATE 2592 COS WEIGHTS FOR JONES INDEX SERIES
	theta<-87.5
	mu<-rep(theta,72)
	for (k in 1:35) {
		theta<-theta-5
		mu<-c(mu,rep(theta,72))
		}	
	mu<-pi*mu/180
	mu<-cos(mu) 
	mu<-mu[temp]

   #ANNUAL AVERAGE
	areal.avg<-function(x,weight=mu) {areal.avg<-sum(weight*x,na.rm=TRUE)/sum(weight*!is.na(x));areal.avg}
	T<-apply(v,1,areal.avg)
	annual.SI<-ts(annavg(T),start=1902,end=1993)
	annual.adj<-annual.SI -mean(annual.SI[1:79]) #center to 1902-1980 zero

   #COMPARE DENSE TO NHMANN
	combine<-ts.union(dense,annual.SI)
	cor(combine,use=use0) #   0.9696518 
	mean(dense[,"instr"]-annual.SI) #  0.1013694
	mean(annual.SI[(1902-1901):(1980-1901)])  # -0.1041084
	annual.adj<-annual.SI-  mean(annual.SI[(1902-1901):(1980-1901)])  #zero to 1902-1980
	
par(mfrow=c(1,2))
	ts.plot(dense[,"instr"],xlab="",ylab="deg C") #dense is centered on 1902-1980
	lines(1902:1993,annual.adj,col="red")
  
	plot(c(dense[1:92,"instr"]),annual.adj,xlab="Archived at NOAA",ylab="Calculated")
	abline(0,1)

############
##SPARSE CALCULATIONS
#############
	mbhscripts<-"http://www.climate2003.com/toolbox/scripts"
	source(file.path(mbhscripts,"MBH98config3.txt"))
	source(file.path(mbhscripts,"treeconfig.functions.txt"))

###########
#data sets for sparse and dense at NOAA
# summary data for NH temperature at NOAA and the new SI

  ##LOAD NOAA VERSIONS
	#no information at NOAA on zero
	#dense is zero on 1902-1980
	url<-"ftp://ftp.ngdc.noaa.gov/paleo/paleocean/by_contributor/mann1998/nhem-sparse.dat"
	g<-read.table(url,skip=1)
	g<-ts(g[,2:3],start=g[1,1],end=g[nrow(g),1])
	dimnames(g)[[2]]<-c("instr","recon")
	sparse<-g
	tsp(sparse) #1854 1993
	round(apply(sparse[(1902-1853):(1980-1853),],2,mean),3)  #-0.003  0.000 
	apply(sparse[(1902-1853):(1993-1853),],2,mean)  #-2.572732e-03  1.009812e-08 


 ###LOAD SPARSE MASK
	load(file.path(mbhdata,"mask.MBH.tab"))
	length(mask.MBH) #219

 ##LOAD INSTRUMENTAL DATA (UNREGULARIZED)
	#what is anomaly: Jones and Briffa 1992. 
	load(file.path(dirname(mbhdata),"MBH98.2004/v.combined.tab"))
	dim(v) # 1680 2592
	temp<-!is.na(match(c(1:2592),mask.MBH)) & c(rep(TRUE,1296),rep(FALSE,1296)) #172 cells #NH cells are in first 1296 cells
	v<-v[,temp] #1680 172

  #CALCULATE 2592 COS WEIGHTS FOR JONES INDEX SERIES
	theta<-87.5
	mu<-rep(theta,72)
	for (k in 1:35) {
		theta<-theta-5
		mu<-c(mu,rep(theta,72))
		}	
	mu<-pi*mu/180
	mu<-cos(mu) 
	mu<-mu[temp]

   #ANNUAL AVERAGE
	areal.avg<-function(x,weight=mu) {areal.avg<-sum(weight*x,na.rm=TRUE)/sum(weight*!is.na(x));areal.avg}
	T<-apply(v,1,areal.avg)/100
	annual.SI<-ts(annavg(T),start=1854,end=1993)
	cor(sparse[,"instr"],annual.SI) #  0.9989637
	mean(sparse[,"instr"]-annual.SI) #  0.08828018
	annual.adj<-annual.SI-mean(annual.SI[(1902-1853):(1980-1853)]) #zero to 1902-1980

  #OTHER DIRECTION #checked to see that values not effected if annual average taken first before areal average
  #PLOTs	

par(mfrow=c(1,1))
	ts.plot(sparse[,"instr"],xlab="",ylab="deg C")
	lines(1854:1993,annual.adj,col="red")

	plot(c(sparse[,"instr"]),annual.adj,xlab="Archived at NOAA",ylab="Calculated")
	abline(0,1)
	

#COMPARISONS
#Jones 2000 1856-2000
url<-"http://cdiac.esd.ornl.gov/ftp/trends/temp/jonescru/nh.dat"
readLines(url)
g<-read.table(url,skip=16,nrow=145,header=TRUE)
NH.2000<-ts(g[,14],start=g[1,1],end=g[nrow(g),1])

	ts.plot( ts.union(NH.2000,NH.NOAA),col=col1,xlim=c(1854,2000))
	#these are highly correlated but there is a difference in mean
	combine<-ts.union(NH.2000,-NH.NOAA)
	plot.ts(apply(combine,1,sum,na.rm=FALSE))
 	mean(apply(combine,1,sum,na.rm=FALSE),na.rm=TRUE) #[1] -0.1138652

	combine<-ts.union(NH.2000,NH.NOAA)
	plot.ts(combine)