How to combine SST contour lines and polygons of reef distribution on a world map showing oceans bathymetry?

I would like to create a world map showing the global distribution of coral reefs in combination with the two 20°C isotherms. I have managed to plot them individually, but not in a single graph. The polygons on the coral reef distribution are stored in cd_polygons, but I couldn't add it to the oce plot using map_polygon(). Conversely, the SST data is stored in levitus from the package ocedata, but I was unable to add it to the ggplot using geom_contour(). I prefer the second approach, as I am more familiar with ggplot2 than with base R. Any ideas?

# Load packages
if (!require("pacman")) install.packages("pacman")
#> Lade nötiges Paket: pacman
p_unload("all")
#> The following packages have been unloaded:
#> pacman
pacman::p_load(
  marmap, oce, ocedata, sf, terra, tidyterra, ncdf4, ggplot2, dplyr, tidyr, tibble, rnaturalearth, here, pipeR
)

# Download sea surface data from: https://data.unep-wcmc.org/datasets/1; file size: ca. 200 MB

# working, but bloats the reprex with console messages on the download status:
#usethis::use_zip("https://datadownload-production.s3.us-east-1.amazonaws.com/WCMC008_CoralReefs2021_v4_1.zip", destdir = here()) 

download.file(
  "https://datadownload-production.s3.us-east-1.amazonaws.com/WCMC008_CoralReefs2021_v4_1.zip",
  destfile = "WCMC008_CoralReefs2021_v4_1.zip",
  mode = "wb",
  cacheOK = F
  )

unzip("WCMC008_CoralReefs2021_v4_1.zip")

unlink("WCMC008_CoralReefs2021_v4_1.zip") # deletes .zip file


# Load data
cd_polygons <- (st_read(here("14_001_WCMC008_CoralReefs2021_v4_1/01_Data/WCMC008_CoralReef2021_Py_v4_1.shp"))) %>>%
  select(geometry) %>>%
  st_as_sf()
#> Reading layer `WCMC008_CoralReef2021_Py_v4_1' from data source 
#>   `C:\Users\Marvin Rds\AppData\Local\Temp\RtmpOKlkJI\reprex-3dec40b172a7-brave-mara\14_001_WCMC008_CoralReefs2021_v4_1\01_Data\WCMC008_CoralReef2021_Py_v4_1.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 17504 features and 18 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -179.9999 ymin: -34.29823 xmax: 179.9999 ymax: 32.51482
#> Geodetic CRS:  WGS 84

data(coastlineWorldMedium)
data(levitus, package = "ocedata")

getNOAA.bathy(-180, 180, -90, 90, keep = T) %>>%
  as.topo() %>>%
  (~bathy)
#> Querying NOAA database ...
#> This may take seconds to minutes, depending on grid size
#> Building bathy matrix ...
#> topo object has data as follows.
#>    longitude[1:5400]: -180.00, -179.93, ..., 179.93, 180.00
#>    latitude[1:2700]: -90.000, -89.933, ..., 89.933, 90.000
#>    z, a 5400x2700 array with value 49.1167679 at [1,1] position


# Plotting attempt with oce
mapPlot(coastlineWorldMedium, col = "lightgray", projection = "+proj=robin", drawBox = F)
mapImage(bathy, col = oceColorsGebco, breaks = seq(-4000, 0, 500)) # bathymetric contours, takes very long to plot
mapGrid()
mapContour(levitus[["longitude"]], levitus[["latitude"]], levitus[["SST"]], levels = 20, lwd = 2)
mapPolygon(coastlineWorldMedium, col = "lightgrey")



## Attempt with ggplot, sf, and tidyterra
# get bathymetric data
getNOAA.bathy(-180, 180, -90, 90, keep = T) %>>%
  as.xyz() %>>%
  rename(lon = V1, lat = V2, depth = V3) %>>%
  filter(depth < .01) %>>%
  as_spatraster(xycols = c(1:2), crs = 4326) %>>% # 4326 == EPSG of NOAA bathy
  project("+proj=robin") %>>%
  (~bat_wld)
#> File already exists ; loading 'marmap_coord_-180;-90;180;90_res_4.csv'
#> class       : SpatRaster 
#> dimensions  : 2700, 5400, 1  (nrow, ncol, nlyr)
#> resolution  : 6298.438, 6298.438  (x, y)
#> extent      : -17005830, 17005737, -8502985, 8502798  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#> source(s)   : memory
#> name        :         depth 
#> min value   : -1.011325e+04 
#> max value   :  5.086064e-03

