Copernicus for Sanctuaries

Code 
librarian::shelf(
  dplyr,
  DT,
  fs,
  ggplot2,
  glue,
  here,
  jsonlite,
  leaflet,
  listviewer,
  mapview,
  ncdf4,
  purrr,
  readr,
  reticulate,
  sf,
  stringr,
  terra,
  tidyr,
  quiet = T
)
options(readr.show_col_types = F)

sanctuaries_geo <- here("data/sanctuaries.geojson")
# copernicusmarine username and password
user <- "bbest1"
pass_txt <- "~/My Drive/private/data.marine.copernicus.eu_bbest1-password.txt"
dir_cm <- here("data/copernicus")
results_csv <- here("data/copernicus.csv")

Sanctuaries

We are extracting Copernicus Marine data from across National Marine Sanctuaries, borrowing polygons from the noaa-onms/climate-dashboard (see Github code).

Code 
if (!file.exists(sanctuaries_geo)) {
  # source: https://github.com/noaa-onms/climate-dashboard/blob/main/data/sanctuaries.rds
  sanctuaries_rds <- here("../climate-dashboard/data/sanctuaries.rds")
  readRDS(sanctuaries_rds) |>
    filter(nms != "TBNMS") |> # exclude Great Lakes
    write_sf(sanctuaries_geo, delete_dsn = T)
}

sanctuaries <- read_sf(sanctuaries_geo)
sanctuaries |>
  st_drop_geometry() |>
  datatable()
Code 
mapView(sanctuaries)

Setup for copernicusmarine in R

Code 
# do once: create virtual enviroment and install copernicusmarine Python module
# virtualenv_create(envname = "CopernicusMarine")
# virtualenv_install(envname = "CopernicusMarine", packages = c("copernicusmarine"))

# TODO: check for CopernicusMarine env with copernicusmarine Python module

# use virtualenv and reticulate::import copernicusmarine Python module
use_virtualenv(virtualenv = "CopernicusMarine", required = TRUE)
cmt <- import("copernicusmarine")

# login
pass <- readLines(pass_txt)
stopifnot(
  cmt$login(user, pass, force_overwrite = T)
) # py_help(cmt$login)
# writes to ~/.copernicusmarine/.copernicusmarine-credentials
Code 
get_dataset_description <- function(dataset_id) {
  cmt$describe(
    dataset_id           = dataset_id,
    disable_progress_bar = T)$json() |>
    fromJSON()
}

get_dataset_times <- function(dataset_description) {
  # dataset_description = ds_times
  services <- pluck(
    dataset_description, "products", "datasets", 1,
    "versions", 1, "parts", 1, "services", 1)
  service_names <- pluck(services, "service_name")
  i_ts <- which(service_names == "arco-time-series")

  coordinates <- pluck(services, "variables", i_ts, "coordinates")
  coord_ids   <- pluck(coordinates, 1, "coordinate_id")
  i_time      <- which(coord_ids == "time")

  time_vals <- pluck(coordinates, 1, "values", i_time) |>
    as.numeric() %>% {./1000}
  time_min <- min(time_vals) |>
    as.POSIXct(origin = "1970-01-01", tz = "UTC")
  time_max <- max(time_vals) |>
    as.POSIXct(origin = "1970-01-01", tz = "UTC")

  stopifnot(
    "get_dataset_times(): time_min is NA" = !is.na(time_min),
    "get_dataset_times(): time_max is NA" = !is.na(time_max))

  list(min = time_min, max = time_max)
}

get_last_date_in_csv <- function(results_csv) {
  if (!file.exists(results_csv))
    return(tibble(
      nms = character(), dataset = character(),
      last_time = as.Date(character())))
  read_csv(results_csv) |>
    mutate(time = as.Date(time)) |>
    group_by(nms, dataset) |>
    summarize(last_time = max(time, na.rm = TRUE), .groups = "drop")
}

Extract Global Ocean Physics Reanalysis

Product > Dataset > Variable

Variables, keep:

  • thetao: Sea water potential temperature [°C]
  • bottomT: Sea water potential temperature at sea floor [°C]
  • so: Sea water salinity [/ 103]
  • mlotst: Ocean mixed layer thickness defined by sigma theta [m]

Variables, skip:

  • siconc: Sea ice area fraction
  • sithick: Sea ice thickness [m]
  • usi: Eastward sea ice velocity [m/s]
  • vsi: Northward sea ice velocity [m/s]
  • uo: Eastward sea water velocity [m/s]
  • vo: Northward sea water velocity [m/s]
  • zos: Sea surface height above geoid [m]

Extract Global Ocean Colour

Setup dataset parameters

