I'm extracting data from OpenStreetMap using osmdata
, and from the query I get a set of points, lines, polygons, and multipolygons. Somehow, when trying to plot multipolygons, leaflet and tmap do not find them, while using plot(st_geometry(...))
works well.
Any ideas as to what is happening?
library(osmdata)
library(sf)
library(tmap)
nkpg_schools <-
opq(bbox = "Norrköping") %>%
add_osm_feature(key = "amenity", value = "school") %>%
osmdata_sf()
plot(st_geometry(nkpg_schools$osm_multipolygons)) #Works well
tmap_mode("view")
tm_basemap("Esri.WorldGrayCanvas") +
tm_shape(nkpg_schools$osm_multipolygons) +
tm_polygons()
#Does not work
This is an interesting problem! My first hunch was that the Norrköping schools multipolygon object was spatially invalid (which does happen with the OSM objects from time to time, unfortunately).
This seems not to be the case, see sf::st_is_valid()
test results.
But! curiouser and curiouser: creating an intermediary object of schools made valid via sf::st_make_valid()
causes the schools to behave themselves.
I confess I can not fully explain what's going there, but if it solves your problem - who are we to judge?
library(osmdata)
library(sf)
library(tmap)
nkpg_schools <-
opq(bbox = "Norrköping") %>%
add_osm_feature(key = "amenity", value = "school") %>%
osmdata_sf()
st_is_valid(nkpg_schools$osm_multipolygons)
#[1] TRUE TRUE TRUE TRUE
intmd <- nkpg_schools$osm_multipolygons %>%
st_make_valid()
tmap_mode("view")
tm_shape(intmd) +
tm_polygons(col = "red") +
tm_basemap("Esri.WorldGrayCanvas")