# Import country data
world_map <- ne_countries(scale = "medium", returnclass = "sf")


# Plot using ggplot and sf
world_map %>>%
  st_transform(crs = "+proj=robin") %>>%
  ggplot() +
  geom_spatraster(data = bat_wld, aes(fill = depth)) +
  scale_fill_gradientn(
    colours = c("#5e24d6", "#22496d", "#042f66", "#054780", "#1074a6", "#218eb7", "#48b5d2", "#72d0e1", "#9ddee7", "#c6edec"),
    breaks = c(0, -2500, -5000, -7500),
    guide = "none",
    na.value = NA
  ) +
  geom_sf() +
  geom_sf(data = cd_polygons, aes(col = "red")) +
  coord_sf(crs = "+proj=robin", expand = F) +
  theme(
    panel.background = element_blank(),
    panel.grid.major = element_line(
      linetype = "solid",
      colour = "grey65",
      linewidth = .25
    ),
    panel.grid.minor = element_line( # does not show up, why?
      linetype = "solid",
      colour = "grey65",
      linewidth = .25
    ),
    panel.ontop = T,
    legend.position = "none"
  )
#> <SpatRaster> resampled to 5e+05 cells.

^{Created on 2024-08-28 with reprex v2.1.1}

Solution

This took a minute to figure out! This method combines the data from levitus to the ggplot you created in your question.

This requires the use of a few additional libraries: reshape2 and ggnewscale.

The first step is to manipulate the levitus data. I have a UDF that I know I found from someone else on SO, but I don't remember where (to give them proper credit). Once the levitus data is processed with the UDF, then it's converted to a SpatRaster (like bat_wld in your question).

sst <- function(dt) {   # convert SST into dataframe for raster plotting
  dimnames(dt$SST) <- list(long = dt$longitude, lat = dt$latitude)
  ret <- melt(dt$SST, value.name = "SST")
  cbind(ret, degF = ret$SST * 9/5 + 32)
}

# prep levitus SST for plotting
lsst <- sst(levitus)   # manipulate data into data frame
# convert to SpatRaster
lsstr <- rast(lsst[, c(1:2, 4)], type = "xyz", crs = "+init=epsg:4326")

I call the plot exactly as present in your question and appended the call for the levitus data at the end. When I call the levitus data, I use tidyterra's geom_spatraster_contour.

world_map %>>%
  st_transform(crs = "+proj=robin") %>>%
  ggplot() +
  geom_spatraster(data = bat_wld, aes(fill = depth)) +
  scale_fill_gradientn(
    colours = c("#5e24d6", "#22496d", "#042f66", "#054780", "#1074a6", "#218eb7", "#48b5d2", "#72d0e1", "#9ddee7", "#c6edec"),
    breaks = c(0, -2500, -5000, -7500),
    guide = "none",
    na.value = NA
  ) +
  geom_sf() +
  geom_sf(data = cd_polygons, aes(col = "red")) +
  coord_sf(crs = "+proj=robin", expand = F) +
  theme(
    panel.background = element_blank(),
    panel.grid.major = element_line(
      linetype = "solid",
      colour = "grey65",
      linewidth = .25
    ),
    panel.grid.minor = element_line( # does not show up, why?
      linetype = "solid",
      colour = "grey65",
      linewidth = .25
    ),
    panel.ontop = T,
    legend.position = "none"                   ### additional layer ###
  ) + new_scale_color() +      # add new color scale & levitus data
                               # breaks set to match plot created with `oce` plot
  tidyterra::geom_spatraster_contour(data = lsstr, breaks = seq(70, 200, 20),
                                     alpha = .7, linewidth = 1, color = "black") +
  coord_sf(crs = "+proj=robin", expand = F)   # maintain Robinson projection

Here's everything to make this plot altogether

(easier copy + paste)

pacman::p_load(
  marmap, oce, ocedata, sf, stringr, terra, tidyterra, ncdf4, ggplot2, 
  dplyr, tidyr, tibble, rnaturalearth, here, pipeR, usethis, reshape2, ggnewscale
)

