mapPlot() should handle coastlines in SFG polygon format

[dankelley]

I'm realizing that this format is pretty handy. An example is in the graph below ... I won't bother with the details, but this relates to my recent edge-of-world work (mentioned in some recent issues). The code stems from private (teaching) material at my private repo PS. the code to make the image in first comment is at private location /Users/kelley/git/teaching_shiny/ortho_test

In the graphs, the tan coastline in the RH panel is something I created by reading the existing data(coastlineWorld) file, and then creating a multi-polygon object.

I'm not suggesting (yet) that we change the structure of the built-in file, but if users are doing the "right" thing and extracting with e.g. coastline[["longitude"]] (as opposed to coastline@data$longitude) then changing the file would not affect user code, because I can easily change what [[ does. I may also add an [[ operation, maybe like coastlineWorld[["polygons"]], to return polygons ... it only takes 0.05s for this resolution.

> class(coastlineWorldPolygon)
[1] "XY"           "MULTIPOLYGON" "sfg"   

Code

library(oce)
library(shiny)
data(coastlineWorld)

dmsg <- function(...) if (debug > 0L) cat(...)

createCoastlinePolygon <- function(coastline)
{
    lon <- c(coastline[["longitude"]], NA)
    lat <- c(coastline[["latitude"]], NA)
    NAs <- which(is.na(lon))
    if (any(1L == diff(NAs))) {
        stop("malformed coastline file (adjacent NA values)")
    }
    npoly <- length(NAs) - 1L
    p <- vector("list", npoly)
    for (i in seq_len(npoly)) {
        look <- seq(NAs[i] + 1L, NAs[i+1] - 1L)
        LON <- lon[look]
        LAT <- lat[look]
        len <- length(look)
        # Guess how to close those polygons that are not closed
        if (LON[1] != LON[len] || LAT[1] != LAT[len]) {
            LON <- c(LON, LON[1])
            LAT <- c(LAT, LAT[1])
        }
        p[[i]] <- list(cbind(LON, LAT))
    }
    sf::st_multipolygon(p)
}

# The first call of the next line takes 0.6s but if I do it again, it takes
# 0.007s.  I am not sure why, exactly, but memory allocation seems like a reason for such a thing.

system.time(coastlineWorldPolygon <- createCoastlinePolygon(coastlineWorld))

Code to go back to vectors

    cwp <- coastlineWorldPolygon # for speed of interactive tests
    class(cwp) # XY MULTIPOLYGON sfg
    methods(class = "sfg")
    system.time({ # 0.012 0.006 0.019 for coastlineWorld (10766 points)
        alon <- NULL
        alat <- NULL
        for (p in cwp) {
            pm <- as.matrix(p[[1]])
            alon <- c(alon, pm[, 1], NA)
            alat <- c(alat, pm[, 2], NA)
        }
    })