# Utility functions for the ship series figures

# Plot a route map
sr.map <-function(lats,longs,times,base_lon=NA,
                  lat_range=NA,
                  lon_range=NA,lon_scale=NA,lat_scale=NA,
                  pch=NA,cex=NA,col=NA,
		  time_range=NA,
                  region=TRUE,pretty=FALSE) {

    # Remove points outside the time range
    if(length(time_range) > 1) {
        for(i in seq(1,length(times))) {
            if(times[i]>time_range[2] || times[i]<time_range[1]) {
	        is.na(lats[i])<-T
		is.na(longs[i])<-T
	    }
	}
    }

    # Guess the base longitude if unspecified
    if(is.na(base_lon)) {
    	range = max(longs,na.rm=T)-min(longs,na.rm=T)
        longs_s = longs
        for(i in seq(1,length(longs_s))) {
            if (!is.na(longs_s[i]) && longs_s[i]< 0) {
                longs_s[i] = longs_s[i]+360
            }
        }
	if(max(longs_s,na.rm=T)-min(longs_s,na.rm=T)<range) base_lon=180
	else base_lon=0
    }
	    
    # Get continental outline data and set to the chosen base longitude
    sr.map.internal.wm <- map('world',interior=FALSE,plot=FALSE)
    is.na(sr.map.internal.wm$x[8836])=T  # Remove Antarctic bug
    longs_s = longs
    if(base_lon != 0) {
        for(i in seq(1,length(sr.map.internal.wm$x))) {
            if(!is.na(sr.map.internal.wm$x[i])) {
                if (sr.map.internal.wm$x[i]< base_lon-180) {
                    sr.map.internal.wm$x[i] = sr.map.internal.wm$x[i]+360
                }
                if (sr.map.internal.wm$x[i]> base_lon+180) {
                    sr.map.internal.wm$x[i] = sr.map.internal.wm$x[i]-360
                }
            }
            if(i>1 && !is.na(sr.map.internal.wm$x[i]) 
                   && !is.na(sr.map.internal.wm$x[i-1])
                   && abs(sr.map.internal.wm$x[i]-sr.map.internal.wm$x[i-1])>50) {
                is.na(sr.map.internal.wm$x[i])=T
            }
        }
        # Correct the data longitudes for the base
        for(i in seq(1,length(longs_s))) {
            if (longs_s[i]< base_lon-180) {
                longs_s[i] = longs_s[i]+360
            }
            if (longs_s[i]> base_lon+180) {
                longs_s[i] = longs_s[i]-360
            }
        }
    }
    # Set default ranges and scales
    if(is.na(pch)) pch="."    # Character to plot
    if(is.na(cex)) cex=1  # Size of character to plot
    if(is.na(col)) col=1    # Colour of character to plot
    if(is.na(lon_range)) lon_range=c(max(base_lon-180,min(longs_s,na.rm=T)-20),
                                     min(base_lon+180,max(longs_s,na.rm=T)+20))
    if(is.na(lat_range)) lat_range=c(max(-90,min(lats,na.rm=T)-20),
                                     min(90,max(lats,na.rm=T)+20))
    if(is.na(lon_scale)) lon_scale = as.numeric(list.nmax.range(seq(-180,360,5),6,
                                                     lon_range[1],lon_range[2]))
    if(is.na(lat_scale)) lat_scale = as.numeric(list.nmax.range(seq(-90,90,5),5,
                                                     lat_range[1],lat_range[2]))
                
    
  mappanel <- function(x,y,...) {
      panel.xyplot(x,y,...)
      llines(sr.map.internal.wm$x,sr.map.internal.wm$y,col="black")
  }

      res<-xyplot(lats ~ longs_s,
         ylab="Latitude",xlab="Longitude",
         xlim=lon_range,
         ylim=lat_range,
         scales=list(x=list(at=lon_scale),
                     y=list(at=lat_scale)
                    ),
         panel=mappanel,
         aspect=(lat_range[2]-lat_range[1])/(lon_range[2]-lon_range[1]),
         region=region,
         pch = pch,
         cex = cex,
         col = col,
         pretty=pretty
      ) 
      
    return(res) 
}


