You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
I'm not sure if this should be an issue or a feature request, but here it goes anyway.
For mapping in the Netherlands, EPSG:28992 is used by default by Dutch government agencies. These agencies also provide WMTS basemaps, both in EPSG:28992 and EPSG:3857. However, it appears that annotation_map_tile is hardcoded to only handle EPSG:3857 tiles.
This raises the following issue: when mixing and matching maps and tiles in different CRSs, it seems to matter whether one transforms the map data to the tile CRS, or whether one first adds an annotation_map_tile before an explicit data argument is mentioned in the ggplot2 list of layers. Reprex follows below. What explains this behavior? Am I doing something unidiomatic here?
Ideally, I'd like to pass data to ggplot() directly without having to transform it to the tile CRS. So as a feature request: can annotation_map_tile (or the underlying rosm machinery) be made to work with arbitrary tile CRSs instead of EPSG:3857 only?
library(assertthat)
library(ggspatial)
library(httr)
library(rosm)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(tidyverse)
# Download official map of the Netherlands# NOTE: uses EPSG:28992url<- parse_url("https://geodata.nationaalgeoregister.nl/cbsgebiedsindelingen/wfs")
url$query<-list(
service="WFS",
version="2.0.0",
request="GetFeature",
typename="cbsgebiedsindelingen:cbs_gemeente_2021_gegeneraliseerd",
outputFormat="application/json"
)
NLD_muni=url %>%
build_url() %>%
st_read()
#> Reading layer `OGRGeoJSON' from data source #> `https://geodata.nationaalgeoregister.nl/cbsgebiedsindelingen/wfs?service=WFS&version=2.0.0&request=GetFeature&typename=cbsgebiedsindelingen%3Acbs_gemeente_2021_gegeneraliseerd&outputFormat=application%2Fjson' #> using driver `GeoJSON'#> Simple feature collection with 352 features and 5 fields#> Geometry type: MULTIPOLYGON#> Dimension: XY#> Bounding box: xmin: 13565.4 ymin: 306846.2 xmax: 278026.1 ymax: 619352.4#> Projected CRS: Amersfoort / RD NewNLD_muni %>% ggplot() + geom_sf()
# Plot 2: pre-transform the data to EPSG:3857 and then add basemap# Southern border correctly runs through the middle of the river & canal.# Northern border correctly runs through the middle of the railroad.muni_321 %>% st_transform(3857) %>% ggplot() +
annotation_map_tile(type="brta_standaard", zoomin=+1) +
geom_sf(alpha=.5)
#> Zoom: 14
# Plot 3: call bare ggplot(), add basemap and pass data to geom_sf()# Southern border correctly runs through the middle of the river & canal.# Northern border correctly runs through the middle of the railroad.
ggplot() +
annotation_map_tile(type="brta_standaard", zoomin=+1) +
geom_sf(data=muni_321, alpha=.5)
#> Zoom: 14
Thanks for raising this issue...normally the CRS of the plot is set according to the first layer, which would be EPSG:3857 by default for annotation_map_tile(). I think that is functioning as intended here but it is very odd that different transformation are getting applied in the first and third plots. I'll have time to look into this on Monday!
You're right that the backend for annotation map tile is hard-coded to get tiles in EPSG 3857...I'd like to improve the backend soon but don't quite have the infrastructure in place to make it happen.
So if I'm understanding this correctly, in the first plot, the data is in 28992, but for the tile requests that data is internally transformed into 3857, the tiles are downloaded and stitched together, and then overlaid on the original data in 28992, which gives small (we're talking about 10 meters here) errors. Ideally, the fact that the tile offsets have been computed in 3857 should be inferred at drawing time. I have no idea how that works internally in {ggspatial}.
In any case, for these RESTful WMTS servers, there is typically an XML document such as this one (notice the same base URL) that has all the required information about the image scales per available zoom level per TileMatrixSet: https://service.pdok.nl/brt/achtergrondkaart/wmts/v2_0/WMTSCapabilities.xml
It would be nice if rosm::register_tile_source() would take base_url, feature and capabilities arguments. Then annotation_map_tile can first query the XML to see if tiles "native" to the requesting layer can be downloaded before computing the various tile offsets. With a fallback to 3857 if none such native tile set is available.
I'm not sure if this should be an issue or a feature request, but here it goes anyway.
For mapping in the Netherlands, EPSG:28992 is used by default by Dutch government agencies. These agencies also provide WMTS basemaps, both in EPSG:28992 and EPSG:3857. However, it appears that
annotation_map_tile
is hardcoded to only handle EPSG:3857 tiles.This raises the following issue: when mixing and matching maps and tiles in different CRSs, it seems to matter whether one transforms the map data to the tile CRS, or whether one first adds an
annotation_map_tile
before an explicitdata
argument is mentioned in theggplot2
list of layers. Reprex follows below. What explains this behavior? Am I doing something unidiomatic here?Ideally, I'd like to pass data to
ggplot()
directly without having to transform it to the tile CRS. So as a feature request: canannotation_map_tile
(or the underlyingrosm
machinery) be made to work with arbitrary tile CRSs instead of EPSG:3857 only?Created on 2021-08-25 by the reprex package (v2.0.1)
The text was updated successfully, but these errors were encountered: