How to add a raster and a polygon dataset to the same map, each with separate fill aesthetics

Question

I am trying to make a map with ocean bathymetry and colored (fill = filled in) polygons by a variable (different color for each Region). I have tried a few different packages and approaches but I am struggling with being able to color/fill by Regions while also having the bathymetry in the map.

I am not sure how to share my data to make reproducible because even if I sample just 1 observation the dput return is extremely lengthy. Can add whatever would be most helpful. In short, it is publicly available data by NOAA.

library(sf)
library(tidyverse)
library(rnaturalearth) #for ne_states()
library(marmap) #for getNOAA.bathy()
library(ggOceanMaps) #for basemap()

Then, for my data I am...

ecomap <- st_read("MYSHAPEFILE.shp") #how to share since dput even for 1 observation is very long
st_crs(ecomap) <- 4326 #set Coordinate Reference System CRS 
ecomap_sf <- st_as_sf(ecomap, wkt = "geom", crs = 4326)

custom_colors <- c("MAB" = "red", 
                   "SNE" = "blue", 
                   "GB" = "green", 
                   "GOM" = "orange")

o

> str(ecomap_sf)
Classes ‘sf’ and 'data.frame':  48 obs. of  7 variables:
 $ Name     : int  1 2 3 4 5 6 7 8 9 10 ...
 $ NumofPoly: int  1 1 1 1 1 1 1 1 1 1 ...
 $ NumofSta : int  1 2 2 1 5 2 2 4 1 3 ...
 $ Region   : chr  "MAB" "MAB" "MAB" "MAB" ...
 $ Area     : num  1540 4496 3190 2510 9620 ...
 $ Type     : chr  "SB" "MS" "IS" "SB" ...
 $ geometry :sfc_POLYGON of length 48; first list element: List of 1
  ..$ : num [1:387, 1:2] -74.8 -74.8 -74.7 -74.7 -74.7 ...
  ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA
  ..- attr(*, "names")= chr [1:6] "Name" "NumofPoly" "NumofSta" "Region" ...


> head(ecomap)
Simple feature collection with 6 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -75.96791 ymin: 35.14219 xmax: -74.28861 ymax: 37.91444
Geodetic CRS:  WGS 84
  Name NumofPoly NumofSta Region     Area Type                       geometry