Code 
datasets <- list(
  # Global Ocean Physics Reanalysis - monthly
  # https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/download?dataset=cmems_mod_glo_phy_my_0.083deg_P1M-m
  "cmems_mod_glo_phy_my_0.083deg_P1M-m" = list(
    vars      = list("thetao", "bottomT", "so", "mlotst"),
    depth_rng = c(0, 0.5)),
  # Plankton, multi-sensor (4 km), monthly
  "cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M" = list(
    vars = list(
      "CHL",
      "CHL_uncertainty",
      "DIATO",
      "DIATO_uncertainty",
      "DINO",
      "DINO_uncertainty",
      "GREEN",
      "GREEN_uncertainty",
      "HAPTO",
      "HAPTO_uncertainty",
      "MICRO",
      "MICRO_uncertainty",
      "NANO",
      "NANO_uncertainty",
      "PICO",
      "PICO_uncertainty",
      "PROCHLO",
      "PROCHLO_uncertainty",
      "PROKAR",
      "PROKAR_uncertainty",
      "flags"),
    depth_rng = c(0, 0.5)),
  # Primary production, multi-sensor (4 km), monthly
  "cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M" = list(
    vars      = list("PP", "PP_uncertainty", "flags"),
    depth_rng = c(0, 0.5))
)

# query Copernicus API to set date_rng (min & max) per dataset
fmt_dt <- \(x) format(x, "%Y-%m-%dT00:00:00")
for (id in names(datasets)) {  # id = "cmems_mod_glo_phy_my_0.083deg_P1M-m"
  ds_times <- get_dataset_description(id) |> get_dataset_times()
  datasets[[id]]$date_rng <- c(fmt_dt(ds_times$min), fmt_dt(ds_times$max))
  message(glue("  {id}: {datasets[[id]]$date_rng[1]} -> {datasets[[id]]$date_rng[2]}"))
}
  cmems_mod_glo_phy_my_0.083deg_P1M-m: 1993-01-01T00:00:00 -> 2026-02-01T00:00:00
  cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M: 1997-09-01T00:00:00 -> 2026-03-01T00:00:00
  cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M: 1997-09-01T00:00:00 -> 2026-03-01T00:00:00
Code 
# summarize last fetched date per dataset from existing results
last_dates <- get_last_date_in_csv(results_csv)
for (id in names(datasets)) {
  last <- last_dates |>
    filter(dataset == id) |>
    pull(last_time)
  datasets[[id]]$date_last_fetched <- if (length(last) > 0)
    as.character(min(last)) else NA
}

datasets |>
  toJSON(auto_unbox = T, pretty = T) |>
  jsonedit() # view in browser

Dataset Variables

Code 
vars_csv <- here("data/copernicus_variables.csv")

get_var_coord_range <- function(coord_df, coord_name) {
  if (!is.data.frame(coord_df)) return(c(NA_real_, NA_real_))
  i <- which(coord_df$coordinate_id == coord_name)
  if (length(i) == 0) return(c(NA_real_, NA_real_))
  vals <- as.numeric(coord_df$values[[i]])
  c(min(vals, na.rm = T), max(vals, na.rm = T))
}

get_dataset_vars_info <- function(dataset_id) {
  desc     <- get_dataset_description(dataset_id)
  product  <- pluck(desc, "products") |> as_tibble() |> slice(1)
  services <- pluck(
    desc, "products", "datasets", 1,
    "versions", 1, "parts", 1, "services", 1)
  i_ts     <- which(services$service_name == "arco-time-series")
  vars     <- pluck(services, "variables", i_ts)
  times    <- get_dataset_times(desc)

  tibble(
    dataset_id    = dataset_id,
    product_id    = product$product_id,
    product_title = product$title,
    time_min      = as.character(times$min),
    time_max      = as.character(times$max),
    variable      = vars$short_name,
    standard_name = vars$standard_name,
    units         = vars$units,
    depth_min     = map_dbl(vars$coordinates, \(c) get_var_coord_range(c, "depth")[1]),
    depth_max     = map_dbl(vars$coordinates, \(c) get_var_coord_range(c, "depth")[2]),
    longitude_min = map_dbl(vars$coordinates, \(c) get_var_coord_range(c, "longitude")[1]),
    longitude_max = map_dbl(vars$coordinates, \(c) get_var_coord_range(c, "longitude")[2]),
    latitude_min  = map_dbl(vars$coordinates, \(c) get_var_coord_range(c, "latitude")[1]),
    latitude_max  = map_dbl(vars$coordinates, \(c) get_var_coord_range(c, "latitude")[2])
  )
}

vars_info <- map_dfr(names(datasets), get_dataset_vars_info)
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Warning in min(vals, na.rm = T): no non-missing arguments to min; returning Inf
Warning in max(vals, na.rm = T): no non-missing arguments to max; returning
-Inf
Code 
write_csv(vars_info, vars_csv)

datatable(vars_info, options = list(scrollX = TRUE))

Iterate over sanctuaries & datasets to subset

Code 
last_dates <- get_last_date_in_csv(results_csv)

