Adds plot_data_availability()

NAMESPACE

@@ -14,6 +14,7 @@ export(plot_abundance_at_age)
 export(plot_biomass)
 export(plot_biomass_at_age)
 export(plot_catch_at_age)
+export(plot_data_availability)
 export(plot_discard)
 export(plot_fishing_mortality)
 export(plot_index)

R/plot_data_availability.R

@@ -0,0 +1,173 @@
+#' Plot the data available to an assessment model
+#'
+#' Plot the years in which each fleet has each type of data (landings,
+#' discards, indices of abundance, and age and length compositions), as in
+#' `r4ss` (`SSplotData()`), `stockassessment` (`dataplot()`), and `spict`
+#' (`plotspict.data()`). Each type of data is a panel, each fleet a row, and
+#' each year with data a square.
+#'
+#' @inheritParams plot_spawning_biomass
+#' @param fleet A string or vector of strings. Fleets to plot.
+#'
+#' Default: NULL (all fleets)
+#'
+#' @returns A `ggplot2` object.
+#'
+#' @details The data are the observations in output from [convert_output()],
+#' found by their labels:
+#'
+#' * Landings: "landings_observed", "landings_observed_weight",
+#'   "landings_observed_numbers", or "catch_observed", or "catch" when the
+#'   output also has "catch_predicted". Only the first of these labels in
+#'   the output is used, so that landings reported in more than one table
+#'   are not shown twice.
+#' * Discards: "discard_observed".
+#' * Index: "index_observed".
+#' * Age and length compositions: "age_composition_observed" and
+#'   "length_composition_observed", or the composition fit tables of Stock
+#'   Synthesis ("all_obs_mean" in the age and length composition modules).
+#' * Compositions: "composition_observed", when the type of composition is
+#'   not reported.
+#'
+#' A year has data when it has an observation greater than zero. All
+#' modules are plotted, so `interactive` is not used; it is accepted for
+#' consistency with the other plotting functions (e.g., in
+#' [save_all_plots()]).
+#'
+#' @export
+#'
+#' @examples
+#' plot_data_availability(stockplotr::example_data)
+plot_data_availability <- function(
+  dat,
+  fleet = NULL,
+  interactive = TRUE,
+  make_rda = FALSE,
+  figures_dir = getwd()
+) {
+  check_data(dat)
+  avail <- data_availability(dat, fleet = fleet)
+
+  plt <- ggplot2::ggplot(avail, ggplot2::aes(x = .data$year, y = .data$fleet)) +
+    ggplot2::geom_point(ggplot2::aes(color = .data$fleet), shape = 15, size = 2.5) +
+    ggplot2::facet_grid(
+      rows = ggplot2::vars(.data$type),
+      scales = "free_y",
+      space = "free_y",
+      switch = "y"
+    ) +
+    ggplot2::scale_y_discrete(limits = rev) +
+    ggplot2::labs(x = "Year", y = NULL) +
+    theme_noaa() +
+    ggplot2::theme(
+      legend.position = "none",
+      strip.placement = "outside",
+      strip.text.y.left = ggplot2::element_text(angle = 0, hjust = 1),
+      panel.grid.major.y = ggplot2::element_line(color = "gray90")
+    )
+
+  if (make_rda) {
+    # Obtain relevant key quantities for captions/alt text
+    data.availability.start.year <- min(avail$year)
+    data.availability.end.year <- max(avail$year)
+    data.availability.n.fleets <- length(unique(avail$fleet))
+    data.availability.types <- data_availability_list(levels(droplevels(avail$type)))
+
+    # calculate & export key quantities
+    export_kqs(
+      data.availability.start.year,
+      data.availability.end.year,
+      data.availability.n.fleets,
+      data.availability.types,
+      dir = figures_dir
+    )
+
+    # Add key quantities to captions/alt text
+    insert_kqs(
+      data.availability.start.year,
+      data.availability.end.year,
+      data.availability.n.fleets,
+      data.availability.types,
+      dir = figures_dir
+    )
+
+    create_rda(
+      object = plt,
+      topic_label = "data_availability",
+      fig_or_table = "figure",
+      dat = dat,
+      dir = figures_dir
+    )
+  }
+  plt
+}
+
+# The years in which each fleet has each type of data: one row per type,
+# fleet, and year, with `type` a factor in the order the types are plotted
+data_availability <- function(dat, fleet = NULL) {
+  if ("model" %in% names(dat) && length(unique(stats::na.omit(dat$model))) > 1) {
+    cli::cli_abort("{.arg dat} has output from more than one model; plot one model at a time.")
+  }
+  label <- dat$label
+  module <- tolower(dat$module_name)
+  # Landings: the first label present, so landings in more than one table
+  # are not shown twice
+  landings_labels <- c(
+    "landings_observed", "landings_observed_weight", "landings_observed_numbers",
+    "catch_observed"
+  )
+  if ("catch_predicted" %in% label) {
+    landings_labels <- c(landings_labels, "catch")
+  }
+  landings_label <- landings_labels[landings_labels %in% label][1]
+  age_fit <- label %in% "all_obs_mean" & module %in% c("fit_age_comps", "age_comp_fit_table")
+  length_fit <- label %in% "all_obs_mean" & module %in% c("fit_len_comps", "len_comp_fit_table")
+  type <- dplyr::case_when(
+    label %in% landings_label ~ "Landings",
+    label %in% "discard_observed" ~ "Discards",
+    label %in% "index_observed" ~ "Index",
+    label %in% "age_composition_observed" | age_fit ~ "Age compositions",
+    label %in% "length_composition_observed" | length_fit ~ "Length compositions",
+    label %in% "composition_observed" ~ "Compositions",
+    TRUE ~ NA_character_
+  )
+  estimate <- suppressWarnings(as.numeric(dat$estimate))
+  keep <- !is.na(type) & !is.na(dat$year) & !is.na(estimate) & estimate > 0
+  if (!any(keep)) {
+    cli::cli_abort("No observed data found in {.arg dat}; see {.fn plot_data_availability} for the labels used.")
+  }
+  avail <- tibble::tibble(
+    type = type[keep],
+    fleet = as.character(dat$fleet[keep]),
+    year = as.numeric(dat$year[keep])
+  )
+  if (any(is.na(avail$fleet))) {
+    cli::cli_alert_info(
+      "{sum(is.na(avail$fleet))} observation{?s} without a fleet {?is/are} shown as fleet {.val unknown}."
+    )
+    avail$fleet[is.na(avail$fleet)] <- "unknown"
+  }
+  if (!is.null(fleet)) {
+    missing <- setdiff(fleet, avail$fleet)
+    if (length(missing) > 0) {
+      cli::cli_alert_info("No observed data for fleet{?s} {.val {missing}}.")
+    }
+    avail <- avail[avail$fleet %in% fleet, , drop = FALSE]
+    if (nrow(avail) == 0) {
+      cli::cli_abort("No observed data for the selected fleets.")
+    }
+  }
+  types <- c("Landings", "Discards", "Index", "Age compositions", "Length compositions", "Compositions")
+  avail <- dplyr::distinct(avail)
+  avail$type <- factor(avail$type, levels = types[types %in% avail$type])
+  dplyr::arrange(avail, .data$type, .data$fleet, .data$year)
+}
+
+# "a", "a and b", or "a, b, and c" (lower case), for captions
+data_availability_list <- function(x) {
+  x <- tolower(x)
+  if (length(x) <= 2) {
+    return(paste(x, collapse = " and "))
+  }
+  paste0(paste(x[-length(x)], collapse = ", "), ", and ", x[length(x)])
+}

