When I extract lake shapes from OpenStreetMap in R I face a problem with the fact that lakes that contain islands are stored both as multipolygons (with islands cut out) and as regular polygons (one for the lake and one for the island).
The consequence is that if I plot solely the polygons, island polygons are not cut out from the surrounding lake, hence having the colour intended for lakes. (Fig 1.)
If I instead plot the multipolygons, I do not get any of the lakes without islands (fig 2.)
If I knew how to do it, I would use only the polygon data, and cut the islands out of the lakes before plotting. I have a gnawing feeling that this should be pretty simple, however, I haven't been able to find out how to do this operation on a collection of polygons where some are overlapping. (The only solutions I have found are based on polygons from one data set being cut out from polygons in another, which is not how the data received from OpenStreetMap is structured).
Is there any straight forward solution to plot all lakes, where islands are holes cut from their surrounding lakes?
library(osmdata) library(ggplot2) library(sf) x_coords <- c(36.60, 36.66) y_coords <- c(67.99, 68.04) # Define bounding box: bounding_box <- matrix(nrow = 2, ncol=2, byrow = T, data = c(x_coords, y_coords), dimnames = list(c("x","y"), c("min","max"))) # Get waterbodies from OpenStreetMap by a "Overpass query": osm_water_sf <- osmdata::opq(bbox = bounding_box) %>% # Limit query to bounding_box osmdata::add_osm_feature(key = 'natural', value = 'water') %>% # Limit query to waterbodies osmdata::osmdata_sf() # Convert to simple features # plot lakes stored as polygons: ggplot(data=osm_water_sf$osm_polygons) + geom_sf(color="blue", fill="lightblue") + xlim(x_coords) + ylim(y_coords) # plot lakes stored as multi polygons: ggplot(data=osm_water_sf$osm_multipolygons) + geom_sf(color="blue", fill="lightblue") + xlim(x_coords) + ylim(y_coords)