Skip to contents

Extracting 3D Environmental Data for a Single Species

Source: vignettes/articles/woa-environmental-extraction-single-species.Rmd
woa-environmental-extraction-single-species.Rmd

Extracting 3D Environmental Data for Sharks and Rays Species

This vignette demonstrates the World Ocean Atlas (WOA) environmental extraction analysis, which summarises temperature and dissolved oxygen conditions across the 3D habitat of shark and ray species. Here we implement this workflow using ocean3d package functions.

Overview

The goal is to produce a single summary table with one row per species, containing statistics (min, max, mean, range) for temperature and dissolved oxygen across each species’ geographic range and depth limits. The workflow has five steps:

  1. Fetch species depth limits from the IUCN Red List API
  2. Load species range polygons
  3. Prepare WOA rasters and the study grid (download, annual means, align to grid)
  4. Extract and summarise environmental conditions per species

This vignette is a scoped-down version of the full WOA environmental extraction workflow that runs for one example species instead of every shark and ray.

The example species is the Banded Wobbegong (Orectolobus halei).

Setup 

library(ocean3d)
library(dplyr)
library(sf)
library(stringr)
library(terra)
library(here)
library(ggplot2)
library(tidyr)
library(tidyterra)

# The single example species.
sp_name <- "Orectolobus halei"

Step 1: Fetch depth limits for the example species

Fetching by species_names returns just this one species’ assessment, so there is no 30-minute group download.

api_key <- Sys.getenv("IUCN_REDLIST_KEY")

# Use IUCN Red List API to get the species depth ranges
depth_table <- fetch_species_assessments(
  api_key = api_key,
  species_names = sp_name
)

# Fill missing depths
depth_table_complete <- depth_table |>
  mutate(fill_missing_depths(upper_depth_limit, lower_depth_limit, genus_name)) |>
  filter(!is.na(upper_depth), !is.na(lower_depth))

# Get the Orectolobus halei depth range
sp_depths <- depth_table_complete[
  depth_table_complete$scientific_name == sp_name, ]
sp_id <- sp_depths$sis_id

Step 2: Load the species range polygon

The IUCN Red List shapefiles have to be downloaded yourself, from this link: https://www.iucnredlist.org/resources/spatial-data-download

Only the one species is read out of the IUCN shapefile via SQL filtering.

# Download the IUCN Red List shapefiles from https://www.iucnredlist.org/resources/spatial-data-download 
# Modify the path to your local save location for the shapefile
iucn_shp <- here("EXAMPLE_DIRECTORY/SHARKS_RAYS_CHIMAERAS/SHARKS_RAYS_CHIMAERAS.shp")
layer_name <- sf::st_layers(iucn_shp)$name[1]

# Read in just Orectolobus halei
sp_range <- sf::st_read(
  iucn_shp,
  query = str_glue(
    "SELECT * FROM \"{layer_name}\" ",
    "WHERE id_no = {sp_id} ",
    "AND presence IN (1) AND origin IN (1, 2) AND seasonal IN (1, 2, 3)"
  ),
  quiet = TRUE
)

Step 3: Prepare the WOA temperature raster and study grid 

# Annual mean temperature (0.25 deg) is both the environmental raster and the
# canonical study grid.
t_annual_path <- woa_download("temperature", period = "annual", resolution = "0.25")
t_annual <- woa_load_nc(t_annual_path, field = "an")

# Crop to the species range (plus a ~1 degree margin) before any heavy work.
# WOA ships globally, but the bathymetry projection, envelope construction, and
# per-depth masking all cost in proportion to grid size. A circumglobal range
# whose bbox spans the antimeridian makes this a no-op.
range_ext <- terra::intersect(
  terra::ext(terra::vect(sp_range)) + 1,
  terra::ext(-180, 180, -90, 90)
)

# as_voxel() re-wraps after the crop: terra operations propagate the class but
# do not re-run its validity rules, so wrapping is how the depth axis is
# re-checked. It is idempotent, so this is cheap.
t_annual <- as_voxel(terra::crop(t_annual, range_ext))

