R/SS_ForeCatch.R

In r4ss/r4ss: R Code for Stock Synthesis

#' Create table of fixed forecast catches
#'
#' Processing values of dead or retained biomass from timeseries output
#' to fit the format required at the bottom of the forecast file.
#' This can be used to map the catches resulting from forecasting with a particular
#' harvest control rule into a model representing a different state of nature.
#' This is a common task for US west coast groundfish but might be useful elsewhere.
#'
#' @template replist
#' @param yrs Range of years in which to fill in forecast catches from timeseries
#' @param average Use average catch over a range of years for forecast
#' (as opposed to using forecast based on control rule)
#' @param avg.yrs Range of years to average over
#' @param total Either single value or vector of annual total forecast catch
#' used to scale values (especially if values are from average catches).
#' For west coast groundfish, total might be ACL for next 2 forecast years
#' @param digits Number of digits to round to in table
#' @param dead TRUE/FALSE switch to choose dead catch instead of retained catch.
#' @param zeros Include entries with zero catch (TRUE/FALSE)
#'
#' @seealso [SS_readforecast()], [SS_readforecast()]
#' @author Ian G. Taylor
#' @examples
#' \dontrun{
#' # create table based on average over past 5 years
#' SS_ForeCatch(base, # object created by SS_output
#'   yrs = 2021:2022, # years with fixed catch
#'   average = TRUE, # catch by fleet from average catch
#'   avg.yrs = 2014:2018
#' ) # use average of catches over past 5 years
#'
#' # create table with pre-defined totals where the first 2 years
#' # are based on current harvest specifications and the next 10 are set to some
#' # new value (with ratio among fleets based on average over past 5 years)
#' SS_ForeCatch(base, # object created by SS_output
#'   yrs = 2021:2022, # years with fixed catch
#'   average = TRUE, # catch by fleet from average catch
#'   avg.yrs = 2016:2020, # use average of catches over past 5 years
#'   total = c(rep(241.3, 2), rep(300, 10))
#' ) # total
#'
#' # create table based on harvest control rule projection in SS
#' # that can be mapped into an alternative state of nature
#' SS_ForeCatch(low_state, # object created by SS_output for low state
#'   yrs = 2021:2032, # forecast period after fixed ACL years
#'   average = FALSE
#' ) # use values forecast in SS, not historic catch
#' }
#'
#' @export

SS_ForeCatch <- function(replist, yrs = 2021:2032,
                         average = FALSE, avg.yrs = 2016:2020,
                         total = NULL, digits = 2,
                         dead = TRUE, zeros = FALSE) {
  # function for creating table of fixed forecast catches
  # based on values in the timeseries output
  timeseries <- replist[["timeseries"]]

  # create new empty object to store stuff
  forecast_catches <- NULL

  if (!all(yrs %in% timeseries[["Yr"]])) {
    warning("Not all requested years are present in timeseries output.")
    yrs <- yrs[yrs %in% timeseries[["Yr"]]]
  }
  # if only one value for total is input, repeat for all years
  if (!is.null(total) & length(total) == 1) {
    total <- rep(total, length(yrs))
  }
  # loop over seasons and areas (only 1 for most models) to subset timeseries
  for (iyr in seq_along(yrs)) {
    y <- yrs[iyr]
    forecast_catches_y <- NULL
    for (iseas in 1:replist[["nseasons"]]) {
      for (iarea in 1:replist[["nareas"]]) {
        if (replist[["SS_versionNumeric"]] < 3.30) {
          fleets_with_catch <- 1:replist[["nfishfleets"]]
        } else {
          fleets_with_catch <- which(replist[["fleet_type"]] == 1)
        }
        for (ifleet in fleets_with_catch) {
          # figure out column name
          if (replist[["catch_units"]][ifleet] == 1) {
            string <- ifelse(dead, "dead(B)", "retain(B):_")
          }
          if (replist[["catch_units"]][ifleet] == 2) {
            string <- ifelse(dead, "dead(N)", "retain(N):_")
          }
          colname <- paste0(string, ":_", ifleet)
          # extract catch
          if (average) {
            catches <- timeseries[
              timeseries[["Area"]] == iarea &
                timeseries[["Seas"]] == iseas &
                timeseries[["Yr"]] %in% avg.yrs,
              colname
            ]
            catch <- mean(catches)
          } else {
            catch <- timeseries[
              timeseries[["Area"]] == iarea &
                timeseries[["Seas"]] == iseas &
                timeseries[["Yr"]] == y,
              colname
            ]
          }
          # create new row for table
          newrow <- data.frame(
            Year = y, Seas = iseas,
            Fleet = ifleet, Catch = catch
          )
          # add new row to previous rows
          forecast_catches_y <- rbind(forecast_catches_y, newrow)
        } # end loop over fleets
      } # end loop over areas
    } # end loop over seasons

    # if requested, scale catches to sum to input total (such as ACL)
    if (!is.null(total)) {
      forecast_catches_y[["Catch"]] <- total[iyr] *
        forecast_catches_y[["Catch"]] / sum(forecast_catches_y[["Catch"]])
    }

    # round values
    forecast_catches_y[["Catch"]] <- round(forecast_catches_y[["Catch"]], digits)

    # add comment on right-hand-side
    forecast_catches_y[["comment"]] <- ""
    forecast_catches_y[["comment"]][1] <- paste0(
      "#sum_for_", y, ": ",
      sum(forecast_catches_y[["Catch"]])
    )
    # add block for this year to other blocks
    forecast_catches <- rbind(forecast_catches, forecast_catches_y)
  } # end loop over years

  # fix up column names
  names(forecast_catches)[1] <- "#Year"
  names(forecast_catches)[4] <- string

  # remove rows with zero catch
  if (!zeros & any(forecast_catches[[string]] == 0)) {
    forecast_catches <- forecast_catches[!forecast_catches[[string]] == 0, ]
  }

  return(forecast_catches)
}