nara_gdalraster.md

  library(gdalraster)
#> GDAL 3.13.0dev-0e4f0cadf6 (released 2026-01-20), GEOS 3.12.1, PROJ 9.7.0
dsn <- "WMTS:https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml,layer=World_Imagery"

ds <- new(GDALRaster, dsn)
d <- read_ds(ds, out_xsize = 1024, out_ysize = 1024, bands = 1:3)
bb <- ds$bbox()
plot(NA, xlim = bb[c(1, 3)], ylim = bb[c(2, 4)], asp = 1, xlab = "", ylab = "")
a <- aperm(array(d, c(1024, 1024, 3)), c(2,1 ,3))
a01 <- a/255
system.time({
rasterImage(a01, bb[1], bb[2], bb[3], bb[4])
})
#>    user  system elapsed 
#>   0.212   0.024   0.238

native_rgb <- function(r, g, b, dm) {
  #res <- bitwOr(bitwOr(bitwOr(r, bitwShiftL(g, 8L)), bitwShiftL(b, 16L)), -16777216L)
  res <- as.integer(r) + as.integer(g) * 256L + as.integer(b) * 65536L - 16777216L
  res <- matrix(res, dm[2L], dm[1L], byrow = TRUE)
  class(res) <- "nativeRaster"
  attr(res, "channels") <- 3L
  res
}

b <- native_rgb(a[,,1L], a[,,2L], a[,,3L], c(1024, 1024))
system.time({
  rasterImage(b, bb[1], bb[2], bb[3], bb[4])
})

#>    user  system elapsed 
#>   0.001   0.000   0.001

^{Created on 2026-01-22 with reprex v2.0.2}

my only thought there (and I had a few half baked responses but had to try a bit of code) is liberal use of -scale and -ot - in just nsidc you immediately see the need for masking values or applying an existing palette, but this gets you a long way for quick and dirty vis. I think it's worth exploring how far we can take it, there's obvious cases of existing Gray/RGB/A Byte, palette, and scaling based on data-range or type-range, and as a minimum I would want application of a colour map like image() provides (cols, or cols + breaks). The actual building of a colour map is an almost orthogonal task I think, but if that can slot in it would be really neat. You absolutely need quantile-like options for "RGB" from 16 bit imagery, even as baked in defaults and I know plot_raster has quite a bit on colour options already.

Does that answer?

## Int16
gebco <- "/vsicurl/https://projects.pawsey.org.au/idea-gebco-tif/GEBCO_2024.tif"
## UInt16 palette image
nsidc <- "/vsicurl/https://noaadata.apps.nsidc.org/NOAA/G02135/south/daily/geotiff/2026/01_Jan/S_20260116_concentration_v4.0.tif"

fit_dims <- function(size = 1024L, wh = c(size, size)) {
  w <- wh[1L]; h <- wh[2L];
  as.integer(ceiling(size * c(w, h) / max(w, h)))
}
read_to_nr <- function(dsn, size = 1024L) {
  ds <- new(gdalraster::GDALRaster, dsn)
  dm <- fit_dims(size, ds$dim()[1:2])
  on.exit(ds$close(), add = TRUE)
  read_to_nativeRaster(ds, out_xsize = dm[1L], out_ysize = dm[2L])
}
make_byte <- function(x, scale = FALSE) {
  x <- sprintf("vrt://%s?ot=Byte", x)
  if (scale) {
    x <- sprintf("%s&scale=true", x)
  }
  x
}

ximage::ximage(read_to_nr(make_byte(gebco)))
ximage::ximage(read_to_nr(make_byte(gebco, scale = TRUE)))

ximage::ximage(read_to_nr(make_byte(nsidc), 512))
ximage::ximage(read_to_nr(make_byte(nsidc, scale = TRUE), 512))

The scale arg in vrt:// can take 1) true/false 2) an output range or 3) an input+output range too. I haven't toyed much with the GDAL colour tools but we can obviously pipeline quite a lot or even extend "vrt://" to provide more control.

Ah, we can easily expand palette (but they won't all be Byte(-limited) so that needs some thought).

make_byte <- function(x, scale = FALSE, expand = "") {
  x <- sprintf("vrt://%s?ot=Byte", x)
  if (scale) {
    x <- sprintf("%s&scale=true", x)
  }
  if (nzchar(expand)) {
    if (!expand %in% c("gray", "rgb", "rgba")) {
      stop("value of 'expand' unsupported. Only gray, rgb or rgba are supported.")
    }
    x <- sprintf("%s&expand=%s", x, expand)
  }
  
  x
}


ximage::ximage(read_to_nr(make_byte(nsidc, expand = "rgb"), 512), asp = 1)

Does that answer?

Yes I'm with you. Thanks for the well-considered reply and demo code. I could have thought that out better before prodding you to walk me through it. I put in a quantile option in plot_raster specifically for 16-bit after all. Appreciate the explanation connecting several pieces.

nativeRaster is more interesting than string encoding every pixel, but it's (now) obviously not so simple extending the data type. Agree it's worth exploring more and glad you added direct read so that gdalraster has that basis now. Will spend more time with the code and pick up on this again.