I don't entirely trust the timings here, but it's faster (maybe 2x faster to plot).
dsn <- "WMTS:https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"
library(gdalraster)
#> GDAL 3.9.0dev-cb4d30f56d, released 2024/04/15 (debug build), GEOS 3.12.1, PROJ 9.3.1
ds_int <- new(GDALRaster, sprintf("vrt://%s?ot=Int32", dsn))
ds_byt <- new(GDALRaster, dsn)
nshape <- c(ds_int$getRasterXSize(), ds_int$getRasterYSize())
shape <- as.integer(round(1024 * c(1, nshape[2]/nshape[1])))
pix_int <- read_ds(ds_int, bands = 1:3, out_xsize = shape[1], out_ysize = shape[2])
## we get Bytes via this PR: https://github.com/USDAForestService/gdalraster/pull/314
pix_byt <- read_ds(ds_byt, bands = 1:3, out_xsize = shape[1], out_ysize = shape[2])
## smaller
object.size(pix_int)
#> 12548360 bytes
object.size(pix_byt)
#> 3138824 bytes
## the encoding to nativeRaster is not complicated, implemented in R here (see {png}, {farver}, and coolbutuseless/{nara})
## https://github.com/wch/r-source/blob/6bbb6d5ac1c4bba97c941436a1beb3e6c1118a3a/src/include/R_ext/GraphicsDevice.h#L869-L870
to_native <- function(red, green, blue) {
as.integer(red + green * 256L + blue * (256L*256L) - 256^3 )
}
nm <- structure(matrix(do.call(to_native, setNames(split(as.integer(pix_byt), rep(1:3, each = prod(attr(pix_byt, "gis")$dim[1:2]))),
c("red", "green", "blue"))),
shape[2]), class = "nativeRaster", channels = 3)
## nativeRaster is faster to plot
system.time(ximage::ximage(nm, extent = attr(pix_byt, "gis")$bbox[c(1, 3, 2, 4)], asp = 1))
#> user system elapsed
#> 0.031 0.004 0.033
system.time(plot_raster(pix_int))
#> user system elapsed
#> 0.584 0.027 0.612
Created on 2024-04-17 with reprex v2.0.2