add-continents.r

## add-continents.r
## Defines function add.continents() which adds a quick, low-level
## outline of the world's continents to an existing plot
## (i.e. no need to fuss around with special map classes,
## but the resulting outline is coarse)

## if use.ndc = TRUE, will draw continents over the entire figure
## device
## -------------------------------------------------------
## Author: Laura Tremblay-Boyer (l.boyer@fisheries.ubc.ca)
## Written on: August 29, 2013
## Time-stamp: <2013-08-29 14:17:45 Laura>

add.continents <- function(..., use.ndc=FALSE, lonlim=c(-180,360), latlim=c(-90,90)) {

## Extract data for continents
## Datasets rmp and gwm from now defunct R package sspline
    if(!exists("rmp")) {
        scd <- getwd()
        setwd("~/Projects/misc-ressources/quick-world-map/")
        rmp <<- read.table("rmp.txt",header=TRUE) #saving objects to global env
        gwm <<- read.table("gwm.txt",header=TRUE)
        setwd(scd)} # back to original WD

    ind1 <- rmp$inc[-1] # rmp gives the indices for individual polygons in gwm
    gwm$ID <- rep(1:length(ind1), ind1) # assign polygon ID to each lat/lon
    gwm.split <- split(gwm, gwm$ID)

    # only plot polygons that are in the current plot frame
    # as defined by par("usr")
    # ... now add lines to current plot
    if(use.ndc) {par(new=TRUE, mfrow=c(1,1))
                 plot(1,type="n",xaxs="i",yaxs="i",xlim=lonlim, ylim=latlim, ann=FALSE, axes=FALSE)}

    check.if.in <- function(poly,plus360=0) {

        pu <- par("usr")

        # include polygon is any edge occurs in the current plot
        p.in <- all(any((poly$lon + plus360) %between% pu[1:2]),
              any(poly$lat %between% pu[3:4]))

        # if the polygon is in, return coordinates and add NA to create
        # a cut when the lines() function is used
        # (could also loop through polygons, but this is WAY faster)
        if(p.in) {return(cbind(c(poly$lon+plus360,NA), c(poly$lat,NA)))}
    }

    # test if polygons are in:
    poly.in <- lapply(gwm.split, check.if.in)
    pmat <- do.call(rbind, poly.in)
    # test if polygons are in using shifted longitudes
    # (thank you South Pacific)
    poly.in360 <- lapply(gwm.split, check.if.in, plus360=360)
    pmat360 <- do.call(rbind, poly.in360)


        lines(pmat, ...)
        lines(pmat360, ...)



}

add.continents.NDC <- function(..., lonlim=c(-180,360), latlim=c(-90,90)) {

    # redefine plot window:
    par(mfrow=c(1,1))
    plot(1,type="n",xaxs="i",yaxs="i",xlim=lonlim, ylim=latlim, ann=FALSE, axes=FALSE)

    # add world map
require(mapdata)

     # EEZ land masses to keep for map --
     eez2keep <- c("American Samoa", "Australia", "Canada","China",
                   "Cook Islands",
"Fiji", "French Polynesia", "Guam", "Hawaii", "Indonesia","Japan",
"Kiribati", "Marshall Islands", "Mexico", "Micronesia",
"Nauru", "New Caledonia", "New Zealand", "Niue", "Northern Mariana Islands",
"Palau", "Papua New Guinea", "Philippines", "Samoa", "Solomon Islands",
"Tokelau", "Tonga", "Tuvalu", "USA", "USSR", "Vanuatu",
                   "Panama","Chile","Argentina", "New Caledonia",
                   "Belize","Nicaragua","Ecuador","Honduras","Costa Rica",
                   "Colombia", "Uruguay", "Brazil", "Peru", "Guatemala",
                   "Venezuela","Bolivia","Paraguay")

#sa <- map('world2Hires', namesonly=T, ylim=c(-50,20), xlim=c(260,300), plot=F)
#    eez2keep <- c(eez2keep, sa)
     data(world2HiresMapEnv) # loads dataset with landlines
     spc.region <- map('world2Hires', namesonly=T, ylim=c(-45,50), plot=F)

     # Keep SPC countries + Indonesia (to show with PNG)
     geo <- lapply(eez2keep,  function(ee) spc.region[grep(ee,spc.region)])

      # ... add land for SPC countries (defined in object "geo")
        map('world2Hires', regions=unlist(geo), add=T,
            wrap=TRUE,fill=TRUE, ...)





}