Pdok WMTS raster tiles downloaden in R

_thomas · 12 februari 2022 om 20:57

Hallo, ik wil graag kaarten maken in R, en zou daar de OpenTopo kaart (link) voor willen gebruiken. Ik heb veel op het internet gezocht maar het is me niet gelukt om iets werkends te reproduceren. Wat ik wel heb kunnen vinden is dat met de package maptiles je in principe raster bestanden van een WMTS service moet kunnen downloaden.
Ik heb het volgende werkende voorbeeld hiervan gevonden:

library(maptiles)
library(sf)
library(raster)

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))

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)

rgb = stack(rgb)
plotRGB(rgb)

Wat ik zou willen is dus exact wat hierboven gebeurt maar dan met de OpenTopo kaarten. Ik hoop dat iemand mij hiermee kan helpen.

Groet, Thomas

wouter.visscher · 14 februari 2022 om 07:02

Hoi @_thomas, volgens mij ben je opzoek naar de volgende URL (laatste stukje aanpassen met de proj+zxy template) →

https://geodata.nationaalgeoregister.nl/tiles/service/wmts/opentopoachtergrondkaart/EPSG:28992/5/12/16.png

https://geodata.nationaalgeoregister.nl/tiles/service/wmts/opentopoachtergrondkaart/EPSG:28992/5/12/16.png

Dit is de RESTful variant van een WMTS request, degene in het NGR (de link in je post) is de KVP variant daarvan. Hopelijk dat het wel met de bovenstaande lukt. Laat je het weten of het wel/niet is gelukt. Dat helpt andere die een WMTS in R willen laden ook weer verder met zoeken.

_thomas · 14 februari 2022 om 07:57

Hoi Wouter,

Bedankt voor de snelle respons! Ik denk dat ik met de link die jij aangaf al een stukje meer in de buurt ben, maar met de volgende code heb ik het nog niet werkend gekregen:

    library(sf)
    library(maptiles)
    p = st_point(c(110784.2,494605.3))
    p = st_sfc(p)
    p = st_set_crs(p, 28992)
    bb = st_bbox(st_buffer(p, 300))

    url <- "https://geodata.nationaalgeoregister.nl/tiles/service/wmts/opentopo/EPSG:28992/{z}/{y}/{x}.png"

    mffp <- list(src="OpenTopo", q=url, sub=NA, cit='')

    t <- get_tiles(bb, provider=mffp, crop= T, forceDownload = T, zoom = 5)

Ik krijg hier de foutmelding “Error in png::readPNG(img) : file is not in PNG format”. Wat nogal curieus is want voor zover ik weet zijn de tiles wel degelijk in png-formaat. Als ik .jpg gebruik krijg ik een HTTP 400 error.

wouter.visscher · 14 februari 2022 om 08:03

Je zit mogelijk “verkeerd” in het Tilegrid
De foutmelding “Error in png::readPNG(img) : file is not in PNG format” is mogelijk dat je XML terugkrijgt

Bijvoorbeeld door een z,x,y waard te geven die te groot is voor het grid. zie response bij https://geodata.nationaalgeoregister.nl/tiles/service/wmts/opentopoachtergrondkaart/EPSG:28992/5/12/10000000006.png (<- extreem voorbeeld)

Als je een ‘fout’ request weet te loggen en deze in een browser gooit kan je verifieren of dat gaande is.

Gezien je het een en ander al werkend had voor EPSG:3857 zou je die mogelijk als start punt kunnen gebruiiken

https://geodata.nationaalgeoregister.nl/tiles/service/wmts/opentopo/EPSG:3857/{z}/{y}/{x}.png

_thomas · 14 februari 2022 om 12:13

Dat dacht ik ook, maar dat lijkt toch niet zo te zijn. Ook als ik https://geodata.nationaalgeoregister.nl/tiles/service/wmts/opentopoachtergrondkaart/EPSG:28992/5/12/16.png als url gebruik en crop=FALSE (dan wordt de bounding box genegeerd) krijg ik dezelfde error. Wat betreft de andere projectie: als ik https://geodata.nationaalgeoregister.nl/tiles/service/wmts/opentopo/EPSG:3857/5/12/16.png opzoek in de browser krijg ik de error unknown tilematrixset: EPSG:3857 dus dat lijkt mij ook geen optie? Misschien moet ik even kijken of ik contact kan krijgen met de beheerders van de maptiles package. Of weet hier iemand een andere package die ik nog zou kunnen proberen?

wouter.visscher · 14 februari 2022 om 12:29

Hmm, sorry dacht dat we opentopo ook als EPSG:3857 hadden (maar dus alleen in EPSG:28992)

Dezelfde python error, of dezelfde error (XML) van de server?

_thomas · 14 februari 2022 om 12:58

Dezelfde error in R.

_thomas · 17 februari 2022 om 10:06

Update: ik heb het min of meer werkend gekregen met de volgende code:
library(tmap)
library(tmaptools)

