###HadISST

###This scapes TRP data

###################################
##  collate and save \\
### download from http://hadobs.metoffice.com/hadisst/data/HadISST1_SST_ 
################
## SOURCES
#################

## http://badc.nerc.ac.uk/data/hadisst/file_format.html#ice
## http://badc.nerc.ac.uk/data/hadisst/file_format.html#sst
## KNMI	http://climexp.knmi.nl/data/ihadisst1_0-360E_-20-20N_na.txt
##
## http://hadobs.metoffice.com/hadisst/data/HadISST_sst.nc.gz"  very large
##
## download from  "http://hadobs.metoffice.com/hadisst/data/HadISST1_SST_2008.txt.gz" and organize to 
## d:/climate/data/gridcell/hadisst/sst",index_had[i],".nc"
## ## SCRAPE HADISST DATA FROM 1970s ON

##http://hadobs.metoffice.com/hadisst/data/HadISST1_prelim.nc


##COLLATE HISTORY

get.hadisst=function(period) {
   	loc=paste("http://hadobs.metoffice.com/hadisst/data/HadISST1_SST_",period,".gz",sep="")
	download.file(loc,"temp.gz",mode="wb")
	handle=gzfile("temp.gz")
     	Data <- readLines(handle,n=-1); 
	close(handle)
	N=length(Data);N #2172
	index=grep("columns",Data) #[1]   1 182
	temp=is.na(match(1:N,index))
	Data=Data[temp]
	Data=gsub("-"," -",Data)
	   range(nchar(Data)) #[1] 2160 2520
	   # 2520 = 360* 7. 
	   # 2160/12 = 180 
	writeLines(Data,"temp.dat")
	x=scan("temp.dat");
	N=length(x);
	N/(360*180) 
	x=array(x,dim=c(360,180,N/360/180)) 
	  #360 180 120 # W to E from dateline; N to S
	temp= (x== -32768)
	x[temp]=NA
	      #image.plot (seq(-180,179,1),seq(-89.5,89.5,1),x[,180:1,1])
		#this is sensible rendering
	temp=  abs( seq(-89.5,89.5,1))<20; sum(temp)  #40

	sst=x[,temp,];dim(sst) # 360 40 12  #long= 360 lat = 40 time= 120 
	if(N/(360*180)>1) {sst=aperm(sst,perm=c(3,1,2));
	      sst=array(sst,dim=c(dim(sst)[1],40*360 ) )} else{
 	        sst=array(sst,dim=c(1,40*360) ) }
	dim(sst) #12 14400 
	return(sst)
	}


## DATA
# gridcell/hadisst/sst/1980.1989.ascii.tab etc: TRP data organized yearx 14400: hour hand lat: minute hand long.
#   collated using make.hadisst- scrapes data and store sst.tab(TRP) in directory 
#
# gridcell/hadisst/hadisst.trp.anom.tab - climatology 
#   #made using function 
# 
#  index_had:  nomenclature of zip files
#
#  preliminary: monthly preliminary data scraped separately
#
#  trp_avg - utility function to calculate TRP average for time period
#
# Update procedure:
#  0. load climatology: already scraped
#  1. load history
#  2. get 2009 over again and collate to history
#  3. get preliminary and collate to history
#  4.  calculate anomaly-average series

###################
###UPDATE STEPS
######################

make.hadisst=function(method="update") {
   	index_had=c(paste(c("1960-1969","1970-1979","1980-1989","1990-1999","2000-2003"),"ascii",sep="."),
            paste(2004:2009,"txt",sep="."));
	K=length(index_had)
  if(method=="update") {
      #this option uses prior collation, refreshes 2009 and updates 2010
	download.file("http://www.climateaudit.info/data/gridcell/hadisst/climatology.trp.tab","temp.dat",mode="wb");load("temp.dat")
    	   #0. load climatology # see history 
	  dim(climatology.trp)
	    #  12 14400
	download.file("http://www.climateaudit.info/data/gridcell/hadisst/hadisst.trp.tab","temp.dat",mode="wb");load("temp.dat")
	   #this is previously collated history
	   tsp(hadisst.trp)
	   	#[1] 1970.000 2009.417   12.000
      #2. update 2009 
 	sst= get.hadisst(period="2009.txt")/100 #saves the scaped value
	   dim(sst) # 12 14400
 	month=(1:nrow(sst))%%12; month[month==0]=12
	    month
	trp_avg=function(x) mean( x[2:length(x)] - climatology.trp[ x[1],],na.rm=T )
	    #function to calculate TRP avg: requires binding month to SST
	update_2009= ts(apply(cbind(month,sst),1,trp_avg),start=2009,freq=12)
	    #updated 2009 
      #3 update prelim 
	sst= get.hadisst(period="update.txt")/100 #saves the scaped value
	  dim(sst) # 12 14400
 	month=(1:nrow(sst))%%12; month[month==0]=12
	   month
	update_new= ts(apply(cbind(month,sst),1,trp_avg),start=2010,freq=12)
    #collate extension to history
	hadisst.trp = ts(c( window(hadisst.trp,end=2008.99), update_2009,update_new),start=tsp(hadisst.trp)[1],freq=12)
    }  #end update method
	#  save(hadisst.trp,file="d:/climate/data/gridcell/hadisst/hadisst.trp.2010_01.tab")
}