# Derive the study grid from the cropped raster so the envelope and the voxel
# stay on the same grid.
study_grid <- t_annual[[1]]
# Bathymetry -> positive depth, aligned to the study grid.
# Project before negating: negating first forces an eager pass over the full
# global 15 arc-second grid, which terra otherwise avoids by reading only the
# window the target grid needs.

# Download the GEBCO bathymetry data from https://www.gebco.net/data-products/gridded-bathymetry-data
bathy <- load_gebco_bathymetry("EXAMPLE_DIRECTORY/gebco_2025_sub_ice_topo/GEBCO_2025_sub_ice.nc")
seafloor <- (terra::project(bathy, study_grid) * -1) |>
  terra::clamp(lower = 0)

Step 4: Build the species’ depth envelope

vect_to_envelope() builds the per-cell depth window. The seafloor goes in as one more depth_max constraint, so the species’ nominal lower limit is clamped to the bed wherever the bed is shallower; per cell the shallowest depth_max wins.

range_rast <- vect_to_envelope(
  polygon = sp_range,
  template = study_grid,
  depth_min = sp_depths$upper_depth,
  depth_max = list(sp_depths$lower_depth, seafloor)
)

# Mask temperature to the cells + depths the species actually occupies.
# mask() is depth-aware for an envelope: it puts the range on the temperature
# stack's own depth levels and masks layer by layer. The grid is already
# scoped to the range, and cells outside the polygon are NA, so no further
# cropping is needed.
masked <- mask(t_annual, range_rast)

Figure 1: Animated GIF across all depths

Each frame is one depth layer, from the surface down to the deepest layer the species occupies. The colour scale is locked to the global temperature range so colours are comparable across frames.

long_df <- terra::as.data.frame(masked, xy = TRUE, na.rm = FALSE) |>
  pivot_longer(cols = -c(x, y), names_to = "layer", values_to = "value") |>
  mutate(depth = depths(layer)) |>
  filter(!is.na(value))

# Lock the legend to the global min/max so colour mapping is identical per frame.
v_limits <- range(long_df$value, na.rm = TRUE)

# Factor with numerically-sorted levels so frames play 0 m -> deepest.
long_df$depth <- factor(long_df$depth, levels = sort(unique(long_df$depth)))

p_anim <- ggplot(long_df, aes(x = x, y = y, fill = value)) +
  geom_raster() +
  geom_sf(data = sp_range, fill = NA, colour = "black",
          linewidth = 0.3, inherit.aes = FALSE) +
  scale_fill_viridis_c(name = "Temp (°C)", limits = v_limits,
                       option = "turbo",
                       na.value = "transparent") +
  coord_sf() +
  labs(subtitle = paste0(sp_name, " — depth: {current_frame} m")) +
  theme_minimal() +
  theme(
    plot.subtitle = element_text(face = "italic"),
    axis.title = element_blank(),
    axis.text = element_blank(),
    panel.grid.major = element_blank()
  ) +
  transition_manual(depth)

anim <- animate(
  p_anim,
  nframes = nlevels(long_df$depth),
  fps = 4, width = 700, height = 400,
  renderer = gifski_renderer()
)

anim_save(here("figures/single_species_depth_animation.gif"), anim)
anim

Figure 1. Annual mean sea temperature (World Ocean Atlas) within the geographic range of the Banded Wobbegong (Orectolobus halei), animated through every standard depth layer the species occupies. Each frame maps one depth, from the sea surface (0 m) down to the deepest layer within the species’ bathymetry-clamped depth limits; the black outline is the IUCN range polygon. The colour scale is fixed across all frames, so warming or cooling with depth is directly comparable frame to frame.

Figure 2: Three-panel static figure (surface / halfway / deepest)

We pick three depths from the layers the species actually occupies: the surface (0 m), the standard WOA depth closest to the halfway point, and the deepest layer with data inside the range.

