Pdok WMTS raster tiles downloaden in R

Toepassingen
#1

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

#2

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.

#3

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.

#4

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)

image
image1067×278 24.7 KB

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
1 like
#5

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?

#6

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?

#7

Dezelfde error in R.

#8

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)

image
image704×552 54.7 KB

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.

#9

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()
2 likes
#10

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
4 likes
#11

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!

#12

Die foutmelding komt doordat Opentopo niet meer in productie is.

Je zult dus een andere achtergrondkaart moeten kiezen.

3 likes