inst/resources/captions_alt_text_template.csv

@@ -48,3 +48,4 @@ catchability,table,Catchability over time for fleet catchability.fleet.,
 sensitivity.runs,table,"Forecasted catch, biomass, spawning biomass, recruitment, fishing mortality, and stock status as projected from the base model configuration.",
 catch_at_age,figure,Fishery age composition (caa.start.year-caa.end.year). The area of the circle is proportional to the catch. Diagonal lines indicated the top 5\\% strongest year classes.,"Bubble plot showing estimated total catch at age. The x axis shows the year, which spans from caa.start.year to caa.end.year. The y axis shows age, which spans from caa.age.min to caa.age.max. The size of the bubbles range from tot.catch.min to tot.catch.max."
 age_to_length_conversion,figure,"Age-to-length conversion shown as the probability distribution of length bins at each age. Density curves show the proportion of fish in each length bin at a given age, with wider curves indicating greater variability in length at age.","A set of density plots showing the age-to-length conversion for species, with one panel for each morph. The x axis shows age in years and the y axis shows length bins (in length.units). At each age, a vertical density curve shows the distribution of lengths. Mean length increases from about length.at.min.age length.units at age min.age to about length.at.max.age length.units at age max.age, after which the distributions stabilize, indicating that growth slows as fish age. "
+data_availability,figure,"Years with data for each fleet, by type of data (data.availability.types).","Grid of squares showing the years in which each fleet has each type of data. The x axis shows the year, which spans from data.availability.start.year to data.availability.end.year. The panels show the types of data (data.availability.types), and each row within a panel is one of data.availability.n.fleets fleets; a square marks a year with data."

inst/resources/key_quantity_template.csv

@@ -29,6 +29,10 @@ cpue.max,,maximum catch per unit effort,,,,
 cpue.min,,minimum catch per unit effort,,,,
 cpue.start.year,,start year of the CPUE time series,,,,
 cpue.units,,units for catch per unit effort,,,,
