R/SSplotYield.R

In r4ss: R Code for Stock Synthesis

#' Plot yield and surplus production.
#'
#' Plot yield and surplus production from Stock Synthesis output. Surplus
#' production is based on Walters et al. (2008).
#'
#'
#' @template replist
#' @param subplots vector controlling which subplots to create
#' Numbering of subplots is as follows:
#' \itemize{
#'   \item 1 yield curve
#'   \item 2 yield curve with reference points
#'   \item 3 surplus production vs. biomass plots (Walters et al. 2008)
#' }
#' @param refpoints character vector of which reference points to display in
#' subplot 2, from the options 'MSY', 'Btgt', and 'SPR'.
#' @param add add to existing plot? (not yet implemented)
#' @param plot plot to active plot device?
#' @param print print to PNG files?
#' @param labels vector of labels for plots (titles and axis labels)
#' @param col line color for equilibrium plot
#' @param col2 line color for dynamic surplus production plot
#' @param lty line type (only applied to equilibrium yield plot at this time)
#' @param lwd line width (only applied to equilibrium yield plot at this time)
#' @param cex.main character expansion for plot titles
#' @param pwidth width of plot
#' @param pheight height of plot
#' @param punits units for PNG file
#' @template res
#' @param ptsize point size for PNG file
#' @param plotdir directory where PNG files will be written. by default it will
#' be the directory where the model was run.
#' @param verbose report progress to R GUI?
#' @author Ian Stewart, Ian Taylor
#' @export
#' @seealso [SS_plots()], [SS_output()]
#' @references Walters, Hilborn, and Christensen, 2008, Surplus production
#' dynamics in declining and recovering fish populations.  Can. J. Fish. Aquat.
#' Sci. 65: 2536-2551
SSplotYield <-
  function(replist,
           subplots = 1:4,
           refpoints = c("MSY", "Btgt", "SPR", "Current"),
           add = FALSE, plot = TRUE, print = FALSE,
           labels = c(
             "Fraction unfished", # 1
             "Equilibrium yield (mt)", # 2
             "Total biomass (mt)", # 3
             "Surplus production (mt)", # 4
             "Yield per recruit (kg)" # 5
           ),
           col = "blue", col2 = "black", lty = 1, lwd = 2, cex.main = 1,
           pwidth = 6.5, pheight = 5.0, punits = "in", res = 300, ptsize = 10,
           plotdir = "default",
           verbose = TRUE) {
    # table to store information on each plot
    plotinfo <- NULL

    # extract quantities from replist
    equil_yield <- replist[["equil_yield"]]
    # column named changed from Catch to Tot_Catch in SSv3.30
    if ("Tot_Catch" %in% names(equil_yield)) {
      equil_yield[["Catch"]] <- equil_yield[["Tot_Catch"]]
    }
    nareas <- replist[["nareas"]]
    nseasons <- replist[["nseasons"]]
    timeseries <- replist[["timeseries"]]
    SSB0 <- replist[["derived_quants"]]["SSB_Virgin", "Value"]

    # function for yield curve
    yieldfunc <- function(refpoints = NULL) {
      if (!add) {
        # empty plot
        plot(0,
          type = "n", xlim = c(0, max(equil_yield[["Depletion"]], 1, na.rm = TRUE)),
          ylim = c(0, max(equil_yield[["Catch"]], na.rm = TRUE)),
          xlab = labels[1], ylab = labels[2]
        )
        abline(h = 0, col = "grey")
        abline(v = 0, col = "grey")
      }

      # add lines for reference points (if requested)
      lines(equil_yield[["Depletion"]], equil_yield[["Catch"]],
        lwd = lwd, col = col, lty = lty
      )
      colvec <- c(4, 2, 3, 1)
      if ("MSY" %in% refpoints) {
        lines(
          x = rep(replist[["derived_quants"]]["SSB_MSY", "Value"] / SSB0, 2),
          y = c(0, replist[["derived_quants"]]["Dead_Catch_MSY", "Value"]),
          col = colvec[1], lwd = 2, lty = 2
        )
      }
      if ("Btgt" %in% refpoints) {
        lines(
          x = rep(replist[["derived_quants"]]["SSB_Btgt", "Value"] / SSB0, 2),
          y = c(0, replist[["derived_quants"]]["Dead_Catch_Btgt", "Value"]),
          col = colvec[2], lwd = 2, lty = 2
        )
      }
      if ("SPR" %in% refpoints) {
        lines(
          x = rep(replist[["derived_quants"]]["SSB_SPR", "Value"] / SSB0, 2),
          y = c(0, replist[["derived_quants"]]["Dead_Catch_SPR", "Value"]),
          col = colvec[3], lwd = 2, lty = 2
        )
      }
      if ("Current" %in% refpoints) {
        which_val <- which(abs(equil_yield[["Depletion"]] - replist[["current_depletion"]]) ==
          min(abs(equil_yield[["Depletion"]] - replist[["current_depletion"]])))[1]
        lines(
          x = rep(replist[["current_depletion"]], 2),
          y = c(0, equil_yield[["Catch"]][which_val]),
          col = colvec[4], lwd = 2, lty = 2
        )
      }
      # legend
      which_lines <- c(
        "MSY" %in% refpoints,
        "Btgt" %in% refpoints,
        "SPR" %in% refpoints,
        "Current" %in% refpoints
      )
      if (any(which_lines)) {
        legend("topright",
          bty = "n", lwd = 2, lty = 2,
          col = colvec[which_lines],
          legend = c("MSY", "B target", "SPR target", "Current")
        )
      }
    }

    if (1 %in% subplots | 2 %in% subplots) {
      # test if data is available
      if (!is.null(equil_yield[[1]][1]) && any(!is.na(equil_yield[[1]]))) {
        # further test for bad values
        # (not sure the circumstances where this is needed)
        if (any(!is.na(equil_yield[["Depletion"]])) &
          any(!is.na(equil_yield[["Catch"]])) &
          any(!is.infinite(equil_yield[["Depletion"]]))) {
          if (1 %in% subplots) {
            # make plot
            if (plot) {
              yieldfunc()
            }
            if (print) {
              file <- "yield1_yield_curve.png"
              caption <- "Yield curve"
              plotinfo <- save_png(
                plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                caption = caption
              )
              yieldfunc()
              dev.off()
            }
          }
          if (2 %in% subplots & !is.null(refpoints)) {
            # make plot
            if (plot) {
              yieldfunc(refpoints = refpoints)
            }
            if (print) {
              file <- "yield2_yield_curve_with_refpoints.png"
              caption <- "Yield curve with reference points"
              plotinfo <- save_png(
                plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                caption = caption
              )
              yieldfunc(refpoints = refpoints)
              dev.off()
            }
          }
        } else {
          cat("Skipped equilibrium yield plots: equil_yield has all NA values\n")
        }
      } else {
        cat("Skipped equilibrium yield plots: no equil_yield results in this model\n")
      }
    }

    # timeseries excluding equilibrium conditions or forecasts
    ts <- timeseries[!timeseries[["Era"]] %in% c("VIRG", "FORE"), ]

    # get total dead catch
    stringB <- "dead(B)"
    catchmat <- as.matrix(ts[, substr(names(ts), 1, nchar(stringB)) == stringB])
    # aggregate catch across fleets
    catch <- rowSums(catchmat)

    # aggregate catch and biomass across seasons and areas
    catch_agg <- aggregate(x = catch, by = list(ts[["Yr"]]), FUN = sum)$x
    Bio_agg <- aggregate(x = ts[["Bio_all"]], by = list(ts[["Yr"]]), FUN = sum)$x

    # number of years to consider
    Nyrs <- length(Bio_agg)

    # function to calculate and plot surplus production
    sprodfunc <- function() {
      sprod <- rep(NA, Nyrs)
      # calculate surplus production as difference in biomass adjusted for catch
      sprod[1:(Nyrs - 1)] <- Bio_agg[2:Nyrs] - Bio_agg[1:(Nyrs - 1)] + catch_agg[1:(Nyrs - 1)]
      sprodgood <- !is.na(sprod)
      Bio_agg_good <- Bio_agg[sprodgood]
      sprod_good <- sprod[sprodgood]
      xlim <- c(0, max(Bio_agg_good, na.rm = TRUE))
      ylim <- c(min(0, sprod_good, na.rm = TRUE), max(sprod_good, na.rm = TRUE))
      # make empty plot
      if (!add) {
        plot(0, ylim = ylim, xlim = xlim, xlab = labels[3], ylab = labels[4], type = "n")
      }
      # add lines
      lines(Bio_agg_good, sprod_good, col = col2)
      # make arrows
      old_warn <- options()$warn # previous setting
      options(warn = -1) # turn off "zero-length arrow" warning
      s <- seq(length(sprod_good) - 1)
      arrows(Bio_agg_good[s], sprod_good[s], Bio_agg_good[s + 1], sprod_good[s + 1],
        length = 0.06, angle = 20, col = col2, lwd = 1.2
      )
      options(warn = old_warn) # returning to old value

      # add lines at 0 and 0
      abline(h = 0, col = "grey")
      abline(v = 0, col = "grey")
      # add blue point at start
      points(Bio_agg_good[1], sprod_good[1], col = col2, bg = "white", pch = 21)
    } # end sprodfunc

    # function to plot time series of Yield per recruit
    YPR_timeseries <- function() {
      # exclude forecast years
      if ("Era" %in% names(sprseries)) {
        sub <- sprseries[["Era"]] != "FORE"
      } else {
        # older versions of SS didn't include the Era column
        sub <- sprseries[["Yr"]] <= replist[["endyr"]]
      }
      # plot a line
      plot(
        x = sprseries[["Yr"]][sub],
        y = sprseries[["YPR"]][sub],
        ylim = c(0, 1.1 * max(sprseries[["YPR"]][sub], na.rm = TRUE)),
        xlab = "Year",
        ylab = labels[5],
        type = "l",
        lwd = 2,
        col = "blue",
        yaxs = "i"
      )
    } # end YPR_timeseries function

    if (3 %in% subplots) {
      if (plot) {
        sprodfunc()
      }
      if (print) {
        file <- "yield3_surplus_production.png"
        caption <-
          paste(
            "Surplus production vs. biomass plot. For interpretation, see<br>",
            "<blockquote>Walters, Hilborn, and  Christensen, 2008,",
            "Surplus production dynamics in declining and",
            "recovering fish populations. <i>Can. J. Fish. Aquat. Sci.</i>",
            "65: 2536-2551.</blockquote>"
          )
        plotinfo <- save_png(
          plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
          pheight = pheight, punits = punits, res = res, ptsize = ptsize,
          caption = caption
        )
        sprodfunc()
        dev.off()
      }
    }

    if (4 %in% subplots) {
      # stuff related to subplot 4 (YPR timeseries)
      sprseries <- replist[["sprseries"]]
      if (is.null(sprseries)) {
        if (verbose) {
          message("Skipping yield per recruit plot because SPR_SERIES not in output")
        }
      } else {
        if (plot) {
          YPR_timeseries()
        }
        if (print) {
          file <- "yield4_YPR_timeseries.png"
          caption <- "Time series of yield per recruit (kg)"
          plotinfo <- save_png(
            plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
            pheight = pheight, punits = punits, res = res, ptsize = ptsize,
            caption = caption
          )
          YPR_timeseries()
          dev.off()
        }
      } # end check for sprseries available
    } # end check for 4 in subplots

    if (!is.null(plotinfo)) plotinfo[["category"]] <- "Yield"
    return(invisible(plotinfo))
  } # end function