for (i in 1:nrow(sanctuaries)) {
  # i = 1

  d_s <- slice(sanctuaries, i)
  bb  <- d_s |>
    st_buffer(9.25 * 1000) |> # buffer by a pixel width 1/12° ~ 9.25 km at equator
    st_bbox()

  dir_out <- here(glue("{dir_cm}/{d_s$nms}"))
  dir.create(dir_out, showWarnings = F, recursive = T)
  message(glue(
    "{i}/{nrow(sanctuaries)}: {d_s$sanctuary} ({d_s$nms}) : beg ~ {format(Sys.time(), '%H:%M:%S')}"
  ))

  for (j in 1:length(datasets)) {
    # j = 1

    ds_id <- names(datasets)[j]
    ds    <- datasets[[j]]

    # check what's already processed
    last      <- last_dates |>
      filter(nms == d_s$nms, dataset == ds_id) |>
      pull(last_time)
    avail_end <- as.Date(ds$date_rng[2])

    last_msg <- if (length(last) == 1) as.character(last) else "none"
    message(glue(
      "  {j}/{length(datasets)}: {ds_id}\n",
      "    last retrieved: {last_msg}, available: {ds$date_rng[1]} to {ds$date_rng[2]}"))

    if (length(last) == 1 && last >= avail_end) {
      message(glue("    -> up to date, skipping"))
      next
    }

    # start one month after last processed, or from the beginning
    start_dt <- if (length(last) == 1) {
      format(seq.Date(last, by = "month", length.out = 2)[2],
             "%Y-%m-%dT00:00:00")
    } else {
      ds$date_rng[1]
    }

    message(glue("    -> fetching {start_dt} to {ds$date_rng[2]}"))
    
    # browser()

    nc <- cmt$subset(
      dataset_id           = ds_id,
      variables            = ds$vars,
      minimum_longitude    = bb$xmin,
      maximum_longitude    = bb$xmax,
      minimum_latitude     = bb$ymin,
      maximum_latitude     = bb$ymax,
      start_datetime       = start_dt,
      end_datetime         = ds$date_rng[2],
      minimum_depth        = ds$depth_rng[1],
      maximum_depth        = ds$depth_rng[2],
      output_directory     = dir_out,
      overwrite            = TRUE,
      service              = "timeseries",
      disable_progress_bar = TRUE)
    nc <- as.character(nc)
    message(glue("  : end ~ {format(Sys.time(), '%H:%M:%S')}"))
  }
}
1/15: Cordell Bank (CBNMS) : beg ~ 16:05:31
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
2/15: Channel Islands (CINMS) : beg ~ 16:05:31
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/15: Chumash Proposed Action (CPNMS) : beg ~ 16:05:31
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
4/15: Flower Garden Banks (FGBNMS) : beg ~ 16:05:31
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
5/15: Florida Keys (FKNMS) : beg ~ 16:05:31
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
6/15: Greater Farallones (GFNMS) : beg ~ 16:05:31
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
7/15: Gray's Reef (GRNMS) : beg ~ 16:05:31
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
8/15: Hawaiian Islands Humpback Whale (HIHWNMS) : beg ~ 16:05:31
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
9/15: Monterey Bay (MBNMS) : beg ~ 16:05:31
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
10/15: Monitor (MNMS) : beg ~ 16:05:31
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
11/15: American Samoa (NMSAS) : beg ~ 16:05:31
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
12/15: Olympic Coast (OCNMS) : beg ~ 16:05:32
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
13/15: Stellwagen Bank (SBNMS) : beg ~ 16:05:32
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
14/15: Monterey Bay - Davidson Seamount (MBNMS-david) : beg ~ 16:05:32
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
15/15: Monterey Bay - Mainland (MBNMS-main) : beg ~ 16:05:32
1/3: cmems_mod_glo_phy_my_0.083deg_P1M-m
last retrieved: 2026-02-01, available: 1993-01-01T00:00:00 to 2026-02-01T00:00:00
    -> up to date, skipping
2/3: cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping
3/3: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
last retrieved: 2026-03-01, available: 1997-09-01T00:00:00 to 2026-03-01T00:00:00
    -> up to date, skipping

Show raster

Grabbing first netcdf file, show full metadata on variables and attributes:

Code 
nms = "FKNMS"
ply <- sanctuaries |>
  filter(nms == !!nms)
# nc <- list.files(glue("data/copernicus/{nms}"), full.names = T)[1:6]
nc <- dir_info(glue("data/copernicus/{nms}")) |>
  arrange(desc(modification_time)) |>
  slice(1) |>
  pull(path)

o <- nc_open(nc)
o
File data/copernicus/FKNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M_multi-vars_83.23W-79.98W_24.23N-25.73N_2026-02-01-2026-03-01.nc (NC_FORMAT_NETCDF4):

     3 variables (excluding dimension variables):
        float PP[longitude,latitude,time]   (Contiguous storage)  
            _FillValue: -999
            valid_max: 10000
            units: mg m-2 day-1
            valid_min: 0
            standard_name: primary_productivity_of_biomass_expressed_as_carbon
            long_name: Primary Productivity - Mean of the binned pixels
            coordinates: time lat lon
            missing_value: -999
        short PP_uncertainty[longitude,latitude,time]   (Contiguous storage)  
            _FillValue: -32768
            valid_max: 32767
            units: %
            valid_min: 0
            long_name: Primary Productivity - Uncertainty estimation
            coordinates: time lat lon
            scale_factor: 0.00999999977648258
            missing_value: -32768
        byte flags[longitude,latitude,time]   (Contiguous storage)  
            standard_name: status_flag
            valid_max: 1
            valid_min: 0
            long_name: Flags
            coordinates: time lat lon

     3 dimensions:
        time  Size:2 
            standard_name: time
            unit_long: Days Since 1900-01-01
            axis: T
            long_name: Time
            units: days since 1900-01-01
            calendar: Gregorian
        latitude  Size:37 
            standard_name: latitude
            axis: Y
            units: degrees_north
            long_name: latitude
        longitude  Size:79 
            standard_name: longitude
            axis: X
            units: degrees_east
            long_name: longitude

    9 global attributes:
        comment: average
        title: cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
        source: surface observation
        institution: ACRI
        history: Created using software developed at ACRI-ST
        Conventions: CF-1.8, ACDD-1.3
        contact: servicedesk.cmems@acri-st.fr
        references: http://www.globcolour.info GlobColour has been originally funded by ESA with data from ESA, NASA, NOAA and GeoEye. This version has received funding from the European Community s Seventh Framework Programme ([FP7/2007-2013]) under grant agreement n. 282723 [OSS2015 project].
        copernicusmarine_version: 2.2.5

