##KAPLAN SST AVERAGES

   ##README
	#http://web1.cdc.noaa.gov/cdc/data.noaa.oisst.v2.html
	# Gridded global SST anomalies from 1856-present derived from UK Met Office SST data which has had sophisticated statistical techniques applied to it to fill in gaps


   ##DOWNLOAD AND SAVE INTO DIRECTORY kaplan.data under name sst.mon.anom.nc
	##ftp://ftp.cdc.noaa.gov/Datasets/kaplan_sst/sst.mon.anom.nc

   ##INSTALL PACKAGE ncdf

   ##FUNCTIONS
	library(ncdf)

	kaplan<-function(lat,long){
		if(long<0) long<-360+long;
		kaplan<-cbind( 1+floor(long/5),18-floor(-lat/5))
		kaplan}
		#simple function to give index from lat-long
 

   ##LOAD DATA
	kaplan.data<-"c:/climate/data/jones/kaplan.sst.v2"		
	Kaplan<-open.ncdf( file.path(kaplan.data,"sst.mon.anom.nc"))
	v1 <- Kaplan$var[[1]]
	sst <- get.var.ncdf( Kaplan, v1 ) 
	dim(sst)
		#  72   36 1808
	    #Spatial coverage:
		 # 87.5S-87.5N (by 5.0), 2.5E-357.5E (by 5.0)
     	    #Temporal coverage:
     		# Monthly anomalies from 01/1856 to 08/2005 (for now)
       	 
	     plot.ts(sst[kaplan(2,150)[1],kaplan(2,150)[2],]) #an example

   #NOW MAKE BASINS FOR WEBSTER ET AL
	flat<-function(lat) { flat<-18-floor(-lat/5);flat}  #gets index from latitude
	flong<-function(long) {temp<-(long<0); long[temp]<- 360+long[temp];flong<-1+floor(long/5);flong}

  	basin.kap<-rep(list(NA),6);names(basin.kap)<-c("NATL","WNP","ENP","SIO","NIO","SWP")
	basin.kap$NATL<-  list(seq(7.5,17.5,5),seq(-87.5, -22.5,5))
	basin.kap$WNP<- list(seq(7.5,17.5,5),seq(122.5, 177.5,5))
	basin.kap$ENP<- list(seq(7.5,17.5,5),seq(-117.5, -92.5,5))
	basin.kap$SIO<- list(seq(-17.5,-7.5,5),seq(52.5, 112.5,5))
	basin.kap$NIO<-list(seq(7.5,17.5,5),seq(57.5, 87.5,5))
	basin.kap$SWP<-list(seq(-17.5,7.5,5),seq(157.5, 177.5,5))

   ##NOW MAKE AVERAGES
	N<-dim(sst)[3]
	X<-array(NA,dim=c(N,6))
	dimnames(X)[[2]]<-names(basin.kap)
	for ( i in 1:6) {
		h<-sst[flong(basin.kap[[i]][[2]]), flat(basin.kap[[i]][[1]]),]
		for (j in 1:N) {X[j,i]<-mean(h[,,j],na.rm=T)}
		}
			#this is an inelegant mean in R but inefficiency is trivial
	X<-ts(X,start=c(1856,1),freq=12)  #dim 1808 6
	plot.ts(X)
	tsp(X) # 1856.000 2006.583   12.000
	write.table(cbind(time(X),X),quote=FALSE,sep="\t",file=file.path(kaplan.data,"basin.average.kaplan.txt"))

  ##PLOT FIGURE 1 # 6-panel plot of averages
	X<-ts(X,start=c(1856,1),freq=12)
	year<-time(X)	

	nf<-layout(array(1:6,dim=c(6,1)),heights=c(1.3,1,1,1,1,1.3))
	par(mar=c(0,4,3,1))
	plot(year,X[,1],type="l",axes=FALSE,ylab=names(basin)[1],main="AVERAGE TEMPERATURES - KAPLAN SST")
	axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0,lty=2)
	par(mar=c(0,4,0,1))
	for (i in 2:5){	plot(year,X[,i],type="l",axes=FALSE,ylab=names(basin)[i])
	axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0,lty=2)}
	par(mar=c(3,4,0,1))
	plot(year,X[,6],type="l",axes=FALSE,ylab=names(basin)[6])
	axis(side=1);axis(side=2,las=1);box();abline(h=0,lty=2)