# R script to show a ship route

library(maps)
library(lattice)
         
sr.map <-function(lats,longs,times,base_lon=0,
                  palette="diverging",ncols=17,levels=NA,lat_range=c(-90,90),
                  lon_range=NA,lon_scale=NA,lat_scale=NA,
                  pch=NA,cex=NA,col=NA,
                  region=TRUE,pretty=FALSE) {

    # Set default ranges and scales
    if(is.na(lon_range)) lon_range=c(base_lon-180,base_lon+180)
    if(is.na(lon_scale)) lon_scale = c(-180,-135,-90,-45,0,45,90,135,180,224,270,315,360)
    if(is.na(lat_scale)) lat_scale = c(-90,-60,-30,0,30,60,90)
    if(is.na(pch)) pch=1    # Character to plot
    if(is.na(cex)) cex=0.2  # Size of character to plot
    if(is.na(col)) col=1    # Colour of character to plot
    
    # Get continental outline data and set to the chosen base latitude
    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
            }
        }
    }
                
    

  mappanel <- function(x,y,...) {
      panel.xyplot(x,y,...)
      llines(sr.map.internal.wm$x,sr.map.internal.wm$y,col="black")
  }

      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="iso",
         region=region,
         pch = pch,
         cex = cex,
         col = col,
         pretty=pretty
      )  

}