Using terra::rast() to read the netcdf file, show first layer:

Code 
r <- rast(nc)
Warning: [rast] unknown calendar (assuming standard): Gregorian
Code 
r
class       : SpatRaster 
size        : 37, 79, 6  (nrow, ncol, nlyr)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -83.25, -79.95833, 24.20833, 25.75  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
sources     : cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M_multi-vars_83.23W-79.98W_24.23N-25.73N_2026-02-01-2026-03-01.nc:PP  (2 layers) 
              cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M_multi-vars_83.23W-79.98W_24.23N-25.73N_2026-02-01-2026-03-01.nc:PP_uncertainty  (2 layers) 
              cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M_multi-vars_83.23W-79.98W_24.23N-25.73N_2026-02-01-2026-03-01.nc:flags  (2 layers) 
varnames    : PP (Primary Productivity - Mean of the binned pixels) 
              PP_uncertainty (Primary Productivity - Uncertainty estimation) 
              flags (Flags) 
names       :         PP_1,         PP_2, PP_un~nty_1, PP_un~nty_2, flags_1, flags_2 
unit        : mg m-2 day-1, mg m-2 day-1,           %,           %,        ,         
time (days) : 2026-02-01 to 2026-03-01 (2 steps)
Code 
i <- 1
cat(glue(
  "
  index: {i}
  name: {names(r)[i]}
  time: {time(r[[i]])}
  varname: {varnames(r[[i]])}
  longname: {longnames(r[[i]])}"
))
index: 1
name: PP_1
time: 2026-02-01
varname: PP
longname: Primary Productivity - Mean of the binned pixels
Code 
# glue("{year}/{var}.tif/{time}_{depth}")

Plot the raster with the sanctuary:

Code 
# plot(r[[i]])
r_1 <- r[[i]] |>
  project(leaflet:::epsg3857, method = "bilinear")
# plot(r_1)

plet(
  # r[[i]],
  r_1,
  # main  = glue("{longnames(r[[i]])}\n{names(r)[i]}\n{time(r[[i]])}"),
  main = glue("{longnames(r[[i]])}\nthetao_depth=0.5\n{time(r[[i]])}"),
  tiles = "Esri.OceanBasemap",
  fill_range = T
) |>
  addPolygons(
    data = ply,
    fillOpacity = 0.2
  )

Summarize raster to sanctuary

Code 
extract_tbl <- function(r, p, stat = "mean") {
  terra::extract(r, p, fun = stat, na.rm = T) |>
    pivot_longer(!ID, names_to = "var_it") |>
    mutate(
      var = str_replace(var_it, "_[0-9]+$", ""),
      stat = !!stat,
      time = time(r)
    ) |>
    select(var, time, stat, value)
}

last_dates <- get_last_date_in_csv(results_csv)
ncs <- dir_info(dir_cm, regexp = "\\.nc$", recurse = TRUE) |> pull(path)

