###This is a single function to read in Data and distribute as needed.

read.extended.list<-function(dir="R:/ike/data/extendedbesttrack",latest="19")                               
{ 
 files<-c("ebtrk1851to1987climoradii",paste("ebtrk",latest,"_atlc",sep=""))
 files<-paste(dir,files,sep="/") 
 files<-paste(files,"txt",sep=".")                             
 out<-NULL  ##empty dataframe
 for(file in files)
  {
   print(file)
   tmp<-read.extended(file)  ##read in new data
   out<-rbind(out,tmp)       ##combine data
  }
outd<-rev(duplicated(out[(nrow(out)):1,c("Basin","Sn","SYear","name","Mo","Da","hr","Yr")]))  
out<-out[!outd,] ###reduce data set
out <- addHour(out,extended=T)
out
}

read.ff<-function(file,fields,as.is=T,...)
  {
   ###This program provides an alternative way to read fixed length files, using a named list
   ###Each named element consists of a single vector listing the elements
   if(any(diff(unlist(fields))<=0)) stop("Fields must be in order")
   lower<-sapply(fields,function(x) x[1])
   upper<-sapply(fields,function(x) rev(x)[1])
   mp<-rep(c(-1,1),length(lower))
   widths<-mp*(diff(c(0,as.vector(rbind(lower,upper))))+mp)
   widths<-widths[widths!=0]
   fieldnames<-names(fields)
   read.fwf(file=file,widths=widths,col.names=fieldnames,header=F,as.is=as.is,...)
  }


read.extended<-function(file, fields=list(Basin=1:2, Sn=3:4, SYear=5:6, name=8:17, Mo=18:19, Da=20:21, hr=22:23,
                           Yr=25:28, lat=29:33, lon=34:39, Wmax=41:43, pmin=45:48, Rmax=50:52,
                               Dey=54:56, poci=58:61, Roci=63:65, 
                               r34.NE=67:69,r34.SE=70:72,r34.SW=73:75,r34.NW=76:78,
                               r50.NE=80:82,r50.SE=83:85,r50.SW=86:88,r50.NW=89:91,
                               r64.NE=93:95,r64.SE=96:98,r64.SW=99:101,r64.NW=102:104,
                               ET=106:106)
                               ,
                          lg=landGrid,extended.track=T
                          )
{
  data<-read.ff(file=file,fields=fields,na.strings=c(" -99","  -99","-99"))
  
  ##Reads in extended data set, and outputs best track data set
   if(extended.track)
    {
  ##Update track data if this is for track data
      data$lon<- -data$lon + 360*(data$lon > 180)
      data$SYear<-data$SYear+100*data$Yr%/%100  ##%ANY% takes precedence over * 
       data$pmin[data$pmin<=0]<-NA
      ### Compute  meanRmax in nautical miles using the regression formula derived using km.
      data$meanRmax = 0.539593*exp(2.636-5.086e-005*(1013-data$pmin)^2 + 0.0394899*data$lat)

      data$M<-is.land(points=data[,c("lon","lat")],grid=lg)
     }
  return(data)
}

