R utilties to convert oce::coastlineWorld and ocedata::coastlineWorld* coastlines to sp::SpatialLines

oce_coastline.R

# coastline.R

bb <- c(-72,-63,39,46)
#' Convert a bbox to sp::polygon, sp::Polygons, or sp::SpatialPolygons object
#' 
#' @export
#' @param bb numeric, a 4-element bbox vector [left, right, bottom, top]
#' @param projstring character sutiable to pass to \code{sp::CRS}, 
#'    by default "+proj=longlat +datum=WGS84"
#' @param id character, the polygon ID, by default 'bbox'
#' @param output_class character, either "SpatialPolygons", "Polygons" or "Polygon"
#' @return sp::SpatialPolygons object or NULL
bbox_to_polygon <- function(bb, 
   projstring = "+proj=longlat +datum=WGS84",
   id = 'bbox',
   output_class = c("SpatialPolygons", "Polygons",  "Polygon")[1]){

   stopifnot(require(rgdal))
   stopifnot(require(sp))

   if(missing(bb)) stop("bounding box, bb, is required")
   stopifnot(length(bb) == 4)
   
   # clockwise to form an 'island'
   bb <- cbind(
      x = c(bb[1], bb[1], bb[2], bb[2], bb[1]),
      x = c(bb[3], bb[4], bb[4], bb[3], bb[3]) )
   # make a Polygon
   Poly <- sp::Polygon(bb)
   if (output_class == 'Polygon') return(Poly)
   # make into Polygons
   Polys <- sp::Polygons(list(Poly), id)
   if (output_class == 'Polygons') return(Polys)
   SpatialPolygons(list(Polys), proj4string = CRS(projstring))
}


#' Divide the provided oce::coastlineWorld* object into segements
#'
#' @export
#' @param x an ocedata::coastlineWorld* object (coarse, medium or fine)
#' @param min_length the minimum number of points in a segment required,
#'    segments short than this are rejected.  Ignore unless greater than 1.
#' @param resolution character, if \code{x} is missing, then try to load the
#'    \code{oce::coastlineWorld},\code{ocedata::coastlineWorldMedium} or 
#'    \code{ocedata::coastlineWorldFine} dataset 
#' @return a list of coastline segments suitable for casting to SpatialLines
split_coastline <- function(x, min_length = 3, 
   resolution = c("coarse", "medium", "fine")[1]){
   
   # convert an oce::coastline object to xy matrix
   coast_to_xy <- function(x){
      x <- as.matrix(cbind(x@data[['longitude']], x@data[['latitude']]))
      colnames(x) <- c("x", "y")
      rownames(x) <- seq_len(nrow(x))
      x
   }
   
   if (missing(x)) {
      x <- switch(tolower(resolution),
         'fine' =    {
                        stopifnot(require(ocedata))
                        data(coastlineWorldFine)
                        coast_to_xy(coastlineWorldFine)
                     },
         'medium'=   {
                        stopifnot(require(ocedata))
                        data(coastlineWorldMedium)
                        coast_to_xy(coastlineWorldMedium)
                     },
                     {
                        stopifnot(require(oce))
                        data(coastlineWorld)
                        coast_to_xy(coastlineWorld)
                     })
   }
   
   if (inherits(x, "coastline")) x <- coast_to_xy(x)
   stopifnot(is.matrix(x))
   stopifnot(all(colnames(x) %in% c("x","y"))) 
   
   r <- rle(is.na(x[,1]))
   ix <- r$values
   r$values <- 1:length(r$values)
   r$values[ix] <- 0
   z <- inverse.rle(r)
   xx <- lapply(split(rownames(x), z), function(i, x) {x[i,]}, x = x)
   xx[["0"]] <- NULL
   if (min_length > 1){
      len <- sapply(xx, nrow)
      xx[len < min_length] <- NULL
   }
   invisible(xx)
}

#' Convert coastline segments as matrices to SpatialLines
#' 
#' @export
#' @param x list of xy martices as per \code{split_coastline}
#'    if missing then the output of \code{split_coastline} is used
#' @param projection character string of projection defintion ala \code{CRS}
#' @param clip_to NULL or Spatial* object the coastline can be clipped to using
#'    \code{rgeos::gIntersection}
#' @param ... further arguments for split_coastline when x is not provided
#' @return a \code{SpatialLines} object
spatial_coastline <- function(x,
   projection = sp::CRS('+proj=longlat +datum=WGS84'),
   clip_to = NULL, ...){
   stopifnot(require(sp))
   if (missing(x)) x <- split_coastline(...)
   coast <- sp::Lines(lapply(x, function(x) sp::Line(x)), 'coast') 
   coast <- sp::SpatialLines(list(coast), proj4string = projection)
   if (!is.null(clip_to) && inherits(clip_to, 'Spatial')){
      stopifnot(require(rgeos))
      coast <- rgeos::gIntersection(clip_to, coast)
   }
   invisible(coast)
}