R/ifcb_extract_biovolumes.R

Defines functions ifcb_extract_biovolumes

Documented in ifcb_extract_biovolumes

utils::globalVariables(c("biovolume", "roi", "roi_number", "Biovolume", "cell_count", "cell_count_resolved"))
#' Extract Biovolumes from IFCB Data and Compute Carbon Content
#'
#' This function reads biovolume data from feature files generated by the `ifcb-analysis` repository (Sosik and Olson 2007)
#' and matches them with corresponding classification results or manual annotations. It calculates biovolume in cubic micrometers and
#' determines if each class is a diatom based on the World Register of Marine Species (WoRMS). Carbon content
#' is computed for each region of interest (ROI) using conversion functions from Menden-Deuer and Lessard (2000),
#' depending on whether the class is identified as a diatom.
#'
#' @param feature_files A path to a folder containing feature files or a character vector of file paths.
#' @param class_files (Optional) A character vector of full paths to classification or manual
#'   annotation files (`.mat`, `.h5`, or `.csv`), or a single path to a folder
#'   containing such files. Supply a single file format per sample: a sample
#'   represented twice (e.g. by both a `.mat` and a `.h5`) would have its ROIs
#'   counted once per file, so this is rejected with an error naming the
#'   affected samples.
#' @param custom_images (Optional) A character vector of image filenames in the format DYYYYMMDDTHHMMSS_IFCBXXX_ZZZZZ(.png),
#'        where "XXX" represents the IFCB number and "ZZZZZ" represents the ROI number.
#'        These filenames should match the `roi_number` assignment in the `feature_files` and can be
#'        used as a substitute for classification files.
#' @param custom_classes (Optional) A character vector of corresponding class labels for `custom_images`.
#' @param class2use_file (Optional) A character string specifying the path to the file containing the `class2use` variable. Only required for manual results (default: NULL).
#' @param micron_factor Conversion factor from microns per pixel (default: 1/3.4).
#' @param diatom_class A character vector specifying diatom class names in WoRMS. Default: `"Bacillariophyceae"`.
#' @param diatom_include Optional character vector of class names that should always be treated as diatoms,
#'        overriding the boolean result of \code{ifcb_is_diatom}. Default: NULL.
#' @param marine_only Logical. If `TRUE`, restricts the WoRMS search to marine taxa only. Default: `FALSE`.
#' @param diatom_equation A character string selecting which Menden-Deuer and Lessard (2000)
#'   carbon-to-volume relationship to apply to diatoms. `"large"` (default) uses the
#'   large-diatom (> 3000 micron^3) equation (`vol2C_lgdiatom`), matching the
#'   `ifcb-analysis` convention. `"all"` uses the all-sizes diatom equation
#'   (`vol2C_diatom`), which assigns more carbon to small cells. `"auto"`
#'   (`vol2C_diatom_auto`) selects between them per volume, using the
#'   large-diatom equation above 3000 micron^3 and the all-sizes equation at or
#'   below it, so each stays inside its calibrated size range. `"auto"` needs no
#'   `cell_count` data and is independent of `carbon_conversion`. Be aware that the
#'   two are discontinuous at that boundary (about 190 pgC against 135 pgC), so
#'   `"auto"` makes predicted carbon fall as a cell grows across it; that is why
#'   it is not the default. `"auto"` selects on whatever volume the conversion
#'   receives, which is the per-cell volume when `carbon_conversion = "cell"`.
#'   Note that biovolume is measured per region of interest (ROI/image), so it is
#'   not a per-cell volume unless `carbon_conversion = "cell"` is used: chains of
#'   small cells register a large ROI biovolume.
#'   Non-diatom protists always use `vol2C_nondiatom` regardless of this setting.
#' @param threshold A character string controlling which classification to use.
#'   `"opt"` (default) uses the threshold-applied classification, where
#'   predictions below the per-class optimal threshold are labeled
#'   `"unclassified"`. Any other value (e.g. `"all"`) uses the raw winning
#'   class without any threshold applied.
#' @param multiblob Logical. If `TRUE`, includes multiblob features. Default: `FALSE`.
#' @param feature_recursive Logical. If `TRUE`, searches recursively for feature files when `feature_files` is a folder. Default: `TRUE`.
#' @param class_recursive Logical. If `TRUE` and `class_files` is a folder, searches recursively for classification files. Default: `TRUE`.
#' @param drop_zero_volume Logical. If `TRUE`, rows where `Biovolume` equals zero (e.g., artifacts such as smudges on the flow cell) are removed. Default: `FALSE`.
#' @param feature_version Optional numeric or character version to filter feature files by (e.g. 2 for "_v2"). Default is NULL (no filtering).
#' @param use_cell_counts Logical. If `TRUE`, reads the optional per-ROI `cell_count`
#'        data stored by the diatom chain counter in `.mat`/`.h5`/`.csv` classification files and
#'        adds `cell_count` (raw) and `cell_count_resolved` (resolved abundance) columns to the
#'        output. Only supported with automated `class_files`; not with manual files
#'        or `custom_images`. Default: `FALSE`.
#' @param single_cell_values Integer vector of `cell_count` values that should be treated
#'        as a single cell when resolving `cell_count_resolved`. Default is `c(-1, 0)`, i.e. ROIs that
#'        were not counted (`-1`) and ROIs where no cells were detected (`0`) each count as one
#'        cell. Values not listed are used verbatim. Only used when `use_cell_counts = TRUE`.
#' @param carbon_conversion A character string controlling how the Menden-Deuer and
#'        Lessard (2000) relationships are applied. `"roi"` (default) applies the
#'        selected equation once to the whole ROI biovolume, matching the
#'        `ifcb-analysis` convention and reproducing previous results exactly.
#'        `"cell"` divides the ROI biovolume by the resolved cell count, applies
#'        the equation to that per-cell volume, and multiplies back by the count,
#'        which is how the relationships are defined (`log pgC cell^-1`). Requires
#'        `use_cell_counts = TRUE`. `carbon_pg` remains a per-ROI total in both
#'        cases, not carbon per cell. ROIs holding a single cell, ROIs that were
#'        never chain-counted, and ROIs from files carrying no `cell_count` data
#'        are all converted as one cell, so only ROIs with `cell_count >= 2`
#'        change.
#' @param use_python Logical. If `TRUE`, attempts to read `.mat` files using a Python-based method (`SciPy`). Default: `FALSE`.
#' @param verbose Logical. If `TRUE`, prints progress messages. Default: `TRUE`.
#' @param mat_folder `r lifecycle::badge("deprecated")`
#'    Use \code{class_files} instead.
#' @param mat_files `r lifecycle::badge("deprecated")`
#'    Use \code{class_files} instead.
#' @param mat_recursive `r lifecycle::badge("deprecated")`
#'    Use \code{class_recursive} instead.
#'
#' @return A data frame containing:
#' - `sample`: The sample name.
#' - `classifier`: The classifier used (if applicable).
#' - `roi_number`: The region of interest (ROI) number.
#' - `class`: The identified taxonomic class.
#' - `biovolume_um3`: Computed biovolume in cubic micrometers.
#' - `carbon_pg`: Estimated carbon content in picograms.
#' - `cell_count`, `cell_count_resolved` (only when `use_cell_counts = TRUE`): the raw per-ROI
#'   cell count and the resolved number of cells used for abundance.
#'
#' @details
#' - **Classification Data Handling:**
#'   - If `class_files` is provided, the function reads class annotations from `.mat`, `.h5`, or `.csv` files.
#'   - If `custom_images` and `custom_classes` are supplied, they override classification file data (e.g. data from a CNN model).
#'   - If both `class_files` and `custom_images/custom_classes` are given, `class_files` takes precedence.
#'
#' - **MAT File Processing:**
#'   - If `use_python = TRUE`, the function reads `.mat` files using `ifcb_read_mat()` (requires Python + `SciPy`).
#'   - Otherwise, it reads `.mat` files with the default R reader.
#'
#' - **Per-cell carbon conversion:**
#'   - The Menden-Deuer and Lessard (2000) relationships are fitted per cell
#'     (`log pgC cell^-1 = log a + b * log V`), but an IFCB biovolume describes a
#'     whole region of interest, which for a chain-forming diatom is the whole
#'     chain. Every one of these relationships has `b < 1`, so applying one to an
#'     aggregated chain volume returns less carbon than applying it per cell and
#'     summing. The two differ by a factor of `n^(1-b)`: about 1.28 for an
#'     8-cell chain and 1.43 for 20 cells under the large-diatom equation.
#'   - `carbon_conversion = "cell"` divides the ROI biovolume evenly among the
#'     counted cells, which assumes the cells in a chain are of similar size.
#'   - It further assumes the ROI biovolume is cell volume. This is weakest for
#'     *Chaetoceros*, whose setae add to the measured ROI biovolume without being
#'     cell material, so the per-cell volume is overestimated and per-cell carbon
#'     is biased high; whole-ROI conversion biases it low instead. Neither is
#'     corrected here.
#'
#' @examples
#' \dontrun{
#' # Using classification results:
#' feature_files <- "data/features"
#' class_files <- "data/classified"
#'
#' biovolume_df <- ifcb_extract_biovolumes(feature_files,
#'                                         class_files)
#'
#' print(biovolume_df)
#'
#' # Using custom classification result:
#' classes <- c("Mesodinium_rubrum",
#'              "Mesodinium_rubrum")
#' images <- c("D20220522T003051_IFCB134_00002",
#'            "D20220522T003051_IFCB134_00003")
#'
#' biovolume_df_custom <- ifcb_extract_biovolumes(feature_files,
#'                                                custom_images = images,
#'                                                custom_classes = classes)
#'
#' print(biovolume_df_custom)
#' }
#'
#' @references Menden-Deuer Susanne, Lessard Evelyn J., (2000), Carbon to volume relationships for dinoflagellates, diatoms, and other protist plankton, Limnology and Oceanography, 45(3), 569-579, doi: 10.4319/lo.2000.45.3.0569.
#' @references Sosik, H. M. and Olson, R. J. (2007), Automated taxonomic classification of phytoplankton sampled with imaging-in-flow cytometry. Limnol. Oceanogr: Methods 5, 204–216.
#' @references Groves, G. J. J., Arthur, G., Bresnan, E., Whyte, C., Arce, P. and Davidson, K. (2026), Automatic enumeration of chains of marine diatoms using "You Only Look Once" - a machine learning approach. Journal of Plankton Research, 48(2), fbaf064, doi: 10.1093/plankt/fbaf064.
#'
#' @export
#'
#' @seealso \code{\link{ifcb_read_features}} \code{\link{ifcb_is_diatom}} \code{\link{ifcb_summarize_cell_counts}} \url{https://github.com/nodc-sweden/ifcb-pytorch-classify} \url{https://www.marinespecies.org/}
ifcb_extract_biovolumes <- function(feature_files, class_files = NULL, custom_images = NULL, custom_classes = NULL,
                                    class2use_file = NULL, micron_factor = 1 / 3.4,
                                    diatom_class = "Bacillariophyceae", diatom_include = NULL, marine_only = FALSE,
                                    diatom_equation = c("large", "all", "auto"),
                                    threshold = "opt", multiblob = FALSE, feature_recursive = TRUE,
                                    class_recursive = TRUE, drop_zero_volume = FALSE,
                                    feature_version = NULL, use_cell_counts = FALSE,
                                    single_cell_values = c(-1, 0), carbon_conversion = c("roi", "cell"),
                                    use_python = FALSE, verbose = TRUE,
                                    mat_folder = deprecated(), mat_files = deprecated(), mat_recursive = deprecated()) {

  diatom_equation <- match.arg(diatom_equation)
  carbon_conversion <- match.arg(carbon_conversion)
  check_carbon_conversion(carbon_conversion, use_cell_counts)

  # Handle deprecated mat_folder argument
  if (lifecycle::is_present(mat_folder)) {
    deprecate_warn("0.7.0", "iRfcb::ifcb_extract_biovolumes(mat_folder = )", "iRfcb::ifcb_extract_biovolumes(class_files = )")
    class_files <- mat_folder
  }

  # Handle deprecated mat_files argument
  if (lifecycle::is_present(mat_files)) {
    deprecate_warn("0.8.0", "iRfcb::ifcb_extract_biovolumes(mat_files = )", "iRfcb::ifcb_extract_biovolumes(class_files = )")
    class_files <- mat_files
  }

  # Handle deprecated mat_recursive argument
  if (lifecycle::is_present(mat_recursive)) {
    deprecate_warn("0.8.0", "iRfcb::ifcb_extract_biovolumes(mat_recursive = )", "iRfcb::ifcb_extract_biovolumes(class_recursive = )")
    class_recursive <- mat_recursive
  }

  if (is.null(class_files) && (is.null(custom_images) || is.null(custom_classes))) {
    cli_abort(c(
      "No classification information supplied.",
      "i" = "Provide either {.arg class_files} or both {.arg custom_images} and {.arg custom_classes}."
    ))
  }

  if (!is.null(class_files) && (!is.null(custom_images) || !is.null(custom_classes))) {
    cli_warn(c(
      "Both {.arg class_files} and {.arg custom_images}/{.arg custom_classes} were provided.",
      "i" = "Using {.arg class_files} and ignoring {.arg custom_images}/{.arg custom_classes}."
    ))
  }

  if (is.character(feature_files)) {
    if (length(feature_files) == 1) {
      if (file.exists(feature_files)) {
        # It's a single file
      } else if (dir.exists(feature_files)) {
        # It's a directory
      } else {
        cli_abort("The specified file or directory does not exist: {.file {feature_files}}")
      }
    } else {
      # feature_files is a vector of files
      if (!all(file.exists(feature_files))) {
        missing <- feature_files[!file.exists(feature_files)]
        cli_abort(c(
          "{length(missing)} of {length(feature_files)} {.arg feature_files} do{?es/} not exist:",
          "x" = "{.file {missing}}"
        ))
      }
    }
  } else {
    cli_abort("{.arg feature_files} must be a character vector of filenames or a single directory path.")
  }

  # Check if feature_files is a single folder path or a vector of file paths
  if (length(feature_files) == 1 && dir.exists(feature_files)) {
    feature_files <- list.files(feature_files, pattern = "D.*\\.csv", full.names = TRUE, recursive = feature_recursive)
  }

  if (!is.null(class_files)) {

    # Check if class_files is a single folder path or a vector of file paths
    if (length(class_files) == 1 && file.info(class_files)$isdir) {
      class_files <- list.files(class_files, pattern = "\\.(mat|h5|csv)$", recursive = class_recursive, full.names = TRUE)
      # A directory may hold non-class .csv files (e.g. dashboard class_scores
      # exports); drop them with a warning before any file is read or its sample
      # name and date are parsed.
      class_files <- drop_invalid_class_csv(class_files)
    }

    if (length(class_files) == 0) {
      cli_abort("No classification files found.")
    }

    # Check if files are manually classified (.h5 and .csv files are never manual)
    is_manual <- tolower(tools::file_ext(class_files[1])) == "mat" &&
      "class2use_manual" %in% ifcb_get_mat_names(class_files[1])

    if (is_manual && is.null(class2use_file)) {
      cli_abort("{.arg class2use_file} must be specified when extracting manual biovolume data.")
    }

    if (use_cell_counts && is_manual) {
      cli_abort(c(
        "{.arg use_cell_counts = TRUE} is not supported for manually annotated files.",
        "i" = "Chain-count data is only stored in automated {.file .mat}, {.file .h5} and {.file .csv} classification files."
      ))
    }
  }

  if (use_cell_counts && !is.null(custom_images)) {
    cli_abort(c(
      "{.arg use_cell_counts = TRUE} cannot be combined with {.arg custom_images}/{.arg custom_classes}.",
      "i" = "Chain-count data is read from {.arg class_files} ({.file .h5} or {.file .csv})."
    ))
  }

  if (!is.null(custom_images)) {
    # Extract date-time from class file paths
    class_date_times <- gsub(".*D(\\d{8}T\\d{6})_.*", "\\1", custom_images)
  } else {
    # Extract date-time from class file paths
    class_date_times <- gsub(".*D(\\d{8}T\\d{6})_.*", "\\1", class_files)
  }

  # Extract date-time from feature file paths
  extracted_dates <- sub(".*D(\\d{8}T\\d{6}).*", "\\1", feature_files)

  # List matching feature files
  feature_files <- feature_files[extracted_dates %in% class_date_times]

  # Read feature files
  features <- ifcb_read_features(feature_files = feature_files,
                                 multiblob = multiblob,
                                 feature_version = feature_version,
                                 biovolume_only = TRUE,
                                 verbose = verbose)

  if (length(features) == 0) {
    cli_abort("No feature data files found.")
  }

  data_list <- vector("list", length(features))
  idx <- 1

  for (file_name in names(features)) {

    file_data <- features[[file_name]]

    if (drop_zero_volume) {
      file_data <- file_data[file_data$Biovolume != 0, , drop = FALSE]
    }

    if (nrow(file_data) > 0) {

      sample_name <- if (multiblob) {
        str_replace(file_name, "_multiblob_v\\d+.csv", "")
      } else {
        str_replace(file_name, "_fe[a-z]*_v\\d+\\.csv", "")
      }

      data_list[[idx]] <- tibble(
        sample = sample_name,
        roi_number = file_data$roi_number,
        biovolume = file_data$Biovolume
      )

      idx <- idx + 1

    } else if (drop_zero_volume) {
      cli_warn("All rows were dropped for file {.file {file_name}} because {.code Biovolume == 0}.")
    }
  }

  biovolume_df <- do.call(rbind, data_list[seq_len(idx - 1)])

  # Stop if the combined data frame is empty
  if (is.null(biovolume_df)) {
    cli_abort("No biovolume data available in feature files.")
  }

  # If custom class list is supplied
  if (!is.null(custom_images) && !is.null(custom_classes)) {
    if (length(custom_images) != length(custom_classes)) {
      cli_abort(c(
        "The number of images does not match the number of class labels.",
        "x" = "{.arg custom_images} has length {length(custom_images)}",
        "x" = "{.arg custom_classes} has length {length(custom_classes)}"
      ))
    }

    image_df <- ifcb_convert_filenames(custom_images)
    class_df <- image_df %>%
      mutate(classifier = NA,
             class = custom_classes) %>%
      select(sample, classifier, roi, class) %>%
      rename(roi_number = roi)

  } else {
    # Find unique samples in biovolume_df
    unique_samples <- unique(biovolume_df$sample)

    # Identify matching class files
    class_samples <- sub(".*(D\\d{8}T\\d{6}_IFCB\\d+).*", "\\1", class_files)
    matching_class_files <- class_files[class_samples %in% unique_samples]

    # Initialize an empty list to store data frames
    tb_list <- list()

    if (is_manual) {

      class_df <- ifcb_count_mat_annotations(matching_class_files,
                                             class2use_file,
                                             sum_level = "roi",
                                             use_python = use_python)

      names(class_df)[2] <- "roi_number"

      class_df$classifier <- NA

    } else {
      n_files <- length(matching_class_files)

      # Guard against one sample resolving to more than one classification file
      # (e.g. a .mat, .h5 and .csv for the same sample in one folder). Every row
      # survives the join below, so that sample's ROIs are duplicated and its
      # counts, biovolume and carbon multiply. This reuses `class_samples`
      # instead of re-deriving names from the file extensions, so the guard
      # cannot disagree with the selection above about which sample a file
      # belongs to.
      matched_samples <- class_samples[class_samples %in% unique_samples]
      dup_samples <- unique(matched_samples[duplicated(matched_samples)])
      if (length(dup_samples) > 0) {
        cli_abort(c(
          "{length(dup_samples)} sample{?s} resolve{?s/} to more than one classification file: {.val {dup_samples}}.",
          "i" = "Supply a single file format per sample (e.g. only {.file .mat} or only {.file .h5}) to avoid double-counting."
        ))
      }

      tb_list <- vector("list", n_files)
      has_chain <- logical(n_files)
      na_gaps <- integer(n_files)
      file_samples <- character(n_files)

      # Set up the progress bar
      if (verbose && n_files > 0) {
        cli_progress_bar("Reading classification files", total = n_files)
      }

      for (i in seq_along(matching_class_files)) {

        if (verbose && n_files > 0) {
          cli_progress_update()
        }

        temp <- suppressWarnings({
          read_class_file(matching_class_files[i], use_python = use_python)
        })

        sample_name <- sub("_class(_v\\d+)?\\.(mat|h5)$", "", basename(matching_class_files[i]))
        # Also handle CSV files, stripping a _class(_vN) suffix there too so a
        # label file named {sample}_class.csv joins under the same sample as
        # {sample}.csv (and as the selection regex above already assumes).
        sample_name <- sub("(_class(_v\\d+)?)?\\.csv$", "", sample_name)

        has_chain[i] <- !is.null(temp$cell_count)
        # A missing value inside a file that does carry cell_count data (a
        # blank CSV cell, an HDF5 NaN, a value that failed to parse) nulls the
        # whole sample-class group below; count the gaps so that can be said
        # out loud rather than surface as an unexplained NA.
        na_gaps[i] <- if (use_cell_counts && has_chain[i]) sum(is.na(temp$cell_count)) else 0L
        file_samples[i] <- sample_name

        tb_list[[i]] <- tibble(
          sample = sample_name,
          classifier = temp$classifierName,
          roi_number = temp$roinum,
          class = if (threshold == "opt")
            unlist(temp$TBclass_above_threshold)
          else
            unlist(temp$TBclass),
          cell_count = if (use_cell_counts) {
            if (is.null(temp$cell_count)) NA_integer_ else temp$cell_count
          } else {
            NULL
          }
        )
      }

      # Close the progress bar
      if (verbose && n_files > 0) cli_progress_done()

      if (use_cell_counts) {
        if (!any(has_chain)) {
          cli_abort(c(
            "{.arg use_cell_counts = TRUE} but none of the classification files contain chain-count data.",
            "i" = "Re-run classification with chain counting enabled to produce a {.code cell_count} dataset."
          ))
        }
        # Not gated on `verbose`: this reports a data-integrity condition that
        # changes the returned numbers, not progress.
        if (!all(has_chain)) {
          cli_warn(c(
            "{sum(!has_chain)} of {n_files} classification file{?s} {qty(sum(!has_chain))}{?does/do} not contain chain-count data.",
            "i" = "ROIs from {qty(sum(!has_chain))}{?this file/these files} are treated as {.code NA} chain counts, so {.field cell_counts} is {.code NA} for the affected samples."
          ))
        }
        if (any(na_gaps > 0)) {
          cli_warn(c(
            "{sum(na_gaps)} ROI{?s} in {sum(na_gaps > 0)} classification file{?s} with chain-count data {qty(sum(na_gaps))}{?has/have} a missing {.code cell_count} value.",
            "i" = "{.field cell_counts} is {.code NA} for the affected sample{?s}: {.val {file_samples[na_gaps > 0]}}."
          ))
        }
      }
    }

    if (!is_manual) {
      # Combine all data frames into one
      class_df <- do.call(rbind, tb_list)
    }
  }

  # Join biovolume_df with class_df
  biovolume_df <- left_join(biovolume_df, class_df, by = c("sample", "roi_number"))

  # Calculate biovolume in cubic microns
  biovolume_df$biovolume_um3 <- biovolume_df$biovolume * (micron_factor ^ 3)

  # Determine if each class is a diatom
  unique_classes <- unique(biovolume_df$class)

  if (verbose) {
    cli_inform("Retrieving WoRMS records...")
  }

  diatom_details <- ifcb_is_diatom(unique_classes,
                                   diatom_class = diatom_class,
                                   marine_only = marine_only,
                                   details = TRUE,
                                   verbose = verbose)

  is_diatom <- tibble(class = unique_classes,
                      worms_class = diatom_details$worms_class,
                      is_diatom = diatom_details$is_diatom)

  # Override diatom classification if diatom_include is provided
  if (!is.null(diatom_include)) {
    matched <- is_diatom$class %in% diatom_include
    if (verbose && any(matched)) {
      cli_alert_info(
        "The following {qty(sum(matched))}class{?es} {?was/were} manually included as diatoms via {.arg diatom_include}:"
      )
      cli_inform("{.val {sort(is_diatom$class[matched])}}")
    }
    is_diatom$is_diatom[matched] <- TRUE
  }

  biovolume_df <- left_join(biovolume_df, is_diatom[, c("class", "is_diatom")], by = "class")

  if (verbose) {
    # Classes treated as diatoms (short, verifiable list - shown in full)
    diatoms <- sort(is_diatom$class[is_diatom$is_diatom])

    # Classes that could not be found in WoRMS (the genuinely ambiguous bucket)
    not_found <- sort(is_diatom$class[is.na(is_diatom$worms_class) & !is_diatom$is_diatom])

    # Classes resolved by WoRMS to a non-diatom class
    non_diatoms <- sort(is_diatom$class[!is_diatom$is_diatom & !is.na(is_diatom$worms_class)])

    if (length(diatoms) > 0) {
      cli_alert_info(
        "{length(diatoms)} of {nrow(is_diatom)} {qty(length(diatoms))}class{?es} {?is/are} treated as diatoms:"
      )
      cli_inform("{.val {cli::cli_vec(diatoms, list('vec-trunc' = Inf))}}")
    }

    if (length(non_diatoms) > 0) {
      cli_alert_info(
        paste("{length(non_diatoms)} {qty(length(non_diatoms))}class{?es} {?is/are} treated as NOT diatoms.",
              "To check for genus homonyms (e.g. Navicula, Actinocyclus, which share names with animals),",
              "run {.code ifcb_is_diatom(details = TRUE)} and inspect the {.field worms_class} column.")
      )
    }

    if (length(not_found) > 0) {
      cli_alert_info(
        "{length(not_found)} {qty(length(not_found))}class{?es} could not be found in WoRMS and {?was/were} assumed NOT diatoms:"
      )
      cli_inform("{.val {cli::cli_vec(not_found, list('vec-trunc' = Inf))}}")
    }
  }

  # Resolve per-ROI cell counts (abundance) from chain counts. Computed here,
  # ahead of the carbon block that may need it, but assigned onto the tibble
  # below so that `cell_count_resolved` keeps its position after `carbon_pg`.
  cells <- if (use_cell_counts) {
    resolve_cell_counts(biovolume_df$cell_count, single_cell_values)
  }

  # Select the diatom carbon conversion function
  vol2C_diatom_fun <- switch(diatom_equation,
    all = vol2C_diatom,
    auto = vol2C_diatom_auto,
    vol2C_lgdiatom
  )

  # Under carbon_conversion = "cell" the equations are applied to the per-cell
  # volume and summed back over the chain; under "roi" the wrapper returns the
  # conversion functions untouched, so the numbers are exactly as before.
  carbon_cells <- if (carbon_conversion == "cell") cells
  diatom_fun <- scale_vol2C_per_cell(vol2C_diatom_fun, carbon_cells)
  nondiatom_fun <- scale_vol2C_per_cell(vol2C_nondiatom, carbon_cells)

  # Calculate carbon content based on diatom classification
  biovolume_df <- biovolume_df %>%
    mutate(carbon_pg = case_when(
      !is.na(is_diatom) & is_diatom ~ diatom_fun(biovolume_um3),
      !is.na(is_diatom) & !is_diatom ~ nondiatom_fun(biovolume_um3),
      is.na(is_diatom) ~ nondiatom_fun(biovolume_um3),
      TRUE ~ NA_real_
    )) %>%
    select(-biovolume, -is_diatom)

  if (use_cell_counts) {
    biovolume_df$cell_count_resolved <- cells
  }

  type_convert(biovolume_df, col_types = cols())
}

Try the iRfcb package in your browser

Any scripts or data that you put into this service are public.

iRfcb documentation built on Aug. 20, 2026, 1:06 a.m.