for (i in 1:length(ncs)) {
  # i = 1
  nc    <- ncs[i]
  nms   <- str_replace(nc, ".*/copernicus/([A-z-]+)/(.+)_multi-vars.*", "\\1")
  ds_id <- str_replace(nc, ".*/copernicus/([A-z-]+)/(.+)_multi-vars.*", "\\2")

  # parse end date from filename to check if already extracted
  nc_end <- str_extract(basename(nc), "[0-9]{4}-[0-9]{2}-[0-9]{2}\\.nc$") |>
    str_remove("\\.nc$") |> as.Date()
  last <- last_dates |>
    filter(nms == !!nms, dataset == !!ds_id) |>
    pull(last_time)
  if (length(last) == 1 && !is.na(nc_end) && nc_end <= last) {
    message(glue("{i}/{length(ncs)}: {nms}/{ds_id} already extracted through {last}, skipping"))
    next
  }

  message(glue("{i}/{length(ncs)}: nms:{nms}, ds:{ds_id}"))

  r <- rast(nc)
  # plet(r[[1]])
  p <- sanctuaries |>
    filter(nms == !!nms)

  d <- bind_rows(
    extract_tbl(r, p, "mean"),
    extract_tbl(r, p, "sd")
  ) |>
    mutate(
      nms = !!nms,
      dataset = !!ds_id
    ) |>
    relocate(nms, dataset)

  if (file.exists(results_csv)) {
    read_csv(results_csv) |>
      bind_rows(d) |>
      group_by(nms, dataset, var, time, stat) |>
      summarize(
        value = last(value),
        .groups = "drop"
      ) |>
      arrange(nms, dataset, var, time, stat) |>
      write_csv(results_csv)
  } else {
    write_csv(d, results_csv)
  }
}
1/165: CBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
2/165: CBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
3/165: CBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
4/165: CBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
5/165: CBNMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
6/165: CBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
7/165: CBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
8/165: CBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
9/165: CBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
10/165: CBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
11/165: CBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
12/165: CINMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
13/165: CINMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
14/165: CINMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
15/165: CINMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
16/165: CINMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
17/165: CINMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
18/165: CINMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
19/165: CINMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
20/165: CINMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
21/165: CINMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
22/165: CINMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
23/165: CPNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
24/165: CPNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
25/165: CPNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
26/165: CPNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
27/165: CPNMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
28/165: CPNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
29/165: CPNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
30/165: CPNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
31/165: CPNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
32/165: CPNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
33/165: CPNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
34/165: FGBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
35/165: FGBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
36/165: FGBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
37/165: FGBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
38/165: FGBNMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
39/165: FGBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
40/165: FGBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
41/165: FGBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
42/165: FGBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
43/165: FGBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
44/165: FGBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
45/165: FKNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
46/165: FKNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
47/165: FKNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
48/165: FKNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
49/165: FKNMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
50/165: FKNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
51/165: FKNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
52/165: FKNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
53/165: FKNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
54/165: FKNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
55/165: FKNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
56/165: GFNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
57/165: GFNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
58/165: GFNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
59/165: GFNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
60/165: GFNMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
61/165: GFNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
62/165: GFNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
63/165: GFNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
64/165: GFNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
65/165: GFNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
66/165: GFNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
67/165: GRNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
68/165: GRNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
69/165: GRNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
70/165: GRNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
71/165: GRNMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
72/165: GRNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
73/165: GRNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
74/165: GRNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
75/165: GRNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
76/165: GRNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
77/165: GRNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
78/165: HIHWNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
79/165: HIHWNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
80/165: HIHWNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
81/165: HIHWNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
82/165: HIHWNMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
83/165: HIHWNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
84/165: HIHWNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
85/165: HIHWNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
86/165: HIHWNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
87/165: HIHWNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
88/165: HIHWNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
89/165: MBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
90/165: MBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
91/165: MBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
92/165: MBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
93/165: MBNMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
94/165: MBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
95/165: MBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
96/165: MBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
97/165: MBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
98/165: MBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
99/165: MBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
100/165: MBNMS-david/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
101/165: MBNMS-david/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
102/165: MBNMS-david/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
103/165: MBNMS-david/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
104/165: MBNMS-david/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
105/165: MBNMS-david/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
106/165: MBNMS-david/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
107/165: MBNMS-david/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
108/165: MBNMS-david/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
109/165: MBNMS-david/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
110/165: MBNMS-david/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
111/165: MBNMS-main/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
112/165: MBNMS-main/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
113/165: MBNMS-main/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
114/165: MBNMS-main/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
115/165: MBNMS-main/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
116/165: MBNMS-main/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
117/165: MBNMS-main/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
118/165: MBNMS-main/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
119/165: MBNMS-main/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
120/165: MBNMS-main/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
121/165: MBNMS-main/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
122/165: MNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
123/165: MNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
124/165: MNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
125/165: MNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
126/165: MNMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
127/165: MNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
128/165: MNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
129/165: MNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
130/165: MNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
131/165: MNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
132/165: MNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
133/165: NMSAS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
134/165: NMSAS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
135/165: NMSAS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
136/165: NMSAS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
137/165: NMSAS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
138/165: NMSAS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
139/165: NMSAS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
140/165: NMSAS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
141/165: NMSAS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
142/165: NMSAS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
143/165: NMSAS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
144/165: OCNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
145/165: OCNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
146/165: OCNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
147/165: OCNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
148/165: OCNMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
149/165: OCNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
150/165: OCNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
151/165: OCNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
152/165: OCNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
153/165: OCNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
154/165: OCNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
155/165: SBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
156/165: SBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
157/165: SBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
158/165: SBNMS/cmems_mod_glo_phy_my_0.083deg_P1M-m already extracted through 2026-02-01, skipping
159/165: SBNMS/cmems_mod_glo_phy_myint_0.083deg_P1M-m already extracted through 2024-12-01, skipping
160/165: SBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
161/165: SBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
162/165: SBNMS/cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
163/165: SBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
164/165: SBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
165/165: SBNMS/cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M already extracted through 2026-03-01, skipping
Code 
results_url <- results_csv |>
  str_replace(
    here(),
    "https://github.com/noaa-onms/eco-indicators/blob/main"
  )

# TODO: inspect how flags might be distributed across the dataset
# d <- read_csv(results_csv)
# d |>
#   filter(var == "flags") |> 
#   pull(value) |> 
#   table()
#                   0 0.00125156445556946 0.00180831826401446 
#               12204                 678                 678 
# 0.00196850393700787 0.00217391304347826 0.00784313725490196 
#                 678                 678                 678 
#  0.0116279069767442   0.035377456883861  0.0425050794246722 
#                 678                 678                 678 
#  0.0443678254708057  0.0465998661296397  0.0883001161535571 
#                 678                 678                 678 
#   0.107517010378544 
#                 678

