R/SSplotMnwt.R

In r4ss/r4ss: R Code for Stock Synthesis

#' Plot mean weight data and fits.
#'
#' Plot mean weight data and fits from Stock Synthesis output. Intervals are
#' based on T-distributions as specified in model.
#'
#'
#' @template replist
#' @param ymax Optional input to override default ymax value.
#' @param subplots Vector of which plots to make (1 = data only, 2 = with fit).
#' If `plotdat = FALSE` then subplot 1 is not created, regardless of
#' choice of `subplots`.
#' @template plot
#' @template print
#' @template plotdir
#' @template fleets
#' @template fleetnames
#' @param datplot Make data-only plot of discards? This can override the choice
#' of `subplots`.
#' @template labels
#' @param col1 first color to use in plot (for expected values)
#' @param col2 second color to use in plot (for observations and intervals)
#' @template pwidth
#' @template pheight
#' @template punits
#' @template res
#' @template ptsize
#' @template cex.main
#' @template verbose
#' @author Ian Taylor, Ian Stewart
#' @export
#' @seealso [SS_plots()], [SS_output()]
SSplotMnwt <-
  function(replist, subplots = 1:2, ymax = NULL,
           plot = TRUE, print = FALSE,
           fleets = "all",
           fleetnames = "default",
           datplot = FALSE,
           labels = c(
             "Year", # 1
             "discard", # 2
             "retained catch", # 3
             "whole catch", # 4
             "Mean individual body weight (kg)", # 5
             "Mean weight in", # 6
             "for"
           ), # 7
           col1 = "blue", col2 = "black",
           pwidth = 6.5, pheight = 5.0, punits = "in", res = 300, ptsize = 10,
           cex.main = 1,
           plotdir = "default", verbose = TRUE) {
    # table to store information on each plot
    plotinfo <- NULL

    # get stuff from replist
    mnwgt <- replist[["mnwgt"]]
    FleetNames <- replist[["FleetNames"]]
    DF_mnwgt <- replist[["DF_mnwgt"]]
    nfleets <- replist[["nfleets"]]
    SS_versionshort <- replist[["SS_versionshort"]]


    if (fleets[1] == "all") fleets <- 1:nfleets
    if (fleetnames[1] == "default") fleetnames <- FleetNames
    if (plotdir == "default") plotdir <- replist[["inputs"]][["dir"]]

    # mean body weight observations ###
    if (!is.na(mnwgt)[[1]][1]) {
      for (ifleet in intersect(fleets, unique(mnwgt[["Fleet"]]))) {
        # usemnwgt is subset of mnwgt for the particular fleet
        usemnwgt <- mnwgt[mnwgt[["Fleet"]] == ifleet & mnwgt[["Obs"]] > 0, ]
        if (SS_versionshort == "3.30") {
          usemnwgt[["Part"]] <- usemnwgt[["Part"]]
        } else {
          usemnwgt[["Part"]] <- usemnwgt[["Mkt"]]
        }
        FleetName <- fleetnames[ifleet]
        for (j in unique(usemnwgt[["Part"]])) {
          yr <- usemnwgt[["Yr"]][usemnwgt[["Part"]] == j]
          ob <- usemnwgt[["Obs"]][usemnwgt[["Part"]] == j]
          cv <- usemnwgt[["CV"]][usemnwgt[["Part"]] == j]
          ex <- usemnwgt[["Exp"]][usemnwgt[["Part"]] == j]
          xmin <- min(yr) - 3
          xmax <- max(yr) + 3
          liw <- -ob * cv * qt(0.025, DF_mnwgt) # quantile of t-distribution
          uiw <- ob * cv * qt(0.975, DF_mnwgt) # quantile of t-distribution
          liw[(ob - liw) < 0] <- ob[(ob - liw) < 0] # no negative limits
          titlepart <- labels[2]
          if (j == 2) titlepart <- labels[3]
          if (j == 0) titlepart <- labels[4]
          ptitle <- paste(labels[6], titlepart, labels[7], FleetName, sep = " ")
          ylab <- labels[5]

          # wrap up plot command in function
          bdywtfunc <- function(addfit) {
            plotCI(
              x = yr, y = ob, uiw = uiw, liw = liw, xlab = labels[1], main = ptitle,
              ylo = 0, col = col2, sfrac = 0.005, ylab = ylab, lty = 1, pch = 21, bg = "white",
              xlim = c(xmin, xmax), cex.main = cex.main, ymax = ymax
            )
            abline(h = 0, col = "grey")
            if (addfit) points(yr, ex, col = col1, cex = 2, pch = "-")
          }

          # make plots
          if (!datplot) {
            subplots <- setdiff(subplots, 1) # don't do subplot 1 if datplot=FALSE
          }
          for (isubplot in subplots) { # loop over subplots (data only or with fit)
            if (isubplot == 1) {
              addfit <- FALSE
            } else {
              addfit <- TRUE
            }
            if (plot) {
              bdywtfunc(addfit = addfit)
            }
            if (print) {
              if (!addfit) {
                file <- paste0("bodywt_data_flt", FleetName, ".png")
              } else {
                file <- paste0("bodywt_fit_flt", FleetName, ".png")
              }
              caption <- ptitle
              plotinfo <- save_png(
                plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                caption = caption
              )
              bdywtfunc(addfit = addfit)
              dev.off()
            }
          } # end loop over subplots
        } # end loop over market categories
      } # end loop over fleets
    }
    if (!is.null(plotinfo)) plotinfo[["category"]] <- "Mnwt"
    return(invisible(plotinfo))
  }