##FUNCTIONS
###############
library(XLConnect)
library("R.matlab")
library(reshape2)

bin=function(A,scale=F,center=T){
  work1= with(A, tapply(sst,list(period,id) ,mean))
  work1=scale(work1,center=center,scale=scale)
  work=melt(work1,id.vars=c("period","id","sst"))
  names(work)=c("period","id","sst")
  work$ocean=info$ocean[match(work$id,info$id)]
  work=work[!is.na(work$sst),]
  work$id=factor(work$id)	
  work$ocean=factor(work$ocean)
  return(work)
}

ext=function(x) c(x,x[length(x)])

#TARGET RECONSTRUCTIONS OF FIGURE 2A
#######################################333
loc="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/pages2k/Ocean2kLR2015.zip"
download.file(loc,"d:/temp/temp.zip",mode="wb")
handle=unz("d:/temp/temp.zip","Ocean2kLR2015/composites_shipped/compositework.mat","rb")
cwork=readMat(handle)
close(handle)


#INFO SUMMARY
###############
loc="http://www.climateaudit.info/data/multiproxy/ocean2k/info_ocean2k.csv"
info=read.csv(loc)

w=c(arc=.034,atl=.183,ind=.113,med=.008, pac=.384,sou=.278)
  #weights

##PROXY DATA
###############
#setwd("d:/climate/data/multiproxy/ocean2k")
#wb=loadWorkbook("Ocean2kLR2015sst.xlsx")
loc="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/pages2k/Ocean2kLR2015sst.xlsx"
dest="d:/temp/temp.dat"
download.file(loc,dest,mode="wb")
wb=loadWorkbook(dest)

K=57
O=NULL
for(i in 1:57) {
work=readWorksheet(wb,sheet=4,startRow=1,startCol=2*i-1,endCol=2*i,colTypes=rep("numeric",2))
 names(work)=c("year","sst")
 work=work[!is.na(work$sst),]
  work$id=info$id[i]
 O=rbind(O,work)
}
O$ocean=info$ocean[match(O$id,info$id)]
O=O[O$year>0&O$year<=2000,]

M=200
O$period= factor( M*floor((O$year-1)/M)+M/2)
O200=O
Bin=Bin200= bin(O200,center=T,scale=F)
Bin200scale= bin(O200,center=T,scale=T)

P=Bin200scale
P=P[!is.na(P$sst),]
P$id=factor(P$id) #two NA-ed 
out=with(P,tapply(sst,list(period,ocean),mean,na.rm=T))
emulation=apply(out,1,function(x) weighted.mean(x,w,na.rm=T ))
X=cbind(emu=emulation,archive=rev(cwork$wavemn))
cor(X) #0.9999759
range(X[,1]-X[,2])
  #[1] -0.003437707  0.005643450