"is.land"<-function(grid=landGrid,points,center=NULL,method="l")
{
###This function determines whether a vector of points is inside hull
###by considering the number of rotations a point on the curve makes about a single point.
###The result counts is this number it should be the number of counter clockwise rotations.
###The sign denotes the curve direction (+ couterclockwise) about an interior point.
###This function was converted from a C program
###The center is chosen so that the grid comprises a simple closed curve
###The result is a vector with the number of track crossings
###This works for all closed curves
###Method=1 use lapply, slower but works on large data sets
###Method=2 use outer3, faster, but works only on single tracks.
  if(!is.null(center))
    grid<-rbind(center,grid)
  xp<-grid[,1]
  yp<-grid[,2]
  x<-points[,1]
  y<-points[,2]
  if(any(colnames(points)[1:2] != colnames(grid)[1:2]))
    warning(paste("points column names:",colnames(points)[1:2],"are not equal to","grid column names:",colnames(grid)[1:2]))
  lx<-length(xp)
  ##Do not add last point, remove first point instead assume that each loop is closed by having the first and last points
  ##duplicated. (use close.curves)
  xpi<-xp[-1]
  xpj<-xp[-lx]
  ypi<-yp[-1]
  ypj<-yp[-lx]
  gridcross<-cbind(xpi,xpj,ypi,ypj)
  xy<-cbind(x,y)
  ##get x crossings
  if(method=="l")
    {
  crossings<-unlist(lapply(1:nrow(xy),
                           function(i,xy,grid)
                           {
                             x<-xy[i,1]
                             y<-xy[i,2]
                             xpi<-grid[,1]
                             xpj<-grid[,2]
                             xi<-x < xpi
                             xj<-x < xpj
### xcross=(xj&(!xi)) | (xi&(!xj) )  ==
                             vi<-(xj-xi)
                             vi[is.na(vi)]<-0 ###Closed Curves separated by NAs
                             vil<-as.logical(vi)
                             if(!any(vil)) return(0)
                             vi<-vi[vil]
                             grid<-grid[vil,]
                             xpi<-grid[,1]     
                             xpj<-grid[,2]
                             ypi<-grid[,3]
                             ypj<-grid[,4]
                             xdelta<-xpi-xpj
                             ydelta<-ypi-ypj
                             sum(vi*(y<(ypi+(x-xpi)*ydelta/xdelta)))
                           },xy=xy,grid=gridcross))
  is.inside<-as.logical(crossings)
}
else
  {##Keep for later, may be faster
  ydelta<-ypj-ypi
  xdelta<-xpj-xpi
   gridcross<-cbind(xdelta,ydelta,xpi,ypi)
   xj<-outer(x,xpj,"<")
   xi<-outer(x,xpi,"<")
  ## xcross=(xj&(!xi)) | (xi&(!xj) )  ==
  xcross<-(xj-xi)
  xcross[is.na(xcross)]<-0 ###Closed Curves separated by NAs
   crossings<-outer3(x=xy,y=gridcross,fun=function(xy,grid,v)
                    {
                      v<-as.vector(v)
                      x<-xy[,1]
                      y<-xy[,2]
                      xdelta<-grid[,1]
                      ydelta<-grid[,2]
                      xp<-grid[,3]
                      yp<-grid[,4]
                      ifelse(v,v*(y<(yp+(x-xp)*ydelta/xdelta)),0)
                    },v=xcross
                    )
  is.inside<-as.logical(rowSums(crossings))
}
#  print(all(is.insideo==is.inside))
  return(is.inside)
}
##C code:
##  int pnpoly(int npol, float *xp, float *yp, float x, float y)
##    {
##      int i, j, c = 0;
##for (i = 0, j = npol-1; i < npol; j = i++) {
##  if ((((xp[i]<=x) && (x<xp[j])) ||
##       ((xp[j]<=x) && (x<xp[i]))) &&
#      
###/*where is the y crossing value?*/
##                   
##      (y < (yp[j] - yp[i]) * (x - xp[i]) / (xp[j] - xp[i]) + yp[i]))
##          c = !c;
##      }
##      return c;
##    }

make.cov<-function(data,months,start,end)
  {
    ###This function takes each data set and finds the mean for the selected months
    ###data is a named list of data sets
    ###start and end are the starting and ending years
    ld<-length(data)
    if(ld!=length(months)) stop("Length mismatch between data and months")
    se<-start:end
    aves<-do.call("data.frame",c(list(se),rep(list(NA),ld)))
    dimnames(aves)[[2]]<-c("Yr",names(data))
    for(i in 1:length(data))
      {
        x<-data[[i]]
        m<-months[[i]]
        y<-as.vector(rowMeans(x[,m,drop=F]))
        yr<-x$Yr
        aves[match(yr,se,nomatch=0),i+1]<-y[yr%in%se]
      }
    return(aves)
  }

readIntensity<-function(file)
{
data<-read.csv(header=F,file)
count<-apply(data,1,function(x) sum(!is.na(x)))
Wmax<-as.vector(t(data[,-1]))
return(data.frame(Yr=rep(data[,1],count-1),Wmax=Wmax[!is.na(Wmax)]))
}
readIntensity2<-function(file,Yr="SYear",Wmax="WmaxS")
{
data<-read.csv(header=T,file)
return(data.frame(Yr=data[[Yr]],Wmax=data[[Wmax]]))
}

