Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Created January 22, 2026 00:15
Show Gist options
  • Select an option

  • Save mdsumner/3f5ed392ad126ad91dd403473c669bdd to your computer and use it in GitHub Desktop.

Select an option

Save mdsumner/3f5ed392ad126ad91dd403473c669bdd to your computer and use it in GitHub Desktop.
  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

@mdsumner
Copy link
Author

mdsumner commented Jan 23, 2026

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)
image

@ctoney
Copy link

ctoney commented Jan 24, 2026

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment