NOmap: Northern Ocean maps

#' Generate a ggplot2 representation of an SOmap object
#'
#' Note: this function is still experimental! Use at your own risk.
#' 
#' @param x SOmap: object as returned by \code{SOmap} or \code{SOmap2}
#'
#' @return A ggplot object
#'
#' @examples
#' \dontrun{
#'   p <- SOmap2(Trim = -45, IWC = TRUE, IWClab = TRUE, Grats = TRUE, fronts = TRUE,
#'               MPA = TRUE, MPAlab = TRUE)
#'   SOgg(p)
#' }
#'
#' @export
SOgg <- function(x) {
    if (!inherits(x, "SOmap")) stop("x must be an SOmap object as generated by SOmap or SOmap2")
    myext <- extent(x$target)
    ## plot bathymetry
    bdf <- raster::as.data.frame(raster::trim(SOmap::latmask(x$bathy$plotargs$x, latitude = x$trim)), xy = TRUE)
    names(bdf)[3] <- "Depth"
    p <- ggplot(data = bdf, aes_string(x = "x", y = "y")) + geom_raster(aes_string(fill = "Depth"))
    p <- p + coord_sf(default = TRUE) ## default = TRUE makes this the default coordinate system for geom_sf objects
    if (!is.null(x$bathy_legend)) {
        p <- p + scale_fill_gradientn(colours = x$bathy$plotargs$col, na.value = "#FFFFFF00")
##        raster::plot(x$bathy_legend$mask$graticule, border = x$bathy_legend$mask$border, col = x$bathy_legend$mask$col, add = TRUE) ## white mask
##        raster::plot(x$bathy_legend$ticks$ticks, add = TRUE, col = x$bathy_legend$ticks$col)
##        raster::plot(x$bathy_legend$legend$legend, lwd = x$bathy_legend$legend$lwd, add = TRUE)
##        raster::plot(x$bathy_legend$legend$legend, border = x$bathy_legend$legend$border, col = x$bathy_legend$legend$col, add = TRUE)
##        suppressMessages(this <- fortify(x$bathy_legend$legend$legend))
##        this$fill <- x$bathy_legend$legend$col[as.numeric(this$id)]
  ## nb, would have to draw each of these polygons in turn, with a fixed (non-aesthetic) fill specification? otherwise it will interfere with the scale_fill_gradientn
##        raster::plot(x$bathy_legend$graticules$graticules, border = x$bathy_legend$graticules$border, col = x$bathy_legend$graticules$col, add = TRUE)
##        text(x$bathy_legend$labels$data, labels = x$bathy_legend$labels$labels, cex = x$bathy_legend$labels$cex, adj = x$bathy_legend$labels$adj)
    } else {
        p <- p + scale_fill_gradientn(colours = x$bathy$plotargs$col, na.value = "#FFFFFF00", guide = NULL)
    }

    ## buffer to use for cropping things back to our extent of interest
    buf <- make_buf(x$trim+2, x$projection)

    if (!is.null(x$coastline)) {
        ## the coastline data has to be trimmed to our northernmost latitude
        ## masking (using e.g. x$bathy_legend$mask$graticule) is likely to be problematic because of z-ordering
        ## TODO check that this trimming is robust
        this <- suppressWarnings(sf::st_intersection(buf, x$coastline$plotargs$x))
        p <- p + geom_sf(data = this, fill = x$coastline$plotargs$col, col = x$coastline$plotargs$border, inherit.aes = FALSE)
    }
#    if (!is.null(x$iwc)) {
#        for (ii in seq_len(length(x$iwc$data))) {
#            this <- as.data.frame(x$iwc$data[[ii]])
#            names(this) <- c("x", "y")
#            p <- p + geom_path(data = this, col = x$iwc$col)
#        }
#        if (!is.null(x$iwc$labels)) {
#            this <- as.data.frame(x$iwc$labels$data)
#            p <- p + geom_text(data = this, aes_string(x = "lon", y = "lat", label = "a"), col = x$iwc$labels$col)##, cex = x$iwc$labels$cex, pos = x$iwc$labels$pos, offset = x$iwc$labels$offset)
#        }
#    }

    ## fronts
    if (!is.null(x$fronts)) {
        this <- x$fronts$plotargs$x##sf::st_as_sf(x$fronts$plotargs$x)
        this <- suppressWarnings(sf::st_intersection(buf, this))
        thiscol <- rep(x$fronts$plotargs$col, ceiling(nrow(this)/length(x$fronts$plotargs$col)))
        thiscol <- thiscol[seq_len(nrow(this))]
        p <- p + geom_sf(data = this, col = thiscol, inherit.aes = FALSE)
    }
    ## Graticule grid
    if (!is.null(x$graticule)) {
        this <- fortify(x$graticule$plotargs$x)
        this <- this[this$long >= myext[1] & this$long <= myext[2] & this$lat >= myext[3] & this$lat <= myext[4], ]
        ##raster::plot(x$graticule[[1]]$data, add = TRUE, col = x$graticule[[1]]$col, lty = x$graticule[[1]]$lty)
        p <- p + geom_path(data = this, aes_string(x = "long", y = "lat", group = "group"), col = x$graticule$plotargs$col, linetype = x$graticule$plotargs$lty)
        if (!is.null(x$graticule$labels)) {
            this <- as.data.frame(x$graticule$labels$plotargs$x)
            this <- this[this$x >= myext[1] & this$x <= myext[2] & this$y >= myext[3] & this$y <= myext[4], ]
            p <- p + geom_text(data = this, aes_string(label = "lab"), parse = TRUE, col = x$graticule$labels$plotargs$col)##, cex = x$graticule$labels$plotargs$cex)
            ## can't use cex on that, it gets translated to a size aesthetic, which is an absolute not relative size
        }
    }

    p <- p + labs() + theme(axis.title = element_blank(),
                            axis.text.x = element_blank(), axis.ticks.x = element_blank(),
                            axis.text.y = element_blank(), axis.ticks.y = element_blank(),
                            panel.border = element_blank(),
                            panel.background = element_blank())

##    plot_research_blocks(x$research_blocks)
##    plot_sprfmo(x$sprfmo_research_blocks)
##    plot_ssru(x$ccamlr_ssru)
##    plot_ssmu(x$ccamlr_ssmu)
##    plot_ccamlr_areas(x$ccamlr_statistical_areas)

    ## NOTE, parse args to geom_text and geom_sf_text seem fragile, need better user control
    if (!is.null(x$eez)) {
        this <- suppressWarnings(sf::st_intersection(buf, sf::st_as_sf(x$eez$plotargs$x)))
        p <- p + geom_sf(data = this, col = x$eez$plotargs$border, fill = x$eez$plotargs$col, inherit.aes = FALSE)
        if (!is.null(x$eez$labels)) {
            this <- x$eez$labels$plotargs$x
            this$lab <- x$eez$labels$plotargs$labels
            this <- suppressWarnings(sf::st_intersection(buf, sf::st_as_sf(this)))
            p <- p + geom_sf_text(data = this, aes_string(label = "lab"), parse = FALSE, col = x$eez$labels$plotargs$col, inherit.aes = FALSE)##, cex = x$eez$labels$cex, pos = x$eez$labels$pos, offset = x$eez$labels$offset)
        }
    }

    if (!is.null(x$mpa)) {
        this <- suppressWarnings(sf::st_intersection(buf, sf::st_as_sf(x$mpa$plotargs$x)))
        p <- p + geom_sf(data = this, col = x$mpa$plotargs$border, fill = x$mpa$plotargs$col, inherit.aes = FALSE)
        if (!is.null(x$mpa$labels)) {
            this <- x$mpa$labels$plotargs$x
            this$lab <- x$mpa$labels$plotargs$labels
            this <- suppressWarnings(sf::st_intersection(buf, sf::st_as_sf(this)))
            p <- p + geom_sf_text(data = this, aes_string(label = "lab"), parse = TRUE, col = x$mpa$labels$plotargs$col, inherit.aes = FALSE)##, cex = x$mpa$labels$cex, pos = x$mpa$labels$pos, offset = x$mpa$labels$offset)
        }
    }

    if (!is.null(x$ccamlr_planning_domains)) {
        this <- suppressWarnings(sf::st_intersection(buf, sf::st_as_sf(x$ccamlr_planning_domains$plotargs$x)))
        ## TODO fix that intersection, is slow because of complexity of coastline
        p <- p + geom_sf(data = this, col = x$ccamlr_planning_domains$plotargs$border, fill = x$ccamlr_planning_domains$plotargs$col, inherit.aes = FALSE)
        if (!is.null(x$ccamlr_planning_domains$labels)) {
            ## this is horrible code
            crds <- as.data.frame(sp::coordinates(x$ccamlr_planning_domains$labels[[1]]$plotargs$x))
            names(crds) <- c("x", "y")
            for (ii in seq_len(length(x$ccamlr_planning_domains$labels))) {
    #            this <- x$ccamlr_planning_domains$labels[[ii]]$plotargs$x
    #            this$lab <- x$ccamlr_planning_domains$labels[[ii]]$plotargs$labels
                crds$lab <- x$ccamlr_planning_domains$labels[[ii]]$plotargs$labels
    #            this <- suppressWarnings(sf::st_intersection(buf, sf::st_as_sf(this)))
                p <- p + geom_text(data = crds, aes_string(x = "x", y = "y", label = "lab"), parse = FALSE, col = x$ccamlr_planning_domains$labels[[ii]]$plotargs$col, inherit.aes = FALSE)##, cex = x$ccamlr_planning_domains$labels[[ii]]$cex, pos = x$ccamlr_planning_domains$labels[[ii]]$pos, offset = x$ccamlr_planning_domains$labels[[ii]]$offset)
            }
        }
    }

    if (!is.null(x$border)) {
        suppressMessages(this <- fortify(x$border$plotargs$x))
        this$col <- (as.numeric(this$id) %% length(x$border$plotargs$col)) + 1
        for (ii in seq_along(x$border$plotargs$col)) {
            p <- p + geom_polygon(data = this[this$col == ii, ], aes_string(x = "long", y = "lat", group = "group"), fill = x$border$plotargs$col[ii], col = "black")
        }
    }
    p
}