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}

Using new nativeRaster support in gdalraster:

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"

ex <- reproj::reproj_extent(c(143, 149, -44, -39), "EPSG:3031", source = "EPSG:4326")

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)))
}
(dm <- fit_dims(1024, wh = diff(ex)[c(1L, 3L)]))
#> [1] 1004 1024

px <- c(diff(ex[1:2])/dm[1L], -diff(ex[3:4])/dm[2L])
gt <-   c(xmin = ex[1L],
          xres = px[1L],
          yskew = 0,
          ymax = ex[4L],
          xskew = 0,
          yres = px[2L])

ds <- create("MEM", 
             tf <- tempfile(fileext = ".mem", tmpdir = "/vsimem"), 
             xsize = dm[1L], ysize = dm[2L], nbands = 3L, dataType = "UInt8", return_obj = TRUE)
ds$setProjection("EPSG:3031")
#> [1] TRUE
ds$setGeoTransform(gt)
#> [1] TRUE

## do the job
warp(dsn, ds, "")
#> 0...10...20...30...40...50...60...70...80...90...100 - done.

## helper function to plot a nativeRaster
## with optional extent to set up plot frame (if add = FALSE) and render to that subwindow
## extent = xmin,xmax,ymin,max
image_nr <- function(x, extent = NULL, add = FALSE) {
  if (is.null(extent)) {
    extent <- c(0, nrow(x), 0, ncol(x))
  }
  if (!add) plot(NA, xlab = "", ylab = "", xlim = extent[1:2], ylim = extent[3:4], asp = 1)
  rasterImage(x, extent[1], extent[3], extent[2], extent[4])
}
image_nr(read_to_nativeRaster(ds), extent = ds$bbox()[c(1, 3, 2, 4)])

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

Very nice. Always appreciate your examples.

Small detail... when using MEM the dataset doesn't actually have a real file name. It's held in a memory allocation managed by the driver and not in /vsimem. I usually pass either an empty string or something like "in-memory raster" for the filename argument. It's possibly better not to pass a real filename to avoid confusion because the string that you pass to create() for example, will appear as the "DSN" value from the object's show() method. That's because I set it here https://github.com/firelab/gdalraster/blob/ebada5f2a1eb17eda8170fb1c9ddbf458eba02a0/src/gdal_exp.cpp#L847 and it then gets set as the GDAL dataset-level Description, which is what show() reports. Normally that's the dataset's filename. MEM is the exception, but minor detail.

Since nativeRaster can be passed to rasterImage(), it seems like gdalraster::plot_raster() could have a code path accepting a nativeRaster input, then we wouldn't need the helper function. I was going to look into that once I have a chance to experiment.

Ah thanks I had meant to change it to GTiff as I knew that was a bit weird, I'll follow that up. I had considered including the $gis metadata on the nativeRaster object for use by plot_raster, just the idea of informal subclassing seemed a bit off (it's probably a harmless shortcut for plot_raster though).

the idea of informal subclassing seemed a bit off

Yeah I agree. It should be S3 classed if we actually implemented something. That could be worthwhile though and start over to use nativeRaster whenever possible. Hugh Graham made several improvements to plot_raster for vrtility. He mentioned at one point that he was trying to keep it as compatible with the original as possible in case it's of interest. I haven't followed up with him but maybe the nativeRaster addition is a good time to consider integrating.

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?

Oh totally that's a very attractive idea. I think it's a very hard thing, to be general and lightweight and abstract as possible but also able to support custom map making - it's too much for one function, but I'm not sure anyone in R has got this really figured out (or anywhere). Think of magick - that's such cool ux but it's entirely incompatible atm. I'm pretty tied to my own hypertidy/ximage, and I like the idea of raster graphics that aren't tied to GDAL and building on that. I'd love to brainstorm it out and I don't think it's out of aspiration range to cut across a lot of other packages too. The core in R is really rasterImage and nativeRaster, so maybe there's a need for plot() method on array/matrix as a basis. I'm in love with ximage(x, extent-optional) so maybe that (or equivalen) is a contribution to R itself.

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.