1    1         1        1    MAB 1540.432   SB POLYGON ((-74.76881 36.5052...
2    2         1        2    MAB 4496.485   MS POLYGON ((-74.81583 36.5052...
3    3         1        2    MAB 3190.414   IS POLYGON ((-75.72364 36.5060...
4    4         1        1    MAB 2510.323   SB POLYGON ((-74.3466 37.56465...
5    5         1        5    MAB 9620.144   MS POLYGON ((-74.63893 37.7118...
6    6         1        2    MAB 2726.496   IS POLYGON ((-75.54484 37.5418...

> unique(ecomap$Region)
[1] "MAB" "SNE" "GB"  "GOM"

I would like to add bathymetry to this map (or similar):

ggplot() +
  borders("world", color = "gray85", fill = "#575151") + 
  geom_sf(data = ecomap_sf1, aes(fill = Region), 
          color = "black", size = 0.5, alpha = 0.5) +    
  scale_fill_manual(values = custom_colors, name = "Region") + 
  coord_sf(xlim = c(-80, -62), ylim = c(34, 48))

I've tried the following:

us <- ne_states(country = "united states of america", returnclass = "sf")

bathy <- getNOAA.bathy(lon1 = -81, lon2 = -62, 
                       lat1 = 24, lat2 = 45, 
                       resolution = 4)
bathy_df <- fortify.bathy(bathy)
bat_xyz <- as.xyz(bathy)


ggplot() +
  geom_sf(data = us) +
  geom_tile(data = bat_xyz, aes(x = V1, y = V2, fill = V3)) +
  geom_sf(data = us)+
  geom_sf(data = ecomap_sf, aes(fill = Region, color = "black"), #the issue is here
          color = "black", 
          size = 0.5, alpha = 0.5) +
  coord_sf(xlim = c(-80, -62), 
           ylim = c(34, 48))

but this gives the following error:

Error in `scale_fill_continuous()`:
! Discrete values supplied to continuous scale.
ℹ Example values: "MAB", "MAB", "MAB", "MAB", and "MAB"

It plots if I use:

aes(color = Region), 
          fill = "transparent"....

coloring the border lines by region (instead of the entire region/polygons) which is good enough for what I need, but I still want to know why filling/coloring in the polygon won't work concurrent with having bathymetry

And I tried using fill=as.factor(Region) and fill=factor(Region) instead, but still get the same error, even when:

> class(ecomap$Region)
[1] "character"

or having it as..

> class(ecomap$Region)
[1] "factor"

I then tried a different approach using ggOceanMaps package

This plots the map (without bathymetry)

basemap(ecomap_sf) + 
  ggspatial::annotation_spatial(ecomap_sf, aes(fill = Region)) + 
  coord_sf(expand = FALSE)

but if I try adding bathymetry=T ...

basemap(ecomap_sf, bathymetry = T) + 
  ggspatial::annotation_spatial(ecomap_sf, aes(fill = Region)) + 
  coord_sf(expand = FALSE)

gives the following error:

Error in `palette()`:
! Insufficient values in manual scale. 12 needed but only 8 provided.

I tried:

basemap(limits = c(-82, -64, 30, 52),
        land.col = "gray21",
        land.border.col = NA,
        bathymetry = T,
        grid.col = NA) +
  geom_sf(data = ecomap_sf, aes(fill = Region), 
          color = "black", size = 0.5, alpha = 0.5) +    
  scale_fill_manual(values = custom_colors, name = "Region") + 
  coord_sf(xlim = c(-80, -62), ylim = c(34, 48))

which plots but "covers" the bathymetry in gray.

I'm guessing the bathymetry colors and the polygon/geometry colors are interacting in a way that does not let me control one over the other..? I don't need a legend for the depths/isobaths. Just a legend for the Regions.

I tried many other things from researching online but cannot figure it out.

This solution from online AI chat kind of works, but it does not seem like I can control the size much because it plots it very small in relation to the limits specified.

base_map1 <- basemap(data = ecomap_sf, 
                     bathymetry = TRUE, bathy.style = "rcb")


base_map1 +
  geom_sf(data = ecomap_sf[ecomap_sf$Region == "GOM", ], 
          fill = custom_colors["GOM"], 
          color = "black", size = 0.5, alpha = 0.2) +
  geom_sf(data = ecomap_sf[ecomap_sf$Region == "GB", ], 
          fill = custom_colors["GB"], 
          color = "black", size = 0.5, alpha = 0.2) +
  geom_sf(data = ecomap_sf[ecomap_sf$Region == "SNE", ], 
          fill = custom_colors["SNE"], 
          color = "black", size = 0.5, alpha = 0.2) +
  geom_sf(data = ecomap_sf[ecomap_sf$Region == "MAB", ], 
          fill = custom_colors["MAB"], 
          color = "black", size = 0.5, alpha = 0.2) +
  coord_sf(xlim = c(-80, -62), ylim = c(34, 48))

To me this seems like an odd approach and not very straightforward

Any suggestions on what is the best way to add bathymetry to a map like this one and be able to fill polygons by a variable?

Answer 1

Your issue arises because your Region values are discrete strings but V3 values are numeric, which are continuous. You can't mix the two. A solution is to use the new_scale_fill() function from the ggnewscale package to reset aes(fill = ) each time it is called. I have combined and adjusted parameters from your various plots, tweaked the bathymetry extent, and edited coord_sf() to avoid whitespace around the bathy_xys data. Modify this to suit your needs:

library(sf)
library(dplyr)
library(rnaturalearth)
library(marmap)
library(ggOceanMaps)
library(ggnewscale)

# Create example sf
ecomap_sf <-  st_as_sfc(st_bbox(c(xmin = -75.96791,
                                  ymin = 35.14219,
                                  xmax = -65.28861,
                                  ymax = 44.91444))) |>
  st_make_grid(n = 4, crs = 4326) |>
  st_as_sf() |>
  mutate(Region = rep(c("MAB", "SNE", "GB",  "GOM"), each = 4))

custom_colors <- c("MAB" = "red", 
                   "SNE" = "blue", 
                   "GB" = "green", 
                   "GOM" = "orange")

us <- ne_states(country = "united states of america", returnclass = "sf")

bathy <- getNOAA.bathy(lon1 = -81, lon2 = -62, 
                       lat1 = 24, lat2 = 47, 
                       resolution = 4)
bathy_df <- fortify.bathy(bathy)
bat_xyz <- as.xyz(bathy)

ggplot() +
  geom_tile(data = bat_xyz, aes(x = V1, y = V2, fill = V3)) +
  geom_sf(data = us) +
  xlab("Longitude (decimal degrees)") + 
  ylab("Latitude (decimal degrees)") +
  guides(fill = guide_legend(title = "Depth(m)")) +
  new_scale_fill() +
  geom_sf(data = ecomap_sf, aes(fill = Region), 
          color = "black",
          size = 0.5,
          alpha = 0.5) +    
  scale_fill_manual(values = custom_colors) +
  coord_sf(xlim = c(-80, -62), 
           ylim = c(34, 46),
           expand = FALSE)