plotting stars proxy object with promises

[edzer]

I've been developing a number of examples here for working with stars proxy objects that have still unvaluated promises. A short reprex is this:

library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.9.0, GDAL 3.2.0, PROJ 7.2.0
tif = system.file("tif/L7_ETMs.tif", package = "stars")
r = split(read_stars(tif, proxy = TRUE))
pc = prcomp(as.data.frame(r)[,-(1:2)]) # based on all data
out = predict(r, pc)
plot(out, breaks = "equal")
plot(merge(out), breaks = "equal", join_zlim = FALSE)
library(tmap)
tm_shape(out) + tm_raster()
# stars_proxy object shown at 349 by 352 cells.
# Error in `[.data.frame`(as.data.frame(shp), shpnames) : 
#   undefined columns selected
# Calls: <Anonymous> ... print_tmap -> mapply -> <Anonymous> -> [ -> [.data.frame
# Execution halted
tm_shape(merge(out)) + tm_raster()
# stars_proxy object shown at 349 by 352 cells.
# Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
# Error in dim.stars(x) : length(d) == length(dim(x[[1]])) is not TRUE
# Calls: <Anonymous> ... do.call -> plot_1 -> nrow -> dim -> dim.stars -> stopifnot
# Execution halted

The problem is that out has different dimensions when it has promises then when after evaluating these:

> dim(out)
   x    y band 
 349  352    6 
> dim(st_as_stars(out))
  x   y 
349 352 

as well as a different number of attributes:

> length(out)
[1] 1
> length(st_as_stars(out))
[1] 6

This is all intentional: before evaluating, we don't know what'll come out. Calling st_as_stars realizes (evaluates) everything, and gives the right object. A challenge for tmap will be that you'll want to pick a good value for downsample in st_as_stars(), and for that you'd want to know the number of attributes / facets to be plotted? Maybe, for now, assume there is only one? Or give a warning that the user should better set something that controls downscale?

How to detect this? is.null(attr(out, "call_list")) until we have a better mechanism.

This looks complicated, but may be a killer feature of stars + tmap when rasters are huge; see the s2 example in the vignette. It is also something that ggplot() will unlikely be able to do automatically, IIUC (as it can't tell the plot size in pixels e.g. to geom_stars()).

pinging @HannaMeyer

[mtennekes]

Thx! Swapping two code chunks did the job, so now it works:-)

Picking a good value for downsample is indeed challenging, because the number of facets is also determined by whether band values are used as RGB values (making one plot instead of 3), and whether other shapes in the same plot call are faceted.

I think a good proxy would be to estimate the number of facets using the dimensions of the unevaluated proxy object. In most cases this would correspond to the realized number of facets. Also in your example, since in both forms (i.e. split to multiple attributes or merged to an additional dimension) result in the same number of facets.

An additional challenge: tm_shape(out) + tm_raster(c("PC2", "PC3")). In this case the variable values "PC2" and "PC3" are not known in out before evaluation. Probably, a good thing to do is to assume that the user didn't make an error, estimate the number of facets to be 2, and check after evaluation whether "PC2" and "PC3" are indeed attribute names.

Pinging #495 since I'm currently breaking my head how to implement data management and faceting in v4.

[przell]

Is this also affecting st_write() with proxy objects? I get Error in dim.stars(x) : length(d) == length(dim(x[[1]])) is not TRUE when I try to st_write a large proxy object (s2 tile). When I decrease the size it works.

[edzer]

@przell I assume you mean write_stars()? Could you provide a reprex, preferably with the S2 tile in package starsdata? (see here)