# Depths of each masked layer, and which layers hold any data for this species.
layer_depths <- depths(masked)
has_data <- vapply(seq_len(terra::nlyr(masked)), function(i) {
  !all(is.na(terra::values(masked[[i]])))
}, logical(1))
occupied_depths <- sort(layer_depths[has_data])

surface_depth <- 0
deepest_depth <- max(occupied_depths)
# Standard depth closest to the midpoint between surface and deepest occupied.
halfway_target <- (surface_depth + deepest_depth) / 2
halfway_depth <- occupied_depths[which.min(abs(occupied_depths - halfway_target))]

panel_depths <- c(surface_depth, halfway_depth, deepest_depth)

The plot_range_at_depth() helper below is defined in this article rather than exported by the package — it is a thin ggplot2 wrapper that draws one map per depth. patchwork stacks them vertically.

# Map view of a species range with environmental values at the depth layer
# nearest `depth`.
plot_range_at_depth <- function(species_range, depth, rast_3d) {
  layer_depths <- depths(rast_3d)
  idx <- which.min(abs(layer_depths - depth))
  layer <- rast_3d[[idx]]
  actual_depth <- layer_depths[idx]

  range_vect <- if (inherits(species_range, "sf")) {
    terra::vect(species_range)
  } else {
    species_range
  }
  if (terra::crs(range_vect) != terra::crs(layer) &&
      terra::crs(range_vect) != "") {
    range_vect <- terra::project(range_vect, terra::crs(layer))
  }

  cropped <- terra::crop(layer, range_vect, mask = TRUE, touches = TRUE)
  range_sf <- if (inherits(species_range, "sf")) {
    species_range
  } else {
    sf::st_as_sf(range_vect)
  }

  ggplot() +
    geom_spatraster(data = cropped) +
    geom_sf(data = range_sf, fill = NA, colour = "black", linewidth = 0.3) +
    scale_fill_viridis_c(name = NULL, na.value = "transparent") +
    coord_sf() +
    labs(subtitle = paste0("Depth: ", actual_depth, " m")) +
    theme_minimal() +
    theme(axis.title = element_blank())
}
# Shared colour limits across the three panels. Without explicit `limits` each
# panel rescales to its own range, so the warm surface and the cold deepest
# layer would render in the same colours and the panels would not be
# comparable. Mirror plot_range_at_depth()'s own layer pick and crop so the
# limits cover exactly the values it draws.
panel_layers <- terra::rast(lapply(panel_depths, function(d) {
  t_annual[[which.min(abs(layer_depths - d))]]
}))
fill_limits <- range(
  terra::values(terra::crop(panel_layers, terra::vect(sp_range),
                            mask = TRUE, touches = TRUE)),
  na.rm = TRUE
)

panels <- lapply(seq_along(panel_depths), function(i) {
  d <- panel_depths[i]
  plot_range_at_depth(species_range = sp_range, depth = d, rast_3d = t_annual) +
    scale_fill_viridis_c(name = "Temp (°C)", na.value = "transparent",
                         option = "turbo", limits = fill_limits) +
    # Fold the panel letter into the subtitle rather than using patchwork's
    # tag_levels, which renders the tag on its own line above this text.
    labs(subtitle = paste0(LETTERS[i], ". Depth: ", d, " m"))
})

fig2 <- (panels[[1]] / panels[[2]] / panels[[3]]) +
  plot_annotation(
    title = sp_name,
    theme = theme(plot.title = element_text(face = "italic"))
  ) +
  plot_layout(guides = "collect")

ggsave(here("figures/single_species_three_panel.png"), fig2,
       width = 6, height = 12, dpi = 300)
fig2

Figure 2. Annual mean sea temperature (World Ocean Atlas) within the range of the Banded Wobbegong (Orectolobus halei) at three depths, stacked from shallow (top) to deep (bottom): (A) the sea surface (0 m), (B) the standard WOA depth nearest the midpoint of the species’ occupied depth range, and (C) the deepest standard depth at which the species occurs (its depth limit clamped to the seafloor). Each panel is masked to the IUCN range polygon and shares a common temperature colour scale, illustrating how the thermal environment available to the species cools with depth.