Lake shapes from OpenStreetMap (OSM) in R. Polygon vs multipolygon issues

Question

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?

Reproducible example:

  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)

Answer 1

As an alternative to the approach taken by elio-diaz I propose the following; the logic is built on on

1. identifying the "lakes on lakes" as polygons touching the multipolygon (the Fig.2 from original post) and
2. erasing the island area from the polygon layer via rmapshaper::ms_erase()

It is interesting that for some unknown (to me) reason sf::st_difference() does not give the expected result, when in theory it should be equivalent to the mapshaper approach.

# identify islands as polygons touching the multipolygon
islands <- st_filter(osm_water_sf$osm_polygons,
                     osm_water_sf$osm_multipolygons,
                     .predicate = st_touches) 

# filter out the islands / via osm_id + erase them from underlying layer
w2 <- osm_water_sf$osm_polygons %>% 
  dplyr::filter(!osm_id %in% islands$osm_id) %>% 
  rmapshaper::ms_erase(islands)

# plot lakes stored as polygons:
ggplot(data=w2) +
  geom_sf(color="blue", fill="lightblue") +
  xlim(x_coords) + ylim(y_coords)

Answer 2

After taking a look at the data, and unless an osm expert points out a better filtering option, I suggest to filter out those polygons that "have" the islands and the islands as well, although they are not stored as polygons with rings (holes), in this case the big lake and the islands are three separate polygons. We may take them out using st_filter

# this returns the lake polygon
p_contains = osm_water_sf$osm_polygons %>% 
  st_filter(osm_water_sf$osm_polygons, .predicate = st_contains_properly) 

# this returns island polygons
# had to use this st_relate() pattern, st_within is TRUE for polygons overlapping themselves

p_contained = osm_water_sf$osm_polygons %>% 
  st_filter(osm_water_sf$osm_polygons, .predicate = st_relate, pattern = "2FF1FF212") 

p_drop = rbind(p_contains, p_contained) %>% as.data.frame() %>% select(osm_id) 

osm_water_sf$osm_polygons %>% 
  filter(!osm_id %in% p_drop$osm_id) %>%
  bind_rows(osm_water_sf$osm_multipolygons) %>% 
  ggplot() + geom_sf()

On st_relate: https://r-spatial.github.io/sf/reference/st_relate.html