meshplot-example.R

scl <- function(x) {
  rg <- range(x, na.rm = TRUE);
  (x - rg[1])/diff(rg)
}

meshplot_discrete <- function(x, crs = NULL, ...) {
  UseMethod("meshplot_discrete")
}
meshplot_discrete.BasicRaster <- function(x, crs = NULL, ...) {
  print("converting to single RasterLayer")
  meshplot_discrete(x[[1]])
}
meshplot_discrete.RasterLayer <- function(x, crs = NULL, ...) {
  qm <- quadmesh::quadmesh(x, na.rm = FALSE)
  ib <- qm$ib
  xy <- t(qm$vb[1:2, ])
  rm(qm)
  if (!is.null(crs) ) {
    if (!raster::isLonLat(x)) {
      xy <- rgdal::project(xy, raster::projection(x), inv = TRUE)
    }
    if (!raster::isLonLat(crs)) xy <- rgdal::project(xy, crs)
  }
  ## we have to remove any infinite vertices
  ## as this affects the entire thing
  bad <- !is.finite(xy[,1]) | !is.finite(xy[,2])
  ## but we must identify the bad xy in the index
  if (any(bad)) ib <- ib[,-which(bad)]
  ## the scale function must not propagate NA
  xx <- scl(xy[c(ib),1])
  yy <- scl(xy[c(ib),2])
  ## we need a identifier grouping for each 4-vertex polygon
  id <- rep(seq_len(ncol(ib)), each  = nrow(ib))
  grid::grid.newpage()

  ## we also have to deal with any values that are NA
  ## because they propagate to destroy the id
  cols <- viridis::viridis(100)[scl(values(x)) * 99 + 1]
  if (any(is.na(cols))) {
    colsna <- rep(cols, each = nrow(ib))
    bad2 <- is.na(colsna)
    xx <- xx[!bad2]
    yy <- yy[!bad2]
    id <- id[!bad2]
    cols <- cols[!is.na(cols)]
  }
  ## turn off col to avoid plotting pixel boundaries
  grid::grid.polygon(xx, yy, id,
                     gp = grid::gpar(col = NA, fill = cols))

}

meshplot_discrete.Spatial <- function(x, crs = NULL,  ...) {
  #if (!is.null(field)) vals <- x[[field]] else vals <- seq_len(nrow(x))
  x <- rangl::rangl(x)

  tt <- rangl:::th3d()

  tt$vb <- t(cbind(x$v$x_, x$v$y_, 0, 1))

  vv <- x$v[, "vertex_"]
  vv$row_n <- seq(nrow(vv))
  pindex <- dplyr::inner_join(dplyr::inner_join(x$o[, c("object_")], x$t), x$tXv)
  vindex <- dplyr::inner_join(x$tXv, vv, "vertex_")
  tt$it <- t(matrix(vindex$row_n, ncol = 3, byrow = TRUE))

  it <- tt$it
  xy <- t(tt$vb[1:2, ])
  rm(tt)
  if (!is.null(crs) ) {
    if (!raster::isLonLat(x)) {
      xy <- rgdal::project(xy, raster::projection(x), inv = TRUE)
    }
    if (!raster::isLonLat(crs)) xy <- rgdal::project(xy, crs)
  }
  ## we have to remove any infinite vertices
  ## as this affects the entire thing
  bad <- !is.finite(xy[,1]) | !is.finite(xy[,2])
  ## but we must identify the bad xy in the index
  if (any(bad)) it <- it[,-which(bad)]
  ## the scale function must not propagate NA
  xx <- scl(xy[c(it),1])
  yy <- scl(xy[c(it),2])
  ## we need a identifier grouping for each 4-vertex polygon
  id <- rep(seq_len(ncol(it)), each  = nrow(it))
  colid <- seq_len(ncol(it))
  grid::grid.newpage()

  ## we also have to deal with any values that are NA
  ## because they propagate to destroy the id
  cols <- viridis::viridis(100)[scl(colid) * 99 + 1]
  grid::grid.polygon(xx, yy, id,
                     gp = grid::gpar(col = NA, fill = cols))

}


library(raster)

## Loading required package: sp

data("wrld_simpl", package = "maptools")
r <- rasterize(wrld_simpl, raster(wrld_simpl, ncol = 720, nrow = 360))
## show the utility and speed of a mesh plot for near-lossless transformation
system.time(meshplot_discrete(r))

##    user  system elapsed 
##   3.375   0.427   3.802

system.time(meshplot_discrete(r, crs = "+proj=laea +datum=WGS84"))

##    user  system elapsed 
##   2.581   0.187   2.768

system.time(meshplot_discrete(wrld_simpl, crs = "+proj=laea +datum=WGS84"))

## Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11 display

## Warning: 'rgl_init' failed, running with rgl.useNULL = TRUE

## Joining, by = "object_"

## Joining, by = "triangle_"

##    user  system elapsed 
##  22.698   0.330  23.072