Add documentation for downsampling of rasters to speed up raster rendering

[paleolimbot]

Just extending discussion from #9 to be its own issue.

The idea is to delay any action involving the original raster object until there is a known height, width, and target resolution, so that large rasters can be passed to annotation_spatial.Raster() without fear of taking forever to render.

Question for @mdsumner: What is the series of events that usually takes place during this operation? I had been thinking something like:

target_res <- calculate_target_res(bounds, height, width, dpi)

st_make_grid(plot_bounds, cellsize = target_res) %>% 
  st_transform(raster_crs) %>% 
  raster::extract(raster, .) %>%
  make_a_raster_grob()

...but maybe it's a bad idea to reverse project a grid from the coord_sf() CRS to the raster CRS and resample that way? The other option would be to raster::projectRaster() %>% raster::resample()? Or maybe you've solved this in the hypertidy/lazyraster package?

Would also be nice to apply this to a rasterGrob() somehow, since these also take forever to render?

[mdsumner]

Definitely don't use make grid, there's no need. It is outlined in lazyraster but I can't expand on this atm. I'll flesh this out some time soon, but using raster functions will give enough to prototype, see what happens when raster applies the maxpixels argument, in various places. That needs to include efficient crop top, which needs work.

[paleolimbot]

Cool...I will set up the Stat/Geom/grob/makeContent() such that it works, and will await your sound guidance on how to make it better.

[mdsumner]

Great! Apologies I can't look now, very interested in this

[paleolimbot]

Whenever you're back, the framework is all set. Currently not faster than a regular rasterGrob(), but would probably be if it were resampling a much larger image and used with annotation_spatial(). Current approach just uses projectRaster() with an output template based on the extent of the panel.

library(ggspatial)
#> Loading required package: ggplot2
#> Loading required package: sf
#> Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3
load_longlake_data()

system.time(
  print(
    ggplot() +
      layer_spatial(longlake_osm, lazy = FALSE) +
      layer_spatial(longlake_depthdf)
  )
)

#>    user  system elapsed 
#>   2.899   0.308   3.430

system.time(
  print(
    ggplot() +
      layer_spatial(longlake_osm, lazy = TRUE) +
      layer_spatial(longlake_depthdf)
  )
)
#> 5.48175224540989 x 4.65488578335662 in (822.262836811483 x 698.232867503493 px) dpi: 150

#>    user  system elapsed 
#>   3.372   0.507   4.104

[paleolimbot]

@edzer ...this is what I was talking about last night. Implemented here: https://github.com/paleolimbot/ggspatial/blob/master/R/layer-spatial-raster.R#L246 , and here: https://github.com/paleolimbot/ggspatial/blob/master/R/layer-spatial-raster.R#L219

[paleolimbot]

Current problem: creating an appropriate destination template for stars::st_warp() (for my own future reference, or if @edzer has a solution that will make this work)

longlake_osm <- stars::read_stars(system.file("longlake/longlake.tif", package = "ggspatial"), proxy = TRUE)
panel_params <- list(
  x_range = c(2330056.03273654, 2332693.06989841), 
  y_range = c(229111.707417625, 231625.212947026)
)

new_dim <- c(x = 500, y = 400, band = 3)

template_raster <- stars::st_as_stars(
  dimensions = stars::st_dimensions(
    x = seq(panel_params$x_range[1], panel_params$x_range[2], length.out = new_dim[1]),
    y = seq(panel_params$y_range[1], panel_params$y_range[2], length.out = new_dim[2]),
    band = seq_len(new_dim[3])
  )
)

plot(longlake_osm)

warped <- stars::st_warp(longlake_osm, dest = template_raster, options = options)
#> Error in CPL_proj_direct(as.character(c(from[1], to[1])), as.matrix(pts)): no arguments in initialization list
plot(warped)
#> Error in plot(warped): object 'warped' not found

^{Created on 2019-01-19 by the reprex package (v0.2.1)}

[edzer]

Thanks! For rgb data you can do

