Как совместить контурные линии ТПО и полигоны распределения рифов на карте мира, показывающей батиметрию океанов?

Это заняло минуту, чтобы понять! Этот метод объединяет данные от levitus до ggplot, которые вы создали в своем вопросе.

Это требует использования нескольких дополнительных библиотек: reshape2 и ggnewscale.

Первый шаг — манипулировать данными levitus. У меня есть UDF, который, как я знаю, я нашел у кого-то еще на SO, но не помню, где (чтобы отдать им должное). Как только данные levitus обрабатываются с помощью UDF, они преобразуются в SpatRaster (например, bat_wld в вашем вопросе).

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")

Я называю сюжет точно так же, как в вашем вопросе, и добавляю в конце вызов данных levitus. Когда я вызываю данные levitus, я использую tidyterrageom_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

Вот все, чтобы сделать этот сюжет вообще

(проще скопировать + вставить)

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