# show table
read_csv(results_csv)
# A tibble: 299,760 × 6
   nms   dataset                           var   time                stat  value
   <chr> <chr>                             <chr> <dttm>              <chr> <dbl>
 1 CBNMS cmems_mod_glo_phy_my_0.083deg_P1… bott… 1993-01-01 00:00:00 mean   5.33
 2 CBNMS cmems_mod_glo_phy_my_0.083deg_P1… bott… 1993-01-01 00:00:00 sd     3.81
 3 CBNMS cmems_mod_glo_phy_my_0.083deg_P1… bott… 1993-02-01 00:00:00 mean   5.30
 4 CBNMS cmems_mod_glo_phy_my_0.083deg_P1… bott… 1993-02-01 00:00:00 sd     3.83
 5 CBNMS cmems_mod_glo_phy_my_0.083deg_P1… bott… 1993-03-01 00:00:00 mean   5.19
 6 CBNMS cmems_mod_glo_phy_my_0.083deg_P1… bott… 1993-03-01 00:00:00 sd     3.70
 7 CBNMS cmems_mod_glo_phy_my_0.083deg_P1… bott… 1993-04-01 00:00:00 mean   4.89
 8 CBNMS cmems_mod_glo_phy_my_0.083deg_P1… bott… 1993-04-01 00:00:00 sd     3.38
 9 CBNMS cmems_mod_glo_phy_my_0.083deg_P1… bott… 1993-05-01 00:00:00 mean   4.76
10 CBNMS cmems_mod_glo_phy_my_0.083deg_P1… bott… 1993-05-01 00:00:00 sd     3.15
# ℹ 299,750 more rows

table of results: copernicus.csv

Generate timeseries plot per sanctuary

Code 
d <- read_csv(results_csv)
# table(d$var)
# bottomT  mlotst  so_depth=0.49402499  thetao_depth=0.49402499
#   12288   12288                12288                    12288
#                 bottomT                     CHL         CHL_uncertainty                   DIATO
#                   11820                   10170                   10170                   10170
#       DIATO_uncertainty                    DINO        DINO_uncertainty                   flags
#                   10170                   10170                   10170                   10170
#                   GREEN       GREEN_uncertainty                   HAPTO       HAPTO_uncertainty
#                   10170                   10170                   10170                   10170
#                   MICRO       MICRO_uncertainty                  mlotst                    NANO
#                   10170                   10170                   11820                   10170
#        NANO_uncertainty                    PICO        PICO_uncertainty                 PROCHLO
#                   10170                   10170                   10170                   10170
#     PROCHLO_uncertainty                  PROKAR      PROKAR_uncertainty     so_depth=0.49402499
#                   10170                   10170                   10170                   11820
# thetao_depth=0.49402499
#                   11820

# table(d$dataset, useNA = "ifany")
#               cmems_mod_glo_phy_my_0.083deg_P1M-m
#                                             47280
# cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M
#                                            213570
#       cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M
#                                             30510

plot_var_facets <- function(nms, data) {
  png <- glue("figures/copernicus/{nms}.png")
  data |>
    filter(
      nms == !!nms
    ) |>
    mutate(
      var = str_replace(var, "_depth=0.49402499", "_0.5m")
    ) |>
    pivot_wider(
      names_from = stat,
      values_from = value
    ) |>
    ggplot(aes(x = time, y = mean, color = var)) +
    facet_wrap(~var) +
    geom_line() +
    geom_ribbon(
      aes(ymin = mean - sd, ymax = mean + sd, fill = var),
      color = NA,
      alpha = 0.2
    ) +
    theme(legend.position = "none") +
    facet_wrap(~var, scales = "free_y") +
    labs(
      title = nms,
      x = "Time",
      y = "Value"
    )
  ggsave(png, width = 8, height = 6)
}

sapply(sanctuaries$nms, plot_var_facets, data = d)
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1098 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 2804 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1440 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1062 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1044 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
                               CBNMS                                CINMS 
      "figures/copernicus/CBNMS.png"       "figures/copernicus/CINMS.png" 
                               CPNMS                               FGBNMS 
      "figures/copernicus/CPNMS.png"      "figures/copernicus/FGBNMS.png" 
                               FKNMS                                GFNMS 
      "figures/copernicus/FKNMS.png"       "figures/copernicus/GFNMS.png" 
                               GRNMS                              HIHWNMS 
      "figures/copernicus/GRNMS.png"     "figures/copernicus/HIHWNMS.png" 
                               MBNMS                                 MNMS 
      "figures/copernicus/MBNMS.png"        "figures/copernicus/MNMS.png" 
                               NMSAS                                OCNMS 
      "figures/copernicus/NMSAS.png"       "figures/copernicus/OCNMS.png" 
                               SBNMS                          MBNMS-david 
      "figures/copernicus/SBNMS.png" "figures/copernicus/MBNMS-david.png" 
                          MBNMS-main 
 "figures/copernicus/MBNMS-main.png"

Show timeseries per sanctuary

Show timeseries per sanctuary, without uncertainty