plot(longlake_osm, rgb = 1:3, maxColorValue=255)

to plot an rgb image.

What you're trying doesn't work because template_raster doesn't have a crs; if I set it to the crs of longlake_osm the two do not overlap; that gives an error, which it shouldn't, but is probably not what you want. The following two work:

(stars::st_warp(longlake_osm, longlake_osm))
(warped <- stars::st_warp(longlake_osm, crs = st_crs(4326)))

Let me know (e.g. here or in a stars issue) if you run into things that you do want, but do not work.

[paleolimbot]

I'm back on this, and thanks to help by @edzer and @thomasp85 at the Tidyverse Dev Day, I'm a bit closer to having this done, but I am stuck on how to get a grob width in pixels at draw time. @thomasp85, if you have a moment, is there a way to do this? grDevices::dev.size("px") seems to always return a value that gives a dpi of 72, even though grDevices::dev.size("in") gives the correct size in inches.

I need an approximate DPI so I can downsample a potentially large raster (like is done in the stars package plot method). I could ask the user for this as a layer parameter, but it seems like there should be a way to query this?

This is a toy reprex that I think gets at what I'm trying to do:

new_custom_grob <- function() {
  grid::gTree(
    cl = "custom_grob"
  )
}

makeContext.custom_grob <- function(x) {

  width_in <- grid::convertWidth(grid::unit(1, "npc"), "in", valueOnly = TRUE)
  height_in <- grid::convertHeight(grid::unit(1, "npc"), "in", valueOnly = TRUE)
  
  # this DPI is not correct (`dev.size("px")`` returns incorrect value such that dpi is always 72)
  dpi <- grDevices::dev.size("px") / grDevices::dev.size("in")

  dpi_text <- sprintf(
    "%s x %s in: dpi: %sx%s", 
    width_in, height_in, dpi[1], dpi[2]
  )

  # add content to the gTree
  grid::setChildren(
    x,
    grid::gList(
      grid::textGrob(dpi_text)
    )
  )
}

grid.custom_grob <- function() {
  grid::grid.newpage()
  grid::grid.draw(new_custom_grob())
}

withr::with_png("test1.png", grid.custom_grob(), height = 300, width = 400, res = 100)
grid::grid.raster(png::readPNG("test1.png"))

withr::with_png("test2.png", grid.custom_grob(), height = 300, width = 400, res = 200)
grid::grid.newpage()
grid::grid.raster(png::readPNG("test2.png"))

^{Created on 2019-01-29 by the reprex package (v0.2.1)}

[paleolimbot]

It looks like the makeContext() method gets called 3 times in the RStudio graphics device: once during the "printing", then twice. The time during the printing doesn't have the same size as the other two times, which are called from the graphics device. Other devices appear only to call makeContext() once. Caching the plot reduces the number of times the raster has to get resampled to 2 (in RStudio)...this is possible by specifying an element of the gTree as an environment.

Log output from some sleuthing:

***Checking for cached raster in x$cache[['raster_925x786']]
vapply(rlang::ctxt_stack(), function(x) as.character(x)[3], character(1))
makeContent.geom_spatial_raster_lazy(x)
makeContent(x)
drawGTree(x)
recordGraphics(drawGTree(x), list(x = x), getNamespace("grid"))
grid.draw.gTree(x$children[[i]], recording = FALSE)
grid.draw(x$children[[i]], recording = FALSE)
drawGTree(x)
recordGraphics(drawGTree(x), list(x = x), getNamespace("grid"))
grid.draw.gTree(x$children[[i]], recording = FALSE)
grid.draw(x$children[[i]], recording = FALSE)
drawGTree(x)
recordGraphics(drawGTree(x), list(x = x), getNamespace("grid"))
grid.draw.gTree(gtable)
grid.draw(gtable)
print.ggplot(ggplot() + layer_spatial(longlake_osm, lazy = TRUE) + layer_spatial(longlake_depthdf) + labs(caption = "Should have a little grey area around the sides, roughly N-S projection"))
print(ggplot() + layer_spatial(longlake_osm, lazy = TRUE) + layer_spatial(longlake_depthdf) + labs(caption = "Should have a little grey area around the sides, roughly N-S projection"))
print(ggplot() + layer_spatial(longlake_osm, lazy = TRUE) + layer_spatial(longlake_depthdf) + labs(caption = "Should have a little grey area around the sides, roughly N-S projection"))
NULL
Calling raster_grob_from_raster()

