1

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.

Improve this question
2
  • That WMTS seems to supply 5 different layers - some are global and some are just a bit of Canada. There's various different coordinate systems. I suspect maptile is expecting everything to be in Google Maps' Web Mercator... Have you tried the image URL for the GoogleMapsCompatibleExt2:epsg:3857 tile set?
    – Spacedman
    Commented Oct 12, 2021 at 11:38
  • Yes, I reached the same conclusion. I tried GoogleMapsCompatibleExt2:epsg:3857 but maptile is expecting lon/lat WGS84 (epsg:4326) so it does not works better in my tests. As I said I could do the math if I were able to figure out the proper inputs (bbox modulo 256px tiles + zoom resolution should be enough I guess) .
    – JRR
    Commented Oct 12, 2021 at 12:17

1 Answer 1

Reset to default
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)

enter image description here

(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:

enter image description here

Improve this answer
1
  • So simple. Thank you for your help
    – JRR
    Commented Oct 12, 2021 at 14:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.