EPSG:28992 tiles?

[gremms1]

Hello,

I am trying to get the following code to work:

library(wmts)
p = st_point(c(110784,494752))
p = st_sfc(p)
p = st_set_crs(p, "EPSG:28992")
bb = st_bbox(st_buffer(p, 300))
url <- "https://services.arcgisonline.nl/arcgis/rest/services/Basiskaarten/Open_Topo/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"
t <- wmts(x=url, loc = bb)

However, this is giving me a lot of errors. Is it even possible to use a custom CRS (EPSG:28992 in this case)? Thanks in advance!

[mdsumner]

what are the errors?

I'd recommend using hypertidy/{gdalio} instead of {wmts}, I'll try an example - thanks for the motivation I should archive this project 🙏

[mdsumner]

this works

##remotes::install_github("hypertidy/gdalio")
url <- "https://services.arcgisonline.nl/arcgis/rest/services/Basiskaarten/Open_Topo/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"

library(gdalio)
ex <- c(110784,494752)[c(1, 1, 2, 2)] + c(-1, 1, -1, 1) * 300
gdalio_set_default_grid(list(extent = ex, 
                             dimension = c(1024, 1024), 
                             projection = "EPSG:28992"))
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))

r <- gdalio_raster(url, bands = 1:3)
library(raster)
plotRGB(r)

HTH, this project {wmts} was an exploration along the road to {gdalio} -it's not perfect but gdalio is way more general and powerful :)

one of the hardest things is getting the right dimension for your given graphics context, you can use dev.size("px") but because you have a data aspect-ratio you need to modify that to suit and I have code somewhere to help with that

the other thing to consider is the resample argument, by default it's just nearest neighbour

https://github.com/hypertidy/gdalio

[gremms1]

Hi, That's amazing! Thank you so much for taking your time to provide a working example! Is it also possible to define the zoom-level somehow? The example you provide is a bit gritty, but there are more zoomed in tiles available as well (see picture)

[mdsumner]

it will always be a bit gritty if you reproject it, you can get the source image exactly by specifying a grid aligned exactly to the source itself, 28992 is a oblique stereographic not Mercator

I would try different resample algs,but the key is requesting the tile at exactly the size you are going to render it (and this is a massive under-supported thing in R itself)

[mdsumner]

oh wait, sorry - the data is in 28992, I'd taken from your OP that it was not

I'll put together an example that gets grid exactly, you still have the problem of rendering that faithfully but I can't help that :)

[gremms1]

Ah yes sorry for not making that exactly clear. Indeed, the tiles are in 28992. Thanks again! I have been trying to find a simple way of making static maps in R with these tiles for ages. This is great.

[mdsumner]

this works, I don't yet have a simpler way of finding the best zoom level to specify the grid target, there's a lot of back and forth here that is hard to explain quickly - but consider that even though this resulting raster is an exact extraction from the given tile level, when you plot it the R device is possibly reshaping that into a different device size (...)

(and this code should be more automated but it's something I've been wondering how to package up)

HTH

ri <- vapour::vapour_raster_info(url)
## work down the stack until we get something near sensible
library(raster)
ov <- matrix(c(ri$dimXY, ri$overviews), ncol = 2, byrow = TRUE)
maxdim <- 2048  ## 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)
plotRGB(r)

[gremms1]

Thank you so much, however, the first line is giving me an error:

Error in raster_info_gdal_cpp(dsn = datasourcename, min_max = min_max) : 
  cannot open dataset
In addition: Warning message:
In raster_info_gdal_cpp(dsn = datasourcename, min_max = min_max) :
  GDAL Error 1: Driver COUCHDB is considered for removal in GDAL 3.5. You are invited to convert any dataset in that format to another more common one .If you need this driver in future GDAL versions, create a ticket at https://github.com/OSGeo/gdal (look first for an existing one first) to explain how critical it is for you (but the GDAL project may still remove it), and to enable it now, set the GDAL_ENABLE_DEPRECATED_DRIVER_COUCHDB configuration option / environment variable to YES

[mdsumner]

hmmm - did you change computers possibly?

try prepending WMTS: i.e.

 url <- "WMTS:https://services.arcgisonline.nl/arcgis/rest/services/Basiskaarten/Open_Topo/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"
ri <-  vapour_raster_info(url)

a reprex:

  url <- "https://services.arcgisonline.nl/arcgis/rest/services/Basiskaarten/Open_Topo/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"

ri <- vapour::vapour_raster_info(url)
## work down the stack until we get something near sensible
library(raster)
#> Loading required package: sp
ov <- matrix(c(ri$dimXY, ri$overviews), ncol = 2, byrow = TRUE)
maxdim <- 2048  ## 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
ex0 <- c(110784,494752)[c(1, 1, 2, 2)] + c(-1, 1, -1, 1) * 300

for (i in 1:nrow(ov)) {
 template <- crop(raster(extent(ri$extent), nrows = ov[i, 2], ncols = ov[i, 1]), ex0)
 dm <- dim(template)[2:1]
 ex0 <- extent(template)
 ex0 <- c(xmin(ex0), xmax(ex0), ymin(ex0), ymax(ex0))
 if (max(dm) <= maxdim) break;
}

library(gdalio)
## 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)
plotRGB(r)

^{Created on 2022-02-22 by the reprex package (v2.0.1)}

I'm thinking about how to expose that raster-align logic, I have it all in packages how I want it but it's still easier to lean on raster for that ...

[gremms1]

Adding "WMTS:" to the url did the trick, thanks! The map it produces is still more gritty than the interactive map i posted earlier though. Can I play with the maxdim variable to find a better looking zoom-level?

[mdsumner]

yes play with that, and try gdalio_raster(url, resample = "cubic") also

it totally depends how you render it too :)

[gremms1]

Yes that makes sense. Thanks again! Cheers