oce_meshplot.R

library(raster)

## Loading required package: sp

plot_raster <- function(x, crs = NULL, add = FALSE, ..., zlim = NULL,  col = NULL) {
  UseMethod("plot_raster")
}

plot_raster.topo <- function(x, crs = NULL, add = FALSE, ..., zlim = NULL, col = NULL)  {
  x <- raster(list(x = x@data$longitude, y = x@data$latitude, z = x@data$z),
              crs = "+proj=longlat +datum=WGS84")
  plot_raster(x, crs = crs, add = add, ..., col = col)
}
plot_raster.BasicRaster <- function(x, crs = NULL, add = FALSE, ..., zlim = NULL,   col = NULL)  {
  if (raster::nlayers(x) > 1) warning("only first layer plotted")
  qm <- quadmesh::quadmesh(x[[1]])
  xy <- t(qm$vb[1:2, ])
  if (!is.null(crs)) {

    src <- raster::projection(x)
    if (!is.na(src) && !identical(src, crs)) {
      if (raster::isLonLat(src)) xy <- xy * pi/180
      xy <- proj4::ptransform(xy, src, crs)
    }
  }
  xy <- as.data.frame(xy)
  names(xy)[1:2] <- c("x", "y")
  quadindex <- qm$ib
  ## propagate NA into the quads
  bad <- !is.finite(xy$x) | !is.finite(xy$y)
  idx <- which(bad)
  bad2 <- quadindex[1,] %in% idx | quadindex[2,] %in% idx |
    quadindex[3,] %in% idx | quadindex[4,] %in% idx
  xy <- xy[!bad, ]
  quadindex <- quadindex[,!bad2]
  ylim <- range(xy$y)
  xlim <- range(xy$x)
  if (!add) {
    grid::grid.newpage()
    asp <- if (raster::isLonLat(crs)) 1/cos(mean(ylim) * pi/180) else 1
    plot.new()
    plot.window(xlim = xlim, ylim = ylim, asp = asp)
  }
  if (is.null(col)) col <- viridis::viridis(64)
  v <- values(x[[1]])
  if (is.null(zlim)) zlim <- range(v, na.rm = TRUE)
  cols <- col[scales::rescale(v, to = c(1, length(col)), from = zlim)]
  sq <- seq_len(ncol(quadindex))[!is.na(cols)]
  # for (i in sq) {
  #  polypath(xy[quadindex[, i], ], col = cols[i], border = NA)
  #}
  vp <- gridBase::baseViewports()
  grid::pushViewport(vp$inner, vp$figure, vp$plot)
  grid::grid.polygon(xy[quadindex, 1], xy[quadindex, 2],
                     rep(seq_len(ncol(quadindex)), each = 4),
                     gp = grid::gpar(fill = cols, col = NA),
                     default.units = "native")
  grid::popViewport()
  invisible(NULL)
}

library(oce)

## Loading required package: gsw

## Loading required package: testthat

data("topoWorld")
data("coastlineWorld")
topo <- topoWorld
lonlim <- c(-70, -50)
latlim <- c(40, 50)
topo <- subset(topo, latlim[1] < latitude & latitude < latlim[2])
topo <- subset(topo, lonlim[1] < longitude & longitude < lonlim[2])


proj <- "+proj=lcc +lat_1=40 +lat_2=50 +lon_0=-60"

## this is very slow to initialize the plot?
mapPlot(lonlim, latlim, projection= proj, type = "n")
plot_raster(topo, crs = proj,
            col = oce.colorsGebco(64), add = TRUE)
mapLines(coastlineWorld)

## can we level-up?

(topo <- aggregate(raadtools::readtopo("gebco_14",
                          xylim = extent(lonlim, latlim)), fact = 4))

## class       : RasterLayer 
## dimensions  : 300, 600, 180000  (nrow, ncol, ncell)
## resolution  : 0.03333333, 0.03333333  (x, y)
## extent      : -70, -50, 40, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : Elevation.relative.to.sea.level 
## values      : -5621.812, 1170  (min, max)

system.time(mapPlot(lonlim, latlim, projection= proj, type = "n"))

##    user  system elapsed 
##  38.934   0.107  39.041

system.time(plot_raster(topo, crs = proj,
            col = oce.colorsGebco(64), add = TRUE))

##    user  system elapsed 
##   3.308   0.000   3.308

mapLines(coastlineWorld)