readCov<-function(file,last=F,header=F,alternate=F,miss=NULL,ma=c("","c","cs","l"))
  {
###This function generates data as needed
###If year is greater than one we generate the past years data as well
###with fields jan.last, feb.last,.....,
###
###Data may have header, but we will NOT use it...
###Input must have Yr field followed by 12 or 13 columns
###ma="", do not center or scale, ma="c", center, ma="cs" both, ma="l", last row is subtracted out

    data<-read.table(header=header,file,fill=T)
    if(alternate)
    {
    data<-data[seq(1,nrow(data),2),]
    }
    if(length(data)>14) stop("Invalid Data")
    body<-data[,2:13]
    if(!is.null(miss))
      body[apply(body,2,function(x,y) x%in%y,y=miss)]<-NA
    Yr<-data[,1]
    nYr<-length(Yr)
    ma<-match.arg(ma)
    if(ma!="")
     { 
      if(ma=="l")
      {
       body<-(body-body[rep(nYr,nYr),])[-nYr,]
       Yr<-Yr[-nYr]
       nYr<-nYr-1       
      }
      else
      {
       if(ma=="cs") scale=T else scale=F
       data<-scale(data,center=T,scale=scale)
      }
    }    
    if(last)
      {
        body<-cbind(rbind(rep(NA,12),body[-nYr,]),body)
        dimnames(body)[[2]]<-c(paste(month.abb,"last",sep="."),month.abb)      
      }
    else
      dimnames(body)[[2]]<-month.abb
    data.frame(Yr=Yr,body)
  }
    
get.data<-function(directory,stem,colstem=c("H","I"),colloc=1:length(colstem)+1,start=1871, end=2005)
  {
    if(length(colstem) != length(colloc)) stop("Colstem length must match number of selected columnes")
    ###colloc, is the column locations, defaults to 2,3 
    ###This function retreives data from point files, and generates columns for each point.
    ###The data format is assumed to be in the form of Year Hurricane and Intense
    #directory is the location of the files
    #stem is the file stem files have names stem_1_counts.txt
    files<-list.files(directory,pattern=paste(stem,".*","[.]txt",sep=""))
    if(length(files)==0) stop("No files in directory")
    startstem<-nchar(stem)
    endstem<-nchar(files)
    subname<-substring(files,startstem+1,endstem)
    values<-regexpr("[0-9]+",subname)
    length<-attr(values,"match.length")
    names<-substring(subname,first=values,last=values+length-1)
    use<-nchar(names)>0
    files<-files[use]
    names<-names[use]
    names<-as.numeric(names)
    files<-files[order(names)]
    fullfiles<-paste(directory,files,sep="/")
    
    names<-paste(rep(paste("P",format(sort(names)),sep="."),each=length(colstem)),rep(colstem,length(names)),sep=".")
    names<-gsub(" ",0,names) ###get correct names...
    data<-lapply(fullfiles,
                 function(x,s,e)
                 {
                   x<-read.table(x,header=F)
                   return(x[x[,1]>=s & x[,1]<=e,colloc,drop=F])
                 },
                 s=start,e=end)    
    testdata<-sapply(data,nrow)
    nrowt<-testdata != end-start+1
    if(any(nrowt))
      stop(paste("File(s):",paste(fullfiles[nrowt],sep=","),"has incorrect number of years, either missing or duplicates",sep=" "))
    testdata<-do.call("cbind",data)
    colIds(testdata)<-names
    testdata
  }

