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}

Another possibility is a corespatial package where we keep a foundation for things like ximage (and grid logic, graphics interaction, recording context (CRS) with a plot, etc). The only rule would be "no dependencies" 😆 But there's always the problem of where S3 definitions belong but if it's just for internal dev use maybe that's ok.

Another possibility is a corespatial package where we keep a foundation for things like ximage (and grid logic, graphics interaction, recording context (CRS) with a plot, etc)

That seems like a good idea. I also like the suggestion that we ponder this some and brainstorm it out. I'm not going to add a heavy dependency, especially on another core R spatial package, just for plotting. At the same time, having basic display capability within the package (or at least with very light dependence) is basically a requirement for a few reasons. But also at the same time, I can't and don't really want to invest heavily in development and maintenance of it. I would likely collaborate with others on that but I can't seeing having the time or sustained interest to do it properly on my own.

But on the other question, I was asking about readToNativeRaster itself:

It doesn't have to limit to 8-bit data does it? Couldn't it read 16-bit multi-band imagery and normalize to 8-bit range before encoding?

(or any data type for that matter, or 16-bit grayscale for example, any 1/3/4 band that wants to be treated as grayscale/rgb/rgba)

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.