+data.availability.end.year,,final year with data in the data availability plot,,,,
+data.availability.n.fleets,,number of fleets with data in the data availability plot,,,,
+data.availability.start.year,,first year with data in the data availability plot,,,,
+data.availability.types,,types of data in the data availability plot,,,,
 F.FMSY.end.yr,,fishing mortality relative to FMSY in the final year,,,,
 F.MSY.terminal,,fishing mortality relative to MSY in the terminal year of fishing mortality,,,,
 F.MSY.terminal.max,,upper estimate of fishing mortality relative to MSY of the 95th percentile confidence interval for the terminal year of fishing mortality,,,,

tests/testthat/test-plot_data_availability.R

@@ -0,0 +1,78 @@
+test_that("plot_data_availability works with each example model output", {
+  est <- readRDS(test_path("fixtures", "fims_output", "fims_estimates.rds"))
+  datasets <- list(
+    example_data = stockplotr::example_data,
+    pollock = stockplotr::pollock,
+    ss3_simple = stockplotr::ss3_simple,
+    vermilion_snapper = stockplotr::vermilion_snapper,
+    fims = suppressMessages(convert_output(est, model = "fims", start_year = 2001))
+  )
+  for (name in names(datasets)) {
+    p <- plot_data_availability(datasets[[name]])
+    expect_s3_class(p, "ggplot")
+    expect_no_error(ggplot2::ggplot_build(p))
+    expect_gt(nrow(p$data), 0)
+  }
+})
+
+test_that("data_availability finds each type of data once per fleet and year", {
+  avail <- data_availability(stockplotr::example_data)
+  expect_equal(
+    levels(avail$type),
+    c("Landings", "Discards", "Index", "Age compositions", "Length compositions")
+  )
+  expect_equal(nrow(avail), nrow(dplyr::distinct(avail)))
+  # Index years match the observed index
+  obs <- stockplotr::example_data[stockplotr::example_data$label %in% "index_observed", ]
+  index <- avail[avail$type == "Index", ]
+  expect_setequal(paste(index$fleet, index$year), unique(paste(obs$fleet, obs$year)))
+  # Landings come from one table only (not also from catch_observed, whose
+  # fleets are numbered)
+  expect_setequal(unique(avail$fleet[avail$type == "Landings"]), c("North", "South"))
+
+  # pollock's catch is observed (it has catch_predicted); its compositions
+  # are of unknown type
+  avail <- data_availability(stockplotr::pollock)
+  expect_true(all(c("Landings", "Compositions") %in% avail$type))
+})
+
+test_that("data_availability selects fleets and reports problems", {
+  avail <- data_availability(stockplotr::vermilion_snapper, fleet = c("rhb", "chl"))
+  expect_setequal(unique(avail$fleet), c("rhb", "chl"))
+  expect_message(
+    data_availability(stockplotr::vermilion_snapper, fleet = c("rhb", "nope")),
+    "No observed data for fleet"
+  )
+  expect_error(
+    suppressMessages(data_availability(stockplotr::vermilion_snapper, fleet = "nope")),
+    "No observed data for the selected fleets"
+  )
+  dat <- stockplotr::vermilion_snapper
+  expect_error(
+    data_availability(dat[!grepl("observed", dat$label), ]),
+    "No observed data found"
+  )
+  # Observations without a fleet are kept, not dropped
+  row <- which(dat$label == "index_observed" & dat$estimate > 0 & !is.na(dat$year))[1]
+  dat$fleet[row] <- NA
+  expect_message(avail <- data_availability(dat), "without a fleet")
+  expect_true("unknown" %in% avail$fleet)
+})
+
+test_that("plot_data_availability saves a caption, alt text, and key quantities", {
+  dir <- tempfile("data_availability")
+  dir.create(dir)
+  on.exit(unlink(dir, recursive = TRUE), add = TRUE)
+  suppressMessages(invisible(utils::capture.output(
+    plot_data_availability(stockplotr::example_data, make_rda = TRUE, figures_dir = dir)
+  )))
+  rda <- NULL
+  load(file.path(dir, "figures", "data_availability_figure.rda"))
+  expect_no_match(c(rda$caption, rda$alt_text), "data.availability.", fixed = TRUE)
+  expect_match(rda$caption, "landings, discards, index, age compositions, and length compositions", fixed = TRUE)
+  expect_match(rda$alt_text, "from 1876 to 2022", fixed = TRUE)
+  kqs <- utils::read.csv(file.path(dir, "key_quantities.csv"))
+  kqs <- stats::setNames(kqs$value, kqs$key_quantity)
+  expect_equal(as.character(kqs[["data.availability.n.fleets"]]), "4")
+  expect_equal(as.character(kqs[["data.availability.start.year"]]), "1876")
+})

tests/testthat/test-save_all_plots.R

@@ -8,6 +8,7 @@ fig_base_temp_files <- c(
   "natural_mortality_figure.rda",
   "biomass_at_age_figure.rda",
   "catch_at_age_figure.rda",
+  "data_availability_figure.rda",
   "recruitment_deviations_figure.rda",
   "recruitment_figure.rda",
   "selectivity_figure.rda",