##simple support functions for extendedBest
vt1<-function(x) apply(x[c("Yr","Mo","Da","hr","lon","lat")],1,paste,collapse=":")  #Function to get the unique observational identifier from each
vt11<-function(x) apply(x[c("Yr","Mo","Da","hr")],1,paste,collapse=":")
vt2<-function(x) apply(x[c("Yr","Mo","Da","hr","Sid")],1,paste,collapse=":")

  
extendedBestSid<-function(x=extended.track,y=use.best)
{
#This function is used to provide a new StormId
#iter1
xx<-vt1(x)
yy<-vt1(y) ##make id vector
xs<-split(xx,x$Sid)
index<-split(1:length(xx),x$Sid)
xm<-lapply(xs,function(x) y[match(x,yy,nomatch=0),"Sid"]) #If no match we dont worry
if(!all(sapply(xm,function(x) all(x==x[1])))) stop("Duplicate storm ID (Sid) in extended.track file iter 1")
Sid1<-sapply(xm,function(x) x[1]) #first ids
id<-which(is.na(Sid1))
y1<-y[!(y$Sid %in% Sid1[-id]),] #remove storms that we have already accounted for 
#iter 2
xx1<-vt11(x)
yy1<-vt11(y1)
xs1<-split(xx1,x$Sid)[id] #Reduce xs1 to "bad" sets 
xm1<-lapply(xs1,function(x) y1[match(x,yy1,nomatch=0),"Sid"])
if(!all(sapply(xm1,function(x) all(x==x[1])))) stop("Duplicate storm ID (Sid) in extended.track file iter 2")
Sid2<-sapply(xm1,function(x) x[1])
if(any(is.na(Sid2))) print(paste(names(xm1)[which(is.na(Sid2))],"Track is missing from best track")) 
Sid1[id]<-Sid2 #Swap in the ids that were zero
Sid<-rep(Sid1,sapply(index,length))
Sid[unlist(index)]<-Sid #There will be some NAs because one storm was removed in 1916
x$Sid<-Sid
xx2<-vt2(x)
yy2<-vt2(y)
xintoy<-match(xx2,yy2,nomatch=0)
x<-x[xintoy>0,]
xold<-x[,c("Sn","lat","lon","Wmax")] #Keep old data
xremove<-c("Sn","SYear","name","Yr","Mo","Da","hr","lon","lat","Wmax","ET","Shour","Sid","Type")
x<-x[,-(match(xremove,colnames(x),nomatch=0))]  
colnames(xold)<-paste(colnames(xold),"old",sep=".")
cbind(y[xintoy,],x,xold)
}

###Support IBTRACS
makeyearmonthday<-function(x){ xx<-as.POSIXlt(x); data.frame(Yr=xx$year+1900,Mo=xx$mo+1,Da=xx$mday,hr=xx$hour)}
stripone<-function(x) substring(x,2,nchar(x))
worldShours<-function(Sid,ISOtime){if(is.character(ISOtime)) ISOtime<-as.POSIXct(ISOtime,tz="UTC")
  hours<-as.numeric(ISOtime)/3600
  timesplit<-split(hours,Sid) 
  times<-unlist(lapply(timesplit,function(x) x-x[1]))
  times[unlist(split(1:length(Sid),Sid))]<-times
   times}
   
createworlduse<-function(x,conv1to10=.88)
{ 
 #Wmax in 10 minute sustained converted to 1 minute sustained (estimated by dividing by .88)
  xx<-x[,c(1,2,3,4,5,6,7,8,9,10,13,20)]
  colnames(xx)<-c("IBtracSid","SYear","Sn","Basin","SubBasin","name","ISOtime","Type","lat","lon","Wmax","Pmin")
  xx$Basin<-stripone(xx$Basin)
  xx$SubBasin<-stripone(xx$SubBasin)
  xx<-cbind(xx,makeyearmonthday(xx$ISOtime))
  xx[xx==-1]<-NA
  xx$Wmax<-xx$Wmax/conv1to10
  xx$Sid<-as.numeric(factor(xx$IBtracSid))
  xx$ISOtime<-as.POSIXct(xx$ISOtime,tz="UTC")
  xx$Shour<-worldShours(xx$Sid,xx$ISOtime)
  xx
}

importLandSeaGrid<-function(con="http://snowdog.larc.nasa.gov/surf/data/water_perc.bin")
{
 waterperc<-readBin(con=con,n=2160*1080,size=1,what="integer",signed=F,endian="Big")
 dim(waterperc)<-c(2160,1080)
 waterperc<-waterperc[1:2160,1080:1]
 longitude<-((1:2160)-.5)/6
 latitude<-((1:1080)-.5)/6-90
 dimnames(waterperc)<-list(longitude=round(longitude,1),latitude=round(latitude,1))
 attr(waterperc,"longitude")<-longitude
 attr(waterperc,"latitude")<-latitude
 waterperc
}


#waterperc<-importLandMask()
#t(waterperc[1650:1660,725:715]) #big bend Florida
#require(fields)
#image.plot(waterperc,col=rev(tim.colors(100)),y=attr(waterperc,"latitude"),x=attr(waterperc,"longitude"))

getLandSeaPerc<-function(x=world.interp,lg=waterperc)
{
#hardcodedd 
index<-cbind(lon=trunc((x$lon*6)%%(6*360))+1,lat=trunc((x$lat+90)*6)+1)
perc<-lg[index]
return(perc)
}