Code 
plot_var_facets_without_uncertainty <- function(nms, data) {
  png <- glue("figures/copernicus/{nms}_without_uncertainty.png")
  data |>
    filter(
      nms == !!nms,
      !str_ends(var, "_uncertainty")
    ) |>
    mutate(
      var = str_replace(var, "_depth=0.49402499", "_0.5m")
    ) |>
    pivot_wider(
      names_from  = stat,
      values_from = value
    ) |>
    ggplot(aes(x = time, y = mean, color = var)) +
    facet_wrap(~var) +
    geom_line() +
    geom_ribbon(
      aes(ymin = mean - sd, ymax = mean + sd, fill = var),
      color = NA,
      alpha = 0.2
    ) +
    theme(legend.position = "none") +
    facet_wrap(~var, scales = "free_y") +
    labs(
      title = nms,
      x     = "Time",
      y     = "Value"
    )
  ggsave(png, width = 8, height = 6)
}

sapply(sanctuaries$nms, plot_var_facets_without_uncertainty, data = d)
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 549 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 2282 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 720 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 531 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 522 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
                                                   CBNMS 
      "figures/copernicus/CBNMS_without_uncertainty.png" 
                                                   CINMS 
      "figures/copernicus/CINMS_without_uncertainty.png" 
                                                   CPNMS 
      "figures/copernicus/CPNMS_without_uncertainty.png" 
                                                  FGBNMS 
     "figures/copernicus/FGBNMS_without_uncertainty.png" 
                                                   FKNMS 
      "figures/copernicus/FKNMS_without_uncertainty.png" 
                                                   GFNMS 
      "figures/copernicus/GFNMS_without_uncertainty.png" 
                                                   GRNMS 
      "figures/copernicus/GRNMS_without_uncertainty.png" 
                                                 HIHWNMS 
    "figures/copernicus/HIHWNMS_without_uncertainty.png" 
                                                   MBNMS 
      "figures/copernicus/MBNMS_without_uncertainty.png" 
                                                    MNMS 
       "figures/copernicus/MNMS_without_uncertainty.png" 
                                                   NMSAS 
      "figures/copernicus/NMSAS_without_uncertainty.png" 
                                                   OCNMS 
      "figures/copernicus/OCNMS_without_uncertainty.png" 
                                                   SBNMS 
      "figures/copernicus/SBNMS_without_uncertainty.png" 
                                             MBNMS-david 
"figures/copernicus/MBNMS-david_without_uncertainty.png" 
                                              MBNMS-main 
 "figures/copernicus/MBNMS-main_without_uncertainty.png"

Notes on PFT Uncertainty Variables

The *_uncertainty variables in the GlobColour plankton dataset exhibit pronounced inflections toward high values across all sanctuaries. These patterns are methodological, not ecological, and reflect the uncertainty estimation framework described below.

What the uncertainty values represent. Uncertainty is expressed as percent error bars (relative uncertainty) on each PFT chlorophyll-a concentration estimate. Values are stored as 2-byte signed integers with a scale factor of 0.01, producing a maximum representable value of 327.67 % (Copernicus Marine Environment Monitoring Service 2025b).

How uncertainty is computed. The per-pixel uncertainty is derived via error propagation combining three independent sources (Xi et al. 2021, Eq. 7):

  1. Remote-sensing reflectance (\(R_{rs}\)) input uncertainty — estimated via Monte Carlo perturbation of the \(R_{rs}\) spectrum, propagated through the EOF-regression retrieval algorithm using a look-up table (LUT) approach;
  2. Model parameter uncertainty — arising from the regression coefficients relating EOF scores to PFT concentrations;
  3. Sea surface temperature (SST) uncertainty — SST is an auxiliary predictor in the EOF-SST hybrid retrieval, and its uncertainty is propagated into the final PFT estimate.

Algorithm origin and sensor dependence. The EOF-SST hybrid algorithm was originally developed for OLCI (Sentinel-3), which has 21 spectral bands, and subsequently extended to the multi-sensor merged GlobColour product (Xi et al. 2020). Because the merged product combines sensors with different spectral configurations, uncertainty is generally higher for sensors with fewer bands (e.g., MODIS Aqua has 9 ocean-colour bands vs. OLCI’s 21).

Sensor transitions drive the inflections. The multi-sensor record combines: SeaWiFS (1997–2010), MODIS Aqua (2002–present), MERIS (2002–2012), VIIRS NPP (2012–present), and OLCI-S3A (2016–present) (Garnesson et al. 2019; Maritorena et al. 2010). Each sensor transition changes the input \(R_{rs}\) characteristics and the number of available spectral bands, producing step-changes in the propagated uncertainty even when the underlying ocean state is stable.

Low concentrations amplify relative uncertainty. PFTs that occur at very low concentrations in temperate and subpolar waters — particularly prokaryotes (PROKAR) and Prochlorococcus (PROCHLO) — yield high percent uncertainties because a small absolute error translates to a large relative error.

Reprocessing discontinuities. The November 2023 reprocessing (NASA R2022, MERIS 4th reprocessing, OLCI baseline collection 3.01) may introduce discontinuities in both concentration and uncertainty time series (Copernicus Marine Environment Monitoring Service 2025b).

Known sensor drifts. NASA sensors exhibit documented radiometric drifts that affect the merged product. The SQO report flags these as a known limitation (Copernicus Marine Environment Monitoring Service 2025a).

