paleolimbot / ggspatial

Enhancing spatial visualization in ggplot2
https://paleolimbot.github.io/ggspatial
368 stars 34 forks source link

annotation_map_tile glitches when working with CRS other than EPSG:3857 #89

Open rhalbersma opened 2 years ago

rhalbersma commented 2 years ago

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:28992
url <- 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 New

NLD_muni %>% ggplot() + geom_sf()

# Register official WMTS server for basemap of the Netherlands
# NOTE: the server provides both EPSG:3857 and EPSG:28992 tiles, rosm only works with the former.
brta_standaard = source_from_url_format(
  url_format = "https://service.pdok.nl/brt/achtergrondkaart/wmts/v2_0/standaard/EPSG:3857/${z}/${x}/${y}.png",
  attribution = "BRT achtergrondkaart Kadaster CC-BY 4.0"
)
register_tile_source(brta_standaard = brta_standaard)
assertthat::assert_that("brta_standaard" %in% osm.types())
#> [1] TRUE

# Home sweet home
muni_321 = NLD_muni %>% filter(statcode == "GM0321")

# Plot 1: zoom in on a particular municipality and add basemap
# Southern border incorrectly runs along the river and canal banks.
# Northern border incorrectly runs just east of the railroad.
muni_321 %>% ggplot() + 
  annotation_map_tile(type = "brta_standaard", zoomin = +1) + 
  geom_sf(alpha = .5)
#> Loading required namespace: raster
#> Zoom: 14
#> Fetching 81 missing tiles
#>   |                                                                              |                                                                      |   0%  |                                                                              |=                                                                     |   1%  |                                                                              |==                                                                    |   2%  |                                                                              |===                                                                   |   4%  |                                                                              |===                                                                   |   5%  |                                                                              |====                                                                  |   6%  |                                                                              |=====                                                                 |   7%  |                                                                              |======                                                                |   9%  |                                                                              |=======                                                               |  10%  |                                                                              |========                                                              |  11%  |                                                                              |=========                                                             |  12%  |                                                                              |==========                                                            |  14%  |                                                                              |==========                                                            |  15%  |                                                                              |===========                                                           |  16%  |                                                                              |============                                                          |  17%  |                                                                              |=============                                                         |  19%  |                                                                              |==============                                                        |  20%  |                                                                              |===============                                                       |  21%  |                                                                              |================                                                      |  22%  |                                                                              |================                                                      |  23%  |                                                                              |=================                                                     |  25%  |                                                                              |==================                                                    |  26%  |                                                                              |===================                                                   |  27%  |                                                                              |====================                                                  |  28%  |                                                                              |=====================                                                 |  30%  |                                                                              |======================                                                |  31%  |                                                                              |======================                                                |  32%  |                                                                              |=======================                                               |  33%  |                                                                              |========================                                              |  35%  |                                                                              |=========================                                             |  36%  |                                                                              |==========================                                            |  37%  |                                                                              |===========================                                           |  38%  |                                                                              |============================                                          |  40%  |                                                                              |=============================                                         |  41%  |                                                                              |=============================                                         |  42%  |                                                                              |==============================                                        |  43%  |                                                                              |===============================                                       |  44%  |                                                                              |================================                                      |  46%  |                                                                              |=================================                                     |  47%  |                                                                              |==================================                                    |  48%  |                                                                              |===================================                                   |  49%  |                                                                              |===================================                                   |  51%  |                                                                              |====================================                                  |  52%  |                                                                              |=====================================                                 |  53%  |                                                                              |======================================                                |  54%  |                                                                              |=======================================                               |  56%  |                                                                              |========================================                              |  57%  |                                                                              |=========================================                             |  58%  |                                                                              |=========================================                             |  59%  |                                                                              |==========================================                            |  60%  |                                                                              |===========================================                           |  62%  |                                                                              |============================================                          |  63%  |                                                                              |=============================================                         |  64%  |                                                                              |==============================================                        |  65%  |                                                                              |===============================================                       |  67%  |                                                                              |================================================                      |  68%  |                                                                              |================================================                      |  69%  |                                                                              |=================================================                     |  70%  |                                                                              |==================================================                    |  72%  |                                                                              |===================================================                   |  73%  |                                                                              |====================================================                  |  74%  |                                                                              |=====================================================                 |  75%  |                                                                              |======================================================                |  77%  |                                                                              |======================================================                |  78%  |                                                                              |=======================================================               |  79%  |                                                                              |========================================================              |  80%  |                                                                              |=========================================================             |  81%  |                                                                              |==========================================================            |  83%  |                                                                              |===========================================================           |  84%  |                                                                              |============================================================          |  85%  |                                                                              |============================================================          |  86%  |                                                                              |=============================================================         |  88%  |                                                                              |==============================================================        |  89%  |                                                                              |===============================================================       |  90%  |                                                                              |================================================================      |  91%  |                                                                              |=================================================================     |  93%  |                                                                              |==================================================================    |  94%  |                                                                              |===================================================================   |  95%  |                                                                              |===================================================================   |  96%  |                                                                              |====================================================================  |  98%  |                                                                              |===================================================================== |  99%  |                                                                              |======================================================================| 100%
#> ...complete!

# 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

Created on 2021-08-25 by the reprex package (v2.0.1)

paleolimbot commented 2 years ago

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.

rhalbersma commented 2 years ago

Thanks for looking into this!

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.