sst <- function(dt) {   # convert SST into dataframe for raster plotting
  dimnames(dt$SST) <- list(long = dt$longitude, lat = dt$latitude)
  ret <- melt(dt$SST, value.name = "SST")
  cbind(ret, degF = ret$SST * 9/5 + 32)
}

#----------- only need to do one time ----------
## Attempt with oce
# Download sea surface data from: https://data.unep-wcmc.org/datasets/1; file size: ca. 200 MB

# working, but bloats the reprex with console messages on the download status:
#use_zip("https://datadownload-production.s3.us-east-1.amazonaws.com/WCMC008_CoralReefs2021_v4_1.zip", destdir = here()) 

# download.file(
#   "https://datadownload-production.s3.us-east-1.amazonaws.com/WCMC008_CoralReefs2021_v4_1.zip",
#   destfile = "WCMC008_CoralReefs2021_v4_1.zip",
#   mode = "wb",
#   cacheOK = F
# )
# 
# unzip("WCMC008_CoralReefs2021_v4_1.zip")
# 
# unlink("WCMC008_CoralReefs2021_v4_1.zip") # deletes .zip file

#------------- unchanged from question --------------
# Load data
cd_polygons <- invisible((st_read(here("14_001_WCMC008_CoralReefs2021_v4_1/01_Data/WCMC008_CoralReef2021_Py_v4_1.shp")))) %>>%
  select(geometry) %>>%
  st_as_sf()
#> Reading layer `WCMC008_CoralReef2021_Py_v4_1' from data source 
#>   `C:\Users\Marvin Rds\AppData\Local\Temp\RtmpOKlkJI\reprex-3dec730f34d7-milky-cub\14_001_WCMC008_CoralReefs2021_v4_1\01_Data\WCMC008_CoralReef2021_Py_v4_1.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 17504 features and 18 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -179.9999 ymin: -34.29823 xmax: 179.9999 ymax: 32.51482
#> Geodetic CRS:  WGS 84

data(levitus, package = "ocedata")

## Attempt with ggplot, sf, and tidyterra
# get bathymetric data
getNOAA.bathy(-180, 180, -90, 90, keep = T) %>>%
  as.xyz() %>>%
  rename(lon = V1, lat = V2, depth = V3) %>>%
  filter(depth < .01) %>>%
  as_spatraster(xycols = c(1:2), crs = 4326) %>>% # 4326 == EPSG of NOAA bathy
  project("+proj=robin") %>>%
  (~bat_wld)

# Import country data
world_map <- ne_countries(scale = "medium", returnclass = "sf") # package: rnaturalearth

#------------------- new content --------------------
# prep levitus SST for plotting
lsst <- sst(levitus)   # manipulate data into data frame
# convert to SpatRaster
lsstr <- rast(lsst[, c(1:2, 4)], type = "xyz", crs = "+init=epsg:4326")

#------------- unchanged from question --------------
world_map %>>%
  st_transform(crs = "+proj=robin") %>>%
  ggplot() +
  geom_spatraster(data = bat_wld, aes(fill = depth)) +
  scale_fill_gradientn(
    colours = c("#5e24d6", "#22496d", "#042f66", "#054780", "#1074a6", "#218eb7", "#48b5d2", "#72d0e1", "#9ddee7", "#c6edec"),
    breaks = c(0, -2500, -5000, -7500),
    guide = "none",
    na.value = NA
  ) +
  geom_sf() +
  geom_sf(data = cd_polygons, aes(col = "red")) +
  coord_sf(crs = "+proj=robin", expand = F) +
  theme(
    panel.background = element_blank(),
    panel.grid.major = element_line(
      linetype = "solid",
      colour = "grey65",
      linewidth = .25
    ),
    panel.grid.minor = element_line( # does not show up, why?
      linetype = "solid",
      colour = "grey65",
      linewidth = .25
    ),
    panel.ontop = T,
    legend.position = "none"             #-------- new layer --------
  ) + new_scale_color() +      # add new color scale & levitus data
                               # breaks set to match plot created with `oce`
  tidyterra::geom_spatraster_contour(data = lsstr, breaks = seq(70, 200, 20),
                                     alpha = .7, linewidth = 1, color = "black") +
  coord_sf(crs = "+proj=robin", expand = F)   # maintain Robinson projection