Plot filled areas for sea/ocean and land mass based on {osmdata} using {ggplot2}

The reprex below shows how I would like to create a map via {osmdata} and {ggplot2} that has sea/ocean in it. I want to color-fill the land and/or sea area. However, it seems unexpectedly difficult to do so. This blog post even claims that it cannot be done.

This vignette of {osmplotr} seems to have to the solution: "Because OpenStreetMap represents coastline as line objects, all coastline data is contained within the $osm_lines object. The osm_line2poly() function can then convert these lines to polygons which can be used to plot filled areas.". Yet, just as in this similar StackOverflow question, the function throws an error as can be seen at the bottom of the reprex. I also found here that the {tigris} package can provide the necessary polygon data - but only for the US.

So how can I get this to work?

library(osmdata)
library(osmplotr)
library(sf)
library(tidyverse)

# define example bbox
bb <- tribble(
  ~xy,  ~min,  ~max,
  "x", 12.00, 12.18,
  "y", 54.08, 54.20
) %>% column_to_rownames("xy") %>% as.matrix()

# get "water"
water <- opq(bb) %>%
  add_osm_feature(key = "natural", value = "water") %>%
  osmdata_sf()

# get "coastline"
coast <- opq(bb) %>%
  add_osm_feature(key = "natural", value = "coastline") %>%
  osmdata_sf()

# ggplot
ggplot() +
  geom_sf(
    data = water$osm_multipolygons,
    fill = "navy",
    color = NA
  ) + 
  geom_sf(
    data = coast$osm_lines,
    fill = "navy",
    color = "blue"
  )

# trying osm_line2poly()
osmplotr::osm_line2poly(coast$osm_lines, bb)
#> Error in FUN(X[[i]], ...): unbenutztes Argument (V = c(3, 1, 6, 7, 2, NA, 5))

^{Created on 2022-09-23 with reprex v2.0.2}

Solution

Thanks to @JindraLacko, I was able to make my reprex work. Basically, we create a rectangle/polygon which is the size of our bbox and then split it via the coastline.

library(lwgeom)
library(osmdata)
library(osmplotr)
library(sf)
library(tidyverse)

### define example bbox
lon_min <- 12.00 # xmin
lon_max <- 12.18 # xmax
lat_min <- 54.08 # ymin
lat_max <- 54.20 # ymax

bb <- get_bbox(c(lon_min, lat_min, lon_max, lat_max))

### get "water" that is not sea as polygons
water <- opq(bb) %>%
  add_osm_feature(key = "natural", value = "water") %>%
  osmdata_sf()

### get sea & land as polygons
# 1. get coastline (as line)
coast <- opq(bb) %>%
  add_osm_feature(key = "natural", value = "coastline") %>%
  osmdata_sf()

# 2. get overall rectangle for bbox
bb_rect <- data.frame(
  lat = c(lat_min, lat_max),
  lon = c(lon_min, lon_max)
) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  st_bbox() %>% 
  st_as_sfc()

# 3. split overall rectangle for bbox via coastline
bb_rect_split <- bb_rect %>% 
  st_split(coast$osm_lines) %>% 
  st_collection_extract("POLYGON")

# 4. extract splitted rectangle parts
land <- bb_rect_split[1]
sea  <- bb_rect_split[2]

### ggplot
ggplot() +
  geom_sf(
    data = land,
    fill = "bisque",
    color = NA
  ) + 
  geom_sf(
    data = sea,
    fill = "navy",
    color = NA
  ) + 
  geom_sf(
    data = water$osm_multipolygons,
    fill = "navy",
    color = NA
  )

^{Created on 2022-09-26 with reprex v2.0.2}