data("NLD_prov")
url <- "https://geodata.nationaalgeoregister.nl/tiles/service/wmts/opentopo/EPSG:28992/{z}/{x}/{y}.jpeg"
map <- tmaptools::read_osm(NLD_prov, type = url, mergeTiles = T, extension="jpeg")

tmap_mode("plot")
     tm_shape(map) +
      tm_rgb() +
     tm_shape(NLD_prov) +
      tm_polygons() +
     tm_graticules(lines=F)

Maar, zoals te zien is gaat het wel ergens mis met de projectie. Als iemand een idee heeft hoe ik dit kan verhelpen hoor ik het graag.

_thomas · 17 februari 2022 om 10:14

Het is mij wel gelukt om een interactieve kaart te produceren met de juiste projectie (de epsg28992 code ben ik op een ander forum tegengekomen en heb de resolutions hier vandaan).

library(tmap)
library(tmaptools)
library(leaflet)

epsg28992 <- leafletCRS(
  crsClass = "L.Proj.CRS"
  ,code = "EPSG:28992"
  ,proj4def = "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs"
  ,origin=c(-285401.92,903401.92)
  ,resolutions = c(3251.206502413005,1625.6032512065026,812.8016256032513,406.40081280162565,203.20040640081282,101.60020320040641, 50.800101600203206,25.400050800101603,12.700025400050801,6.350012700025401,3.1750063500127004,1.5875031750063502,0.7937515875031751,0.39687579375158755,0.19843789687579377,0.09921894843789689,0.04960947421894844)
)


url <- "https://geodata.nationaalgeoregister.nl/tiles/service/wmts/opentopo/EPSG:28992/{z}/{x}/{y}.jpeg"
data("World")
tmap_mode("view")
tm_shape(World) +
  tm_polygons() +
  tm_basemap(NULL) +
  tm_view(projection = epsg28992, set.view = c(5.092098,52.093992,4)) +
  tm_tiles(server = url) +
  tm_mouse_coordinates()

_thomas · 22 februari 2022 om 19:29

Met wat hulp van derden is het inmiddels ook gelukt om een statische kaart te maken.
Hieronder de code voor de geïnteresseerden:

#remotes::install_github("hypertidy/gdalio")
library(gdalio)
library(stars)
library(vapour)
library(raster)
library(tmap)
library(sf)


ex <- st_point(c(4.7,52)) %>%
  st_sfc() %>%
  st_set_crs(4326) %>%
  st_transform(crs=28992) %>%
  st_buffer(10000) %>%
  st_as_sf()


url <- "WMTS:https://geodata.nationaalgeoregister.nl/tiles/service/WMTS/1.0.0/WMTSCapabilities.xml,layer=opentopo"
ri <-  vapour_raster_info(url)
## work down the stack until we get something near sensible
ov <- matrix(c(ri$dimXY, ri$overviews), ncol = 2, byrow = TRUE)
maxdim <- 750  ## note that we want the first overview that *crops* to this maximum
## and we also need the aligned extent for this, to get exact pix
for (i in 1:nrow(ov)) {
  template <- crop(raster(extent(ri$extent), nrows = ov[i, 2], ncols = ov[i, 1]), ex)
  dm <- dim(template)[2:1]
  ex0 <- extent(template)
  ex0 <- c(xmin(ex0), xmax(ex0), ymin(ex0), ymax(ex0))
  if (max(dm) <= maxdim) break;
}

## now we know the right dm for the right extent to get it exact
gdalio_set_default_grid(list(extent = ex0, 
                             dimension = dm, 
                             projection = "EPSG:28992"))
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
r <- gdalio_raster(url, bands = 1:3)

map <- st_as_stars(r)

tmap_mode("plot")

map_static <- 
  tm_shape(map, raster.downsample = T) +
  tm_rgb() +
  tm_graticules(lines=F, n.x = 3, n.y = 3)

map_static

nvdoormaal · 18 december 2023 om 15:35

Hoi @_thomas ,

Gebruik je de bovenstaande code nog? Ik ben ook aan het worstelen om WMS / WMTS als statische achtergrond te gebruiken voor kaartjes in R.
Wanneer ik je code gebruik, krijg ik helaas de volgende foutmelding:
“gdalinfo - unable to open ‘WMTS:https://geodata.nationaalgeoregister.nl/tiles/service/WMTS/1.0.0/WMTSCapabilities.xml,layer=opentopo’.
ERROR 1: HTTP error code : 404
ERROR 1: Missing Capabilities.Contents element
Error in vapour_raster_info(url) : GDAL was unable to open ^^”

Daarom was ik benieuwd of de code bij jou nog wel werkt of ondertussen hebt aangepast / verder verfijnd.

Ik hoor het graag!

sbjager · 19 december 2023 om 07:27

Die foutmelding komt doordat Opentopo niet meer in productie is.

Je zult dus een andere achtergrondkaart moeten kiezen.