###########################
####HADCRU
#####################

#GLB, NH, SH: use get_hadcru
#grid: use get_hadcru_grid
#zonal: use subsets of gridded data
#http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.3.0.0.annual_nh.txt
####################################
# NH,SH, GLB from Had Composites
######################################

#12 HAdcru3c.glb #potential target 1961-1990
#http://hadobs.metoffice.com/hadcrut3/diagnostics/index.html: list 
#Global
#    * Mean of northern and southern hemisphere averages (Recommended for general use).
#    * Simple area weighted average.
#    * Comparison of Global Land, Ocean and Blended series.

#Hemispheric
#    * Northern hemisphere.
#    * Southern hemisphere.
#    * Northern-Southern hemisphere difference.

#Regional
#    * Tropics 30S-30N. http://hadobs.metoffice.com/hadcrut3/diagnostics/regional/30-30/monthly
#                      http://hadobs.metoffice.com/hadcrut3/diagnostics/regional/north_30n/annual
#    * Northern extratropics (north of 30N).
#    * Southern extratropics (south of 30S).
#    * Tropics and mid-latitudes 60S-60N.
#    * Central England.


#http://hadobs.metoffice.com/hadcrut3/diagnostics/time-series.html
 #   * Column 1 is the date.
  #  * Column 2 is the best estimate anomaly.
   #  * Columns 11 and 12 are the upper and lower 95% uncertainty ranges from the combined effects of all the uncertainties.

get_hadcru=function(hemi="gl",dset="hadcrut4",version=3,period="monthly",verbose="default")  { 
 if(period=="monthly") {
  if(dset=="hadcrut3") {
	#http://hadobs.metoffice.com/hadcrut3/diagnostics/hemispheric/northern/monthly
	index=c(gl="global/nh+sh",nh="hemispheric/northern",sh="hemispheric/southern",trp="regional/30-30",
	  nh_ext="regional/north_30n", sh_ext="regional/south_30s" )
 	url2=file.path("http://hadobs.metoffice.com/hadcrut3/diagnostics",index[paste(hemi)],"monthly" )  
  	D=read.table(url2) #dim(D) #1928 12 #start 1850
    	name0=c("year","anom","u_sample","l_sample","u_coverage","l_coverage","u_bias","l_bias","u_sample_cover","l_sample_cover",
    	   "u_total","l_total")
    	names(D)=name0
    	hadcru.glb=ts(D[,2],start=c(1850,1),freq=12)	
	    #this is one month more recent than cru version on June 18,2008
	if(verbose=="verbose") return(D)  else return(hadcru.glb)
   }

 if(dset=="hadcrut4") {
	#http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/download.html
	index=c(gl="ns_avg",nh="nh",sh="sh",trp="30S_30N" )
	loc=paste("http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.2.0.0.monthly_",
	    index[paste(hemi)],".txt" ,sep="")  
  	D=read.table(loc) #dim(D) #1928 12 #start 1850
    	name0=c("year","anom","u_sample","l_sample","u_coverage","l_coverage","u_bias","l_bias","u_sample_cover","l_sample_cover",
    	   "u_total","l_total")
    	names(D)=name0
    	hadcru.glb=ts(D[,2],start=c(1850,1),freq=12)	
	    #this is one month more recent than cru version on June 18,2008
	if(verbose=="verbose") return(D)  else return(hadcru.glb)
   }

   if(dset=="cru") {
	url<-paste("http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3v",hemi,".txt",sep="")  #1850 2008 
	h<-scan(url)
	start0<-h[1]
	N<-length(h)
	h<-array(h,dim=c(27,N/27))
	h<-t(h)
	hadcru3v.ann<-ts(h[,14],start=start0) #1850 - 
	hadcru3v=ts(c(t(h[,2:13])),start=c(1850,1),freq=12)
	temp=(hadcru3v==0)&(time(hadcru3v)>=tsp(hadcru3v.ann)[2])
	hadcru3v[temp]=NA
	hadcru3v=window(hadcru3v,end=max( time(hadcru3v)[!is.na(hadcru3v)]) )
	return(hadcru3v)
   }
  } # monthly
 if(period=="annual") {
  if(dset=="hadcrut3") {
	#http://hadobs.metoffice.com/hadcrut3/diagnostics/hemispheric/northern/monthly
	index=c(gl="global/nh+sh",nh="hemispheric/northern",sh="hemispheric/southern",trp="regional/30-30",
	  nh_ext="regional/north_30n", sh_ext="regional/south_30s" )
 	url2=file.path("http://www.metoffice.gov.uk/hadobs/hadcrut3/diagnostics",index[paste(hemi)],"annual" )  
  	D=read.table(url2) #dim(D) #158 12 #start 1850
	}
  if(dset=="hadcrut4") {
	index=c(gl="ns_avg",nh="nh",sh="sh",trp="30S_30N" )
loc=paste("http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.",version,
	".0.0.annual_",  index[paste(hemi)],".txt" ,sep="")  
    	D=read.table(loc) #dim(D) #158 12 #start 1850
	}
    	name0=c("year","anom","u_sample","l_sample","u_coverage","l_coverage","u_bias","l_bias","u_sample_cover","l_sample_cover",
    	   "u_total","l_total")
    	names(D)=name0
    	hadcru.glb=ts(D[,2],start=1850	)
	    #this is one month more recent than cru version on June 18,2008
	if(verbose=="verbose") return(D)  else return(hadcru.glb)
   }#annual

} 


	#  hadcru.glb=get.hadcru()

 #http://hadobs.metoffice.com/hadsst2/diagnostics/global/nh+sh/monthly

