Using new nativeRaster support in gdalraster from github, pull a reprojected image from a tile server direct to R's internal image encoding"

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}