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?