#http://hadobs.metoffice.com/hadcrut3/diagnostics/time-series.html
 #   * Column 1 is the date.
  #  * Column 2 is the best estimate anomaly.
  #  * Columns 3 and 4 are the upper and lower 95% uncertainty ranges from the station and grid-box sampling uncertainties.
  #  * Columns 5 and 6 are the upper and lower 95% uncertainty ranges from the coverage uncertainties.
  #  * Columns 7 and 8 are the upper and lower 95% uncertainty ranges from the bias uncertainties.
  #  * Columns 9 and 10 are the upper and lower 95% uncertainty ranges from the combined station and grid-box sampling, and coverage uncertainties.
  #  * Columns 11 and 12 are the upper and lower 95% uncertainty ranges from the combined effects of all the uncertainties.




##########################
##GET GRID
##########################
#FUNCTION TO DOWNLOAD NCDF FROM MET OFFICE 
#dsets HadCRUT4, CRUTEM4
get_grid=function(dset= "HadCRUT4") {

if(dset=="HadCRUT4") {
#loc="http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/gridded_fields/HadCRUT.4.2.0.0.median_netcdf.zip"
#download.file(loc,"d:/temp/temp.zip")
#download.file(loc, "tempnc.zip", "auto", quiet = FALSE, cacheOK = TRUE);
#handle <- unz("d:/temp/tempnc.zip","HadCRUT.4.2.0.0.median.nc" ,"rb") 
#	data <- readBin(handle, "raw", 1024*1024);
#	close(handle);
#	writeBin(data, "d:/temp/temp.nc");
#	rm(data)
#	v <- open.ncdf("d:/temp/temp.nc");
#	instr <- get.var.ncdf( v, v$var[[1]]) # 1850 2006
#	dim(instr)# [1]  72   36 1965
#	return(instr)	

	#loc="http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/gridded_fields/HadCRUT.4.2.0.0.median_ascii.zip"
	loc="http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/gridded_fields/HadCRUT.4.3.0.0.median_ascii.zip"
	download.file(loc,"d:/temp/temp.zip")
	handle=unz("d:/temp/temp.zip","HadCRUT.4.3.0.0.median.txt","r")
	x=scan(handle, what="",sep="\n")
	close(handle)
	n=nchar(x)
	temp= n==122
	x=x[!temp]
	writeLines(x, "d:/temp/temp.dat")
	x= scan("d:/temp/temp.dat")
	N=length(x)
	had=array(x, dim=c(72*36, N/2592) )
	had[had < -99]=NA
	had=ts(t(had),start=1850,freq=12)
	return(had) 
} #N to S


     if(dset=="CRUTEM4") { #N to S 
	loc="http://www.metoffice.gov.uk/hadobs/crutem4/data/gridded_fields/CRUTEM.4.2.0.0.anomalies.txt.gz"
dest="d:/temp/temp.gz"
	download.file(loc,dest,mode="wb")
	handle=gzfile(dest)
	x=readLines(handle); 
	close(handle)
	n=nchar(x)
	M=min( sort(unique(n)))
	x=x[! n==M]
	#x=gsub("-32768","   NA",x)
	writeLines(x, "d:/temp/temp.dat")
	x= scan("d:/temp/temp.dat")
	N=length(x); N/2592
	work=array(x, dim=c(72,36, N/2592) )
	work[work< -999]=NA
	return(work)
}
}


get_hadcru_trp=function(lat=20) {
	require("ncdf")
	hadcru=get_grid(dset="HadCRUT4")
  	  dim(hadcru) #n TO s
	temp= abs(seq(87.5,-87.5,-5))<=lat
	y=apply(hadcru[,temp,],3,mean,na.rm=T)
	y=ts(y,start=1850,freq=12)
	return(y)	
}