wmts_gdal.md

dsn <- "WMTS:https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"

ex <- c(-1, 1, -1, 1) * 5000 ## 5km
pt <- cbind(147, -42)         ## lon,lat anywhere you want
res <- c(10, 10) ## 10m pixels, obvsly update this inline with the buffer above to get the nx,ny dimension you actually want
## use a local projection, because this works anywhere you only need modify ex and res (or dim: ncols/nrows)
## we can use longlat or mercator or whatever, just set ex and res/dim
crs <- sprintf("+proj=laea +lon_0=%f +lat_0=%f +datum=WGS84", pt[1], pt[2])

im <- vapour::gdal_raster_nara(dsn, target_ext = ex, target_crs = crs, target_res = res)
ximage::ximage(im, asp = 1)

## equivalent to (need a workaround with sds because of https://github.com/rspatial/terra/issues/1735)
library(terra)
src <- sds(dsn, 1)[[1]]  ## don't work with this, it's volatile and potentially huge
## src <- rast(dsn)  
template <- rast(ext(ex), crs = crs, res = res, nlyrs = nlyr(src))
im <- project(src, template, by_util = TRUE, method = "cubic")
plotRGB(im)