Weakest retrieval skill. Prokaryotes and Prochlorococcus have the lowest retrieval skill among all PFTs (\(R^2\) = 0.42–0.51 against in-situ validation) (Copernicus Marine Environment Monitoring Service 2025a; Xi et al. 2021), meaning their uncertainty estimates are both the largest and least well-constrained.

VIIRS-SNPP post-2017 exclusions. Data from VIIRS-SNPP after 2017 have been excluded from some products due to identified spurious trends in the merged record (Xi et al. 2025).

Using PFT Data for Multivariate Ordination

Despite the uncertainty caveats above, the PFT concentration data are well-validated and suitable for multivariate analyses. Below are recommended approaches for detecting ecosystem-level shifts across sanctuaries.

  1. Use concentration variables only (exclude _uncertainty). The PFT chlorophyll-a concentrations themselves are the validated retrieval products; the _uncertainty fields are informational metadata describing retrieval confidence. They should not be included as ordination input variables alongside the concentrations.

  2. Log-transform concentrations. PFT concentrations span orders of magnitude. Log-transformation stabilizes variance and matches the scale on which the retrieval algorithm was trained and validated (Xi et al. 2021). Use \(\log_{10}(\text{Chl-}a + \epsilon)\) where \(\epsilon\) is a small constant to handle zeros.

  3. Exclude low-skill PFTs or aggregate them. Consider dropping PROKAR and PROCHLO from ordination (weakest \(R^2 < 0.55\)), or combine them into a single “prokaryote” group to improve signal-to-noise ratio.

  4. PCA / RDA on a sanctuary \(\times\) time matrix. Structure the data as: rows = sanctuary-month combinations, columns = PFT variables (CHL, DIATO, DINO, GREEN, HAPTO, MICRO, NANO, PICO) plus physical variables (thetao, so, mlotst, bottomT, PP). Redundancy analysis (RDA) can then partition biological community variation explained by environmental drivers.

  5. Use uncertainty as weights. In weighted ordination methods (e.g., weighted PCA), use \(1 / \text{uncertainty}\) as observation weights to downweight unreliable observations rather than discarding them entirely.

  6. Temporal anomaly approach. Compute anomalies relative to a monthly climatology per sanctuary to remove the dominant seasonal signal, then ordinate on the anomaly matrix to detect regime shifts and long-term trends.

  7. Non-metric MDS (NMDS). NMDS is robust to the non-normal distributions common in ocean-colour data. Use Bray-Curtis dissimilarity on the PFT concentration matrix (after log-transformation) to preserve compositional distances.

  8. Rolling window analysis. Apply ordination in sliding temporal windows (e.g., 5-year windows advancing by 1 year) to track community composition shifts over time, overlaid with environmental drivers (SST, salinity, mixed layer depth).

References

Copernicus Marine Environment Monitoring Service. 2025a. “Ocean Colour Thematic Assembly Centre Scientific Quality Overview.” CMEMS-OC-SQO, Issue 6.0. CMEMS.
———. 2025b. “Ocean Colour Thematic Assembly Centre Product User Manual.” CMEMS-OC-PUM, Issue 6.0. CMEMS.
Garnesson, Philippe, Antoine Mangin, Odile Fanton d’Andon, Julien Demaria, and Marine Bretagnon. 2019. “The CMEMS GlobColour Chlorophyll a Product Based on Satellite Observation: Multi-Sensor Merging and Flagging Strategies.” Ocean Science 15 (3): 819–30. https://doi.org/10.5194/os-15-819-2019.
Maritorena, Stéphane, Odile Hembise Fanton d’Andon, Antoine Mangin, and David A. Siegel. 2010. “Merged Satellite Ocean Color Data Products Using a Bio-Optical Model: Characteristics, Benefits and Issues.” Remote Sensing of Environment 114 (8): 1791–1804. https://doi.org/10.1016/j.rse.2010.04.002.
Xi, Hongyan, Marine Bretagnon, Ehsan Mehdipour, Julien Demaria, Antoine Mangin, and Astrid Bracher. 2025. “Consistent Long-Term Observations of Surface Phytoplankton Functional Types from Space.” State of the Planet 6-osr9: 7. https://doi.org/10.5194/sp-6-osr9-7-2025.
Xi, Hongyan, Svetlana N. Losa, Antoine Mangin, Philippe Garnesson, Marine Bretagnon, Julien Demaria, Mariana A. Soppa, Odile Hembise Fanton d’Andon, and Astrid Bracher. 2021. “Global Chlorophyll a Concentrations of Phytoplankton Functional Types with Detailed Uncertainty Assessment Using Multisensor Ocean Color and Sea Surface Temperature Satellite Products.” Journal of Geophysical Research: Oceans 126 (5): e2020JC017127. https://doi.org/10.1029/2020JC017127.
Xi, Hongyan, Svetlana N. Losa, Antoine Mangin, Mariana A. Soppa, Philippe Garnesson, Julien Demaria, Yangyang Liu, Odile Hembise Fanton d’Andon, and Astrid Bracher. 2020. “Global Retrieval of Phytoplankton Functional Types Based on Empirical Orthogonal Functions Using CMEMS GlobColour Merged Products and Further Extension to OLCI Data.” Remote Sensing of Environment 240: 111704. https://doi.org/10.1016/j.rse.2020.111704.