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.)

enter image description here

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:

  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"),
  # 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)

2 Answers 2


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,
                     .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) %>% 

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

lakes with islands cut out

  • Thanks! I was trying the st_difference as well, getting frustrated when it did not work as I thought it would. I will flag your answer as the "accepted" answer since you also shed some light on the cutting problem. (Although I suppose Elio Diaz' method is more computationally efficient.)
    – Smerla
    Commented Jul 14, 2021 at 9:11
  • 1
    @Smerla Thanks for your kind words! I did not really mean to compete with Elio Diaz, but to suggest a possible alternative. The root issue here is a messy geometry, and messy geometries have a kind of an Anna Karenina issue - all clean datasets are pure, but each messy datasets is messy in a special way, and it pays to have a variety of tools at your disposal :) Commented Jul 14, 2021 at 10:40
  • 1
    Indeed, I believe that the diversity in solutions is one of the real strong aspects of Stackexchange. However, the "competition" part is an involuntary effect of the voting system that only allows one accepted answer
    – Smerla
    Commented Jul 14, 2021 at 11:45

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

enter image description here

  • 1
    Thanks! That's a great answer. I started explore some filtering solution on my own, but did not manage to come up with anything good. I have troubles between choosing to mark your or Jindra Lackos as the "accepted answer" since I like both.
    – Smerla
    Commented Jul 14, 2021 at 9:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.