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}

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.