***Checking for cached raster in x$cache[['raster_670x569']]
vapply(rlang::ctxt_stack(), function(x) as.character(x)[3], character(1))
makeContent.geom_spatial_raster_lazy(x)
makeContent(x)
drawGTree(x)
recordGraphics(drawGTree(x), list(x = x), getNamespace("grid"))
grid.draw.gTree(x$children[[i]], recording = FALSE)
grid.draw(x$children[[i]], recording = FALSE)
drawGTree(x)
recordGraphics(drawGTree(x), list(x = x), getNamespace("grid"))
grid.draw.gTree(x$children[[i]], recording = FALSE)
grid.draw(x$children[[i]], recording = FALSE)
drawGTree(x)
doTryCatch(return(expr), name, parentenv, handler)
tryCatchOne(expr, names, parentenv, handlers[[1]])
tryCatchList(expr, classes, parentenv, handlers)
tryCatch(.Call("rs_GEcopyDisplayList", fromDevice), error = function(e) warning(e))
(function (fromDevice) 
{
    tryCatch(.Call("rs_GEcopyDisplayList", fromDevice), error = function(e) warning(e))
})(1)
NULL
Calling raster_grob_from_raster()

***Checking for cached raster in x$cache[['raster_670x569']]
vapply(rlang::ctxt_stack(), function(x) as.character(x)[3], character(1))
makeContent.geom_spatial_raster_lazy(x)
makeContent(x)
drawGTree(x)
recordGraphics(drawGTree(x), list(x = x), getNamespace("grid"))
grid.draw.gTree(x$children[[i]], recording = FALSE)
grid.draw(x$children[[i]], recording = FALSE)
drawGTree(x)
recordGraphics(drawGTree(x), list(x = x), getNamespace("grid"))
grid.draw.gTree(x$children[[i]], recording = FALSE)
grid.draw(x$children[[i]], recording = FALSE)
drawGTree(x)
doTryCatch(return(expr), name, parentenv, handler)
tryCatchOne(expr, names, parentenv, handlers[[1]])
tryCatchList(expr, classes, parentenv, handlers)
tryCatch(.Call("rs_GEcopyDisplayList", fromDevice), error = function(e) warning(e))
(function (fromDevice) 
{
    tryCatch(.Call("rs_GEcopyDisplayList", fromDevice), error = function(e) warning(e))
})(1)
NULL
Using cached raster

[paleolimbot]

Cooking with fire! It turns out it is very important when resampling not to create a rasterGrob() with more pixels than is needed, since grid's rendering of the rasterGrob takes a bit.

This involves capping the maximum number of pixels in the x and y direction in the final projection space before draw time. This also makes more effective use of the draw cache, since the same size rasterGrob might be relevant as the plot size increases.

The draw-time resampling (for rgb, raster package) is now faster than the build-time resampling!

library(ggspatial)
#> Loading required package: ggplot2
load_longlake_data(which = c("longlake_osm", "longlake_depthdf"), raster_format = "raster")

system.time(
  print(
    ggplot() +
      layer_spatial(longlake_osm, lazy = FALSE) +
      layer_spatial(longlake_depthdf)
  )
)

#>    user  system elapsed 
#>   1.914   0.328   2.477

system.time(
  print(
    ggplot() +
      layer_spatial(longlake_osm, lazy = TRUE) +
      layer_spatial(longlake_depthdf)
  )
)

#>    user  system elapsed 
#>   1.717   0.301   2.250

^{Created on 2019-02-03 by the reprex package (v0.2.1)}