R – Download Tiles Using Custom WMTS Provider

Question

I'm using a WMTS layer in QGIS. Everything works fine. Now I would like to get access to the images in R both using maptiles and mapview. I do not know the WMTS standard and consequently I hardly understand the layer capability metadata but I eventually built a valid url from that xml file. The code do download something but it is always in a somewhat random location. The XML capability file is here.

library(maptiles)
library(sf)

p = st_point(c(63500, 624900))
p = st_sfc(p)
p = st_set_crs(p, 32198)
p = st_transform(p, 4326)
bb = st_bbox(st_buffer(p, 300))

# Good
rgb = get_tiles(bb, provider = "Esri.WorldImagery")
library(terra)
plot(rgb)

mffp <- list(
  src = 'Aerial photography',
  q = 'https://servicesmatriciels.mern.gouv.qc.ca:443/erdas-iws/ogc/wmts/Inventaire_Ecoforestier/Inventaire_Ecoforestier/default/OGC:1.0:GlobalCRS84Pixel/{z}/{x}/{y}.jpg',
  sub = NA,
  cit = ''
)

# In the ocean (maybe)
rgb = get_tiles(bb, provider = mffp)
plot(rgb)

Notice that I'm keen to bypass maptile and do the math myself to retrieve the tiles if necessary but at this stage I don't even know what to compute. Also I did not find yet how to add a provider in mapview but if I can work with maptile it's already good enough.

Answer 1

The x and y are the other way round for the Google Maps Compatible layer:

url="https://servicesmatriciels.mern.gouv.qc.ca:443/erdas-iws/ogc/wmts/Inventaire_Ecoforestier/Inventaire_Ecoforestier/default/GoogleMapsCompatibleExt2:epsg:3857/{z}/{y}/{x}.jpg"
mffp = list(src="Aerial Photography", q=url, sub=NA, cit='')
rgb = get_tiles(bb, provider = mffp)
plot(rgb)

complains because there's some NAs in there:

Error in grDevices::rgb(RGB[, 1], RGB[, 2], RGB[, 3], alpha = alpha, maxColorValue = scale) : 
  color intensity NA, not in 0:255

whateven is a SpatRaster object anyway? Let's make it a raster package stack:

library(raster)
rgb = stack(rgb)
plotRGB(rgb)

(ah, SpatRaster is from the terra package. I guess there's a way to deal with NAs but I dunno yet. The raster package seems happy enough for now with NAs in there, and clearly the crux of the problem is solved.)

Figured this out by putting the URL into QGIS as an XYZ Tiles layer - this is what you get if the {x} and {y} are the wrong way round: