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