layer_spatial.stars produces much lower resolution output than layer_spatial.Raster on same data

[MilesMcBain]

I am struggling to understand why the output I get for {stars} RGB rasters with layer_spatial looks so grainy when compared with {raster} RGB rasters? Is there some kind of interpolation happening?

An example follows. Note the difference in image quality produced by ggsave using the same parameters on raster vs stars.

Code in the as_stars_tile function here comes from @edzer in :https://github.com/r-spatial/stars/issues/315. It depends on the dev version of stars.

library(snapbox) ## remotes::install_github("anthonynorth/snapbox")
library(sf)
library(stars)
library(ggspatial)
library(ggplot2)

bbox <- st_bbox(
   c(xmin = 147, ymin = -43, xmax = 147.7, ymax = -42.65),
   crs = 4326
)

map_img <- snapbox:::get_map_image(
    bbox = bbox,
    map_style = mapbox_dark(),
    mapbox_api_access_token = Sys.getenv("MAPBOX_ACCESS_TOKEN"),
    width = 1280,
    height = 873,
    retina = TRUE,
    mapbox_logo = TRUE,
    attribution = TRUE,
    purge_cache = TRUE
  )

as_stars_tile <- function(map_img, map_img_bbox) {
  stars_img_transposed <-
    stars::st_as_stars(.x = aperm(map_img, c(2, 1, 3))[, nrow(map_img):1, ])
  img_bbox <-
    sf::st_bbox(map_img_bbox[c("xmin", "xmax", "ymin", "ymax")],
    crs = 3857
  )
  stars_tile <- stars::st_set_bbox(stars_img_transposed, img_bbox)
  names(attr(stars_tile, "dimensions"))[3] <- "band"
  stars_tile
}

as_raster_tile <- function(map_img, map_img_bbox) {
raster_tile <-
  raster::brick(map_img) %>%
  raster::setExtent(raster::extent(
    map_img_bbox$xmin,
    map_img_bbox$xmax,
    map_img_bbox$ymin,
    map_img_bbox$ymax
  ))
raster::crs(raster_tile) <- "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
raster_tile
}

stars_tile  <- as_stars_tile(map_img, snapbox:::get_mercator_bbox.bbox(bbox))

raster_tile <- as_raster_tile(map_img, snapbox:::get_mercator_bbox.bbox(bbox))

p1 <-
ggplot() + 
  layer_spatial(stars_tile)

p2 <-
ggplot() + 
  layer_spatial(raster_tile)

ggsave("stars.png", p1,  width = 12, height = 8, units = "in")
ggsave("raster.png", p2,  width = 12, height = 8, units = "in")

stars.png

raster.png

[MilesMcBain]

Simpler reprex with the long lake osm data:

  library(ggplot2)
  library(stars)
  library(ggspatial)
  library(raster)

  stars_rgb <- stars::read_stars(system.file("longlake/longlake.tif", package = "ggspatial"))
  ggspatial:::load_longlake_data(which = c("longlake_osm"), raster_format = "raster")

  ggplot() + 
    layer_spatial(longlake_osm)

  ggsave("raster.png", width = 12, height = 8, units = "in")

  ggplot() +
    layer_spatial(stars_rgb)
  ggsave("stars.png",  width = 12, height = 8, units = "in")

raster.png

stars.png

[paleolimbot]

Thanks for reprexing this so nicely! I'll do a more thorough investigation in a few, but I think if you set interpolate = TRUE for both, you should get the same result. I do some guessing to try to figure out whether or not to use interpolation, and it looks like my guessing is incorrect somewhere!

[MilesMcBain]

Setting interpolate = TRUE makes a marginal improvement but so marginal that at first I thought it was having no effect.

I think it's because st_warp uses an inferior interpolation method by default, compared with the bilinear interpolation that gets used in the raster case.

I tried to turn on bilinear interpolation in the st_warp call in raster_grob_from_stars however the template_raster that is created is not compatible with the dest argument expected by st_warp for this case. It seems that it needs to have the data, so that dest[[1]] can work.

[paleolimbot]

Ooof. Thanks for putting the time in to diagnose this...I'll fix in the next few days!

[paleolimbot]

I think this is fixed...I couldn't get stars::st_warp() to work for this, so I went straight to sf::gdal_util("warp", ...). Let me know if this fixes the problem (or creates new ones...)! I need to improve my testing for stars objects...this issue reminded me that they need serious updating.

[MilesMcBain]

This is a vast improvement. Interestingly to my eyes raster still somehow has the edge. See the latest stars.png:

It looks like raster is hand rolling it's own bilinear interpolation. The png file size it produces is slightly larger.

But I think this is good enough. Thankyou!

[edzer]

Maybe I'm naive, but why would you do a bilinear interpolation on an image with a categorical variable?

[MilesMcBain]

That's not what's happening here. This long lake data is just an example.

The real use is for map images that have text, hillshading, roads, glyphs etc. Could also have satellite images as base layer.

[MilesMcBain]

But the comment does make me think that the interpolate arg may have a confusing name.

Interpolation always occurs to match the raster size to the plot size, however when interpolate = TRUE bilinear interpolation is done instead of nereast neighbour. So maybe you're better off having an 'interpolation_method` arg? 🤔

[paleolimbot]

@edzer You don't! The default for interpolate is FALSE except for RGB (although my detection of RGB isn't perfect). I changed this after you brought it up in a previous thread (thanks!).

I'll investigate the lower quality and add some tests for stars objects...the code isn't quite as I remember it and likely has a few other issues.

interpolate = TRUE/FALSE is was named from the current geom_raster(), which I think takes the name from rasterGrob(). Really, there's two interpolations happening...one if/when the raster has to be reprojected (operates on values) and one when the colours get resampled to match the screen resolution (operates on colours). Ideally there would only be one (each screen pixel would get resampled + coloured once for each window resize), but when I tried this a few years ago I ran into some limitations. Ideally I'd like to solve those issues together...I think interpolate = TRUE or interpolate = FALSE does the job most of the time.