# Plot one or two time-series variables
SeriesPlot <- function(xData,yData,
                       xrange=NA,xtics=NA,xaxis=T,xlabel=NA,
                       yrange=NA,yaxis=T,ylabel=NA,ygp=NA,
                       y2Data=NA,y2range=NA,y2axis=F,y2label=NA,y2gp=NA) {
    
   # Calculate xtics if not defined
   if(!is.list(xtics) && is.na(xtics)) {
       if(!is.list(xrange) && is.na(xrange)) {
           xtics = auto.xtics(min(xData,na.rm=T),max(xData,na.rm=T))
       }
       else {
           xtics = auto.xtics(min(xrange,na.rm=T),max(xrange,na.rm=T))
       }
   }
   
   # Calculate xrange if not defined
   if(!is.list(xrange) && is.na(xrange)) {
       xrange = c(min(xData,na.rm=T),max(xData,na.rm=T))
       if(is.list(xtics) || !is.na(xtics)) {
           xrange[1] = min(xrange[1],min(xaxis.label.to.at(xtics)))
           xrange[2] = max(xrange[2],max(xaxis.label.to.at(xtics)))
       }
   }
   # Discard data outside xrange
   else {
       xmin = min(xrange)
       xmax = max(xrange)
       for(i in seq(1,length(xData))) {
           if(xData[i]<xmin || xData[i]>xmax) is.na(xData[i])=T
       }
   }
   
   # Set default yrange in case all data missing
   if(!is.list(yrange) && is.na(yrange)) {
       if(!is.finite(min(yData,na.rm=T))) { yrange=c(0,1) }
       else {
           yrange = c(min(yData,na.rm=T),max(yData,na.rm=T))
           if(yrange[2]==yrange[1]) yrange[2] = yrange[2]+1
       }
   }
   # Set default graphics parameters for y line
   if(is.na(ygp)) {
       ygp=gpar(col=rgb(0,0,0,1))  #  Black
   }
   
   # If second y axis selected 
   if(is.list(y2Data) || !is.na(y2Data)) {
       # set default y2range
       if(!is.list(y2range) && is.na(y2range)) {
           y2range = c(min(y2Data,na.rm=T),max(y2Data,na.rm=T))
           if(!is.finite(min(y2range,na.rm=T))) { y2range=c(0,1) }
           if(y2range[2]==y2range[1]) y2range[2] = y2range[2]+1
       }
       # Set default graphics parameters
       if(is.na(y2gp)) {
           y2gp=gpar(col=rgb(0,0,1,1))  # Blue
       }
       
   }
      
   # Plot the data
    if(xaxis) {
        pushViewport(plotViewport(c(4,4,1,5)))
    }
    else {
        pushViewport(plotViewport(c(1,4,1,5)))
    }
    pushViewport(dataViewport(xrange,yrange))
    grid.lines(x=unit(xData,"native"),y=unit(yData,"native")
                ,gp=ygp)
    if(xaxis) {
        if(is.list(xtics) || !is.na(xtics)) {
            grid.xaxis(at=xaxis.label.to.at(xtics),label=xtics)
        }
        else {
            grid.xaxis()
        }
    }
    if(!is.na(xlabel)) {
        grid.text(xlabel,y=unit(-3,"lines"))
    }
    if(yaxis) {
        grid.yaxis(main=T)
    }
    if(!is.na(ylabel)) {
        grid.text(ylabel,x=unit(-3,"lines"), rot=90,gp=ygp)
    }
    popViewport()  # dataViewport
   # Add the second y axis data, if selected
    if(is.list(y2Data) || !is.na(y2Data)) {
        pushViewport(dataViewport(xrange,y2range))
        grid.lines(x=unit(xData,"native"),y=unit(y2Data,"native")
                    ,gp=y2gp)
        if(y2axis) {
            grid.yaxis(main=F)
        }
        if(!is.na(y2label)) {
            grid.text(y2label,x=unit(1,"npc")+unit(3,"lines"), rot=90,gp=y2gp)
        }
        popViewport()  # dataViewport    
    }
    popViewport()  # plotViewport
    
}

# Days to the beginning of each month
ltt.private.dim <- c(0,31,59,90,120,151,181,212,243,273,304,334)
ltt.private.dim.leap <- c(0,31,60,91,121,152,182,213,244,274,305,335)

# Convert a string formatted as yyyy/mm/dd to a fraction of a year
date.string.to.year.fraction <- function(Dstring) {
    year <- as.integer(substr(Dstring,1,4))
    month <- 1
    if(nchar(Dstring)>=7) month <- as.integer(substr(Dstring,6,7))
    day <- 1
    if(nchar(Dstring)>=10) day <- as.integer(substr(Dstring,9,10))
    fraction <- 0
    if(is.leap.year(year)) { 
        fraction <- year + (ltt.private.dim.leap[month]+day-0.5)/366
    }
    else {
        fraction <- year + (ltt.private.dim[month]+day-0.5)/365
    }
    return(fraction)
}

is.leap.year <- function(year) {
    if(year/400-as.integer(year/400)==0) return(T)
    if(year/4-as.integer(year/4)==0 && 
       year/100-as.integer(year/100)) return(T)
    return(F)
}

# Calculate a set of locations from a set of labels
xaxis.label.to.at <- function(label) {
    at <- rep(0.0,length(label))
    for(i in seq(1,length(label))) { at[i] <- date.string.to.year.fraction(label[i]) }
    return(at) 
}

# Functions for auto date axis
nearest.year <- function(year) {  # Integer from floating point year
    return(as.integer(year+0.5))
}
nearest.month <- function(year) {  # Integer month from floating point year
    fraction <- year-as.integer(year)
    month = as.integer((fraction+0.0417)*12)+1
    if(month>12) month==12
    return(month)
}
list.nmax <- function(llong,n) {  # decimate a list to have <=n entries
    while(length(llong)>n) {
       lnew = rep('a',0)
       for(i in seq(1,length(llong),by=2)) {
       	  lnew<-append(lnew,llong[i])
       }
       llong <- lnew
    }
    return(llong)
}
list.nmax.range <- function(llong,n,start,end) {
                                  # decimate a list to have <=n entries
                                  # within a set range
    lnew = rep('a',0)
    for(i in seq(1,length(llong))) {
        if(llong[i]<start || llong[i]>end) next
	lnew<-append(lnew,llong[i])
    }
    return(list.nmax(lnew,n))
}
series.years <- function(start,end) { # start and end are floating point years
   series = rep('a',0)
   for(i in seq(nearest.year(start),nearest.year(end))) {
      series <- append(series,sprintf("%04d",i))
   }
   return(series)
}
series.months <- function(start,end) {
   series = rep('a',0)
   for(y in seq(as.integer(start),as.integer(end))) {
      for(m in seq(1,12)) {
          if((y==as.integer(start) && m < nearest.month(start))||
	      y==as.integer(end) && m > nearest.month(end)) next
          series <- append(series,sprintf("%04d/%02d",y,m))
      }
   }
   return(series)
}
auto.xtics <- function(start,end) { # start and end are floating point years
    if(nearest.year(end)-nearest.year(start)>=3) {
    	return(list.nmax(series.years(start,end),5))
    }
    else {
        return(list.nmax(series.months(start,end),5))
    }
}