##SEAICE FUNCTIONS
  #the hole is coded as 251; 
  #the hole function counts the number of cells with code 251
  #the area function counts cells at 250 and under, thus EXCLUDING the hole;
    #requires a matrix Cells separately collated for cell areas
	 #multiplies Cell area by percentage cover where pct exceeds 15
	 #discontinuity when "hole" changes in 1986-87
  #the extent function counts cells with code 251 and under, with more than 15%
	#thus INCLUDES hole at 100%

##COMMENT the extent in 1987 is nearly identical to the extent in 1986 with changed hole
 #but the area increased by 5139-4584= 555 which is approx amount of discontinuity perhaps 

# tapply(Arctic$area,Arctic$year,min,na.rm=TRUE)
# 1979  1980  1981  1982  1983  1984  1985  1986  1987  1988  1989  1990  1991  1992  1993  1994  1995  1996  1997  1998  1999 
#4.530 4.644 4.193 4.294 4.491 3.887 4.148 4.584 5.139 4.999 4.698 4.492 4.331 4.906 4.314 4.665 4.277 5.123 4.793 4.124 4.074 
# 2000  2001  2002  2003  2004  2005  2006  2007  2008 
#4.030 4.405 3.888 3.982 4.146 3.920 3.907 2.793 4.243 
#  tapply(Arctic$extent,Arctic$year,min,na.rm=TRUE)
# 1979  1980  1981  1982  1983  1984  1985  1986  1987  1988  1989  1990  1991  1992  1993  1994  1995  1996  1997  1998  1999 
#7.844 8.396 7.816 8.009 8.168 7.238 7.457 8.070 8.079 8.193 7.920 7.053 7.300 8.244 7.322 8.057 6.961 8.170 7.432 7.173 6.689 
# 2000  2001  2002  2003  2004  2005  2006  2007  2008 
#6.934 7.618 6.498 6.875 6.537 6.151 6.774 4.969 7.149 #
# tapply(Arctic$hole,Arctic$year,min,na.rm=TRUE)
#  1979   1980   1981   1982   1983   1984   1985   1986   1987   1988   1989   1990   1991   1992   1993   1994   1995   1996 
#1.1926 1.1926 1.1926 1.1926 1.1926 1.1926 1.1926 1.1926 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 
#  1997   1998   1999   2000   2001   2002   2003   2004   2005   2006   2007   2008 
#0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 

##CELLS 

##POLAR STEORTGTAPH http://nsidc.org/data/grids/ps_grid.html
#The polar stereographic projection specifies a projection plane or grid tangent to the Earth at 70 degrees. The planar grid is designed 
#so that the grid cells at 70 degrees latitude are 25 km by 25 km.
 #url="ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/tools/pss25area_v2.dat issomething different

update.nsidc=function() {
	 N0=1+max( (1:nrow(Arctic))[!is.na(Arctic$hole)]);N0
	 N1=(1:nrow(Arctic))[Arctic$julian==(julian(Sys.Date())-1)   ];N1
	 for (k in N0:N1) {
	   Data=scrape.arctic(day=Arctic$dd[k],month=Arctic$mm[k],year=Arctic$year[k],sat=Arctic$sat[k])
	   dest.file=paste(Arctic$year[k],Arctic$mm[k],Arctic$dd[k],Arctic$sat[k],sep="_")
	   save(Data,file=file.path("d:/climate/data/seaice/nsidc",paste(dest.file,"tab",sep=".") ))
   	   if ( !class(Data)=="logical") {
	     Arctic$hole[k]=hole(Data) 
	     Arctic$extent[k]=extent_seaice(Data) 
	     Arctic$area[k]= area_seaice(Data) 
	    }
	  } #k
	update.nsidc=Arctic
	update.nsidc
	}
	
	#Arctic=update.nsidc()

 make.Cells=function(){
	url="ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/tools/psn25area_v2.dat"
	download.file(url,"temp.bin",mode="wb")
	handle=file("temp.bin","rb")
   	X=readBin(handle, "int", 448*304, size=4, signed = FALSE)
	close(handle)
	make.Cells=array(X,dim=c(304,448))/1000
	make.Cells=make.Cells[,448:1]
	make.Cells
	}

	Cells=make.Cells() #304 448
 	#save(Cells,file="d:/climate/data/seaice/Cells.tab") #dim 304 448

	#par(mar=c(3,3,2,1))
	#image.plot (1:304,1:448,Cells)

#Distortion in the grid increases as the latitude decreases, because more of the Earth's surface falls 
#into any given grid cell, which can be quite significant at the edge of the northern SSM/I grid 
#where distortion reaches 31 percent. 

  extent_seaice = function(X,method="N") {	
	temp=(X>251);X[temp]=NA
	temp2=(X>=.15)&!is.na(X) #this includes hole at 100%
	extent_seaice=sum( Cells[temp2])/1E6
	extent_seaice
	}
  #    extent_seaice(Data) #  10.0122


  hole = function(X,method="N") {	
	temp=(X==251);
	hole=sum( Cells[temp])/1E6
	hole
	}
#	hole(Data) # ] 0.3107784



   area_seaice = function(X,method="N") {	
	temp=(X>250);X[temp]=NA
	X[!temp]=X[!temp]/250
	temp2=(X>=.15)&!is.na(X)
	area_seaice=sum( (Cells*X)[temp2])/1E6
	area_seaice
	}
#	area_seaice(Data) # 6.93482


 scrape.arctic=function(day,month,year=2008,sat="override") {
	if(sat=="override") sat=sat_select(DD)
	if (year==2008) {loc= "ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north" ;
	  orig.file=paste("nt_2008",pastedate(month),pastedate(day),"_",sat,"_nrt_n.bin",sep="") 
	url=file.path(loc,orig.file) } else {
	  loc=file.path("ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/final-gsfc/north/daily",year)
	  orig.file=paste("nt_",year,pastedate(month),pastedate(day),"_",sat,"_v01_n.bin",sep="")
	  url=file.path( loc,orig.file)}
	dest.file=paste(year,month,day,sat,sep="_")
	test =try( url(url,"rb"))
	if( length(class(test))==1) scrape.arctic=NA else {
 	  download.file(url,"d:/temp/temp.bin",mode="wb")
	  unlink(url);
	  handle=file("d:/temp/temp.bin","rb")
	  header= readBin(handle, "int", 1*300, size=1, signed = FALSE, endian = "big")
	  X = readBin(handle, "int", 448*304, size=1, signed = FALSE, endian = "big")
	  unlink(handle);close(handle)
	  X=array(X,dim=c(304,448))
	  scrape.arctic=X[,448:1]
	  }
	scrape.arctic
	}

 pastedate=function(x) {if (nchar(x)==1) pastedate=paste("0",x,sep="") else pastedate=paste(x);pastedate}  
	#day=Arctic$dd[k];month=Arctic$mm[k];year=Arctic$year[k];sat=Arctic$sat[k]