Binding several sf together only plots first layer

Question

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

Allan Cameron · Accepted Answer · 2022-05-17T12:24:07.313

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)

Thanks Allan, that works well. I came across this "..." statements a few times but never dared trying. A whole new world opens up to me :) The only thing that does not work with your solution is that I'd like to be able to show a legend with the different depths and corresponding colors, which works with my initial try but not with your solution. Do you see how we could solve this too, now that each layer is plotted separately? — Fred-LM, May 17 '22 at 12:13
@Fred-LM I'd probably switch to `scale_fill_manual` in that case. See my update. — Allan Cameron, May 17 '22 at 12:24

Binding several sf together only plots first layer

1 Answers1