How to add a raster and a polygon dataset to the same map, each with separate fill aesthetics
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?
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)
- grid.arrange's relative height argument didn't works well in RMarkdown
- How apply a function to every element of a matrix?
- Polar plot in R using ggplot2
- Find total length of time in multiple values of a dataset that overlaps with a given period of time
- Faster than R's `simplify2array`?
- What is the most simple way to label a data.table frequency chart in R?
- How to show dates in a secondary x-axis while the primary x-axis is numerical?
- online compiler for r-markdown document
- Unecessarily large column widths in HTML tables published by R Quarto when the table header contains special characters
- Bootstrapping a vector of results, by group in R
- Speeding Up Monte Carlo Simulations for Weather-Dependent Claims
- Is there an R class which allows to list 2 elements but shows only the first of them by default when called?
- Convert day of week number to weekday name in R
- How do you use "<<-" (scoping assignment) in R?
- Using reactivePoll to monitor and process file changes
- How to find out which graph is most similar to another graph?
- ggplot2: rectangular background behind axis text in circular bar plot
- How to simulate button click using rvest
- CMIP6, automating by month
- Covert dummy variables to single categorical in R?
- Animated reducing bar chart with gganimate
- Create an object class matrix/array from data frame
- How to add a vector of names to `sort_by()` in base R
- Among many subfolders, find CSVs starting with a given string and within ZIPs or not, and merge them while adding their names in a new column with R
- How to add a greek letter to an axis label to a 3D plot in plotly?
- bigrquery + dplyr error when using functions and sourcing from separate R script
- Shiny App Error: /v1/applications/ 400 - Validation Error Execution halted
- Rupee Unicode replacement not working in R
- Mutate multiple variable to create multiple new variables
- How to use the Stata equivalent of a foreach loop in R with mutate()