Binding several sf together only plots first layer

I'm trying to create a script to generate generic maps with bathymetry, and I'm struggling to get it to work. My issue is that depending on the map I want to make, I will call different bathymetric layers, but when I bind them together, only the first one is plotted. I don't get where the issue comes from. Any idea?

I created a function to load my bathymetric shapefiles from the 'rnaturalearth' package:

load_shapefile <- function(file, url, shapefile){ 
  if (!file.exists(file.path(tempdir(), file))) {
    url <- url
    download.file(url, file.path(tempdir(), file))
    unzip(file.path(tempdir(), file), exdir = tempdir())
  }
  st_read(dsn = tempdir(), layer = shapefile,
          quiet = TRUE)
}

I then load a few shapefiles:

bathy_200 <- load_shapefile("ne_10m_bathymetry_K_200.zip",
                            "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_bathymetry_K_200.zip",
                            "ne_10m_bathymetry_K_200")
bathy_1000 <- load_shapefile("ne_10m_bathymetry_J_1000.zip",
                             "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_bathymetry_J_1000.zip",
                             "ne_10m_bathymetry_J_1000")
bathy_2000 <- load_shapefile("ne_10m_bathymetry_I_2000.zip",
                             "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_bathymetry_I_2000.zip",
                             "ne_10m_bathymetry_I_2000")

Let's say I only want to plot the 200 and 2000 isobaths, but still want to create a generic function that would also easily allow me to plot 1000 and 2000, or 200 and 1000 instead:

plot_map <- function(bathy){
ggplot() +
  geom_sf(data = bathy %>%
            left_join(tibble(depth = c(200, 1000, 2000), 
                             fill = c("#E2EFF6", "#B7D7EA", "#8DBEDC")),
                      by = "depth"),
          aes(fill = fill),
          color = NA) +
  coord_sf(xlim = c(-20, 20),
           ylim = c(-20, 20),
           expand = FALSE) +
  scale_fill_identity()
}

When I plot my data, only the 200m shapefile is plotted...

plot_map(bind_rows(bathy_200, bathy_2000))

What dit I miss? I've tried several solutions, e.g. arranging data etc. but the outcome remains the same...

I think your main problem is with overplotting, which you can only really see by setting strong fill colors and an alpha value:

plot_map <- function(bathy){
ggplot() +
  geom_sf(data = bathy %>%
            left_join(tibble(depth = c(200, 1000, 2000), 
                             fill = c("red", "green", "blue")),
                      by = "depth"),
          aes(fill = fill),
          color = NA, alpha = 0.5) +
  coord_sf(xlim = c(-20, 20),
           ylim = c(-20, 20),
           expand = FALSE) +
  scale_fill_identity()
}

plot_map(bind_rows(bathy_200, bathy_2000))

This problem would disappear if you could plot each data frame in a separate layer. To control this without having to specify which layer gets plotted first, I would write the function like this:

plot_map <- function(...){
  
  dfs <- list(...)
  
  dfs <- dfs[order(sapply(dfs, function(x) mean(x$depth)))]
  
  cols <- c("200" = "#E2EFF6", "1000" = "#B7D7EA", "2000" = "#8DBEDC")
  
  ggplot() +
    lapply(dfs, function(df) {
      geom_sf(data = df, aes(fill = factor(depth)), color = NA)
    }) +
    coord_sf(xlim = c(-20, 20),
             ylim = c(-20, 20),
             expand = FALSE) +
    scale_fill_manual(values = cols) +
    theme_bw() +
    labs(fill = 'depth') +
    theme(panel.background = element_rect(fill = '#A0A978'),
          panel.grid = element_line(colour =  '#C0C0C080'))
}

This allows you to pass in the data frames separately instead of binding them. The shallowest layer will always be plotted first, with deeper layers on top

plot_map(bathy_200, bathy_1000, bathy_2000)