Nothing
#' Plot Dynamic B0
#'
#' Plots the spawning output with and without fishing mortality
#'
#' @template replist
#' @param ylab Y-axis label. Default is "Spawning biomass (mt)" which is replaced
#' by "Spawning output" for models with `replist[["SpawnOutputUnits"]] == "numbers"`
#' @param equilibrium Show equilibrium in plot? Applies whether "yrs" is specified
#' or not.
#' @param forecast Show forecast years in plot? Only applies if yrs = "all".
#' @param yrs Which years to include. Default "all" will show startyr to endyr + 1
#' modified by the arguments `forecast`.
#' @param plot Plot to active plot device?
#' @param print Print to PNG files?
#' @param plotdir Directory where PNG files will be written. By default it will
#' be the directory where the model was run.
#' @template verbose
#' @param uncertainty Show 95% uncertainty intervals around point estimates?
#' These intervals will only appear when uncertainty in the dynamic B0
#' estimates is available via the control file settings for
#' "read specs for more stddev reporting".
#' @template legend
#' @param legendlabels Character vector with labels for the unfished
#' equilibrium point (if `equilibrium = TRUE`) and the two lines showing
#' spawning biomass or output without and with fishing.
#' @template legendloc
#' @param col Optional vector of colors to be used for the two lines
#' (single value will apply to both lines).
#' @param lty Optional vector of line types to be used for the two lines
#' (single value will apply to both lines).
#' @param lwd Optional vector of line widths to be used for the two lines
#' (single value will apply to both lines).
#' @param add add to existing plot
#' @template pwidth
#' @template pheight
#' @template punits
#' @template ptsize
#' @template res
#' @template mainTitle
#' @template mar
#' @author Ian G. Taylor
#' @export
#' @seealso [SSplotTimeseries()]
SSplotDynamicB0 <- function(replist,
ylab = "Spawning biomass (mt)",
equilibrium = TRUE,
forecast = FALSE,
yrs = "all",
plot = TRUE, print = FALSE, plotdir = "default", verbose = TRUE,
uncertainty = TRUE,
legend = TRUE,
legendlabels = c("equilibrium", "without fishing", "with fishing"),
legendloc = "bottom",
col = c("blue", "red"),
lty = 1,
lwd = 2,
add = FALSE,
pwidth = 6.5, pheight = 5.0, punits = "in", res = 300, ptsize = 10,
mainTitle = FALSE,
mar = NULL) {
Dynamic_Bzero <- replist[["Dynamic_Bzero"]]
if (is.null(Dynamic_Bzero)) {
warning('No element "Dynamic_Bzero" in replist input')
return()
}
if (!"SSB_nofishing" %in% names(replist[["Dynamic_Bzero"]])) {
warning("Skipping Dynamic B0 plot (not yet working for this model configuration)")
return()
}
# check if spawning output rather than spawning biomass is plotted
if (is.null(replist[["SpawnOutputUnits"]]) ||
is.na(replist[["SpawnOutputUnits"]]) ||
replist[["SpawnOutputUnits"]] == "numbers") { # quantity from test in SS_output
if (ylab == "Spawning biomass (mt)") {
ylab <- "Spawning output"
}
}
# set default plot margins
if (is.null(mar)) {
if (mainTitle) {
main <- "Dynamic B0"
mar <- c(5, 4, 4, 2) + 0.1
} else {
main <- ""
mar <- c(5, 4, 2, 2) + 0.1
}
}
# repeat inputs for line style if needed
if (length(col) == 1) {
col <- rep(col, 2)
}
if (length(lty) == 1) {
lty <- rep(lty, 2)
}
if (length(lwd) == 1) {
lwd <- rep(lwd, 2)
}
plotinfo <- NULL
# set plot directory
if (plotdir == "default") {
plotdir <- replist[["inputs"]][["dir"]]
}
# set default range of years
if (yrs[1] == "all") {
yrs <- replist[["startyr"]]:(replist[["endyr"]] + 1)
if (forecast) {
yrs <- replist[["startyr"]]:max(Dynamic_Bzero[["Yr"]])
}
}
# equilibrium year shifted by 1 because INIT point isn't shown
sub_equil <- Dynamic_Bzero[["Era"]] == "VIRG"
equil_yr <- Dynamic_Bzero[["Yr"]][sub_equil] + 1
# set ymax
ymax <- max(ifelse(equilibrium,
Dynamic_Bzero[["SSB"]][sub_equil],
NA
),
Dynamic_Bzero[["SSB"]][Dynamic_Bzero[["Yr"]] %in% yrs],
Dynamic_Bzero[["SSB_nofishing"]][Dynamic_Bzero[["Yr"]] %in% yrs],
na.rm = TRUE
)
# get subsets of derived quantities and extract year value
quants <- replist[["derived_quants"]]
quants_SSB <- quants[grep("SSB_", quants[["Label"]]), ]
quants_SSB_nofishing <- quants[grep("Dyn_Bzero_", quants[["Label"]]), ]
if (uncertainty & nrow(quants_SSB_nofishing) == 0) {
uncertainty <- FALSE
warning(
"Dynamic B0 not found in derived quantities, ",
"changing uncertainty to FALSE. ",
"To get uncertainty, modify control file under ",
"'read specs for more stddev reporting'."
)
}
if (uncertainty) {
# calculate intervals
# add columns to store info
Dynamic_Bzero[["SSB_lo"]] <- NA
Dynamic_Bzero[["SSB_hi"]] <- NA
Dynamic_Bzero[["SSB_nofishing_lo"]] <- NA
Dynamic_Bzero[["SSB_nofishing_hi"]] <- NA
# get year for each row
quants_SSB[["Yr"]] <- substring(quants_SSB[["Label"]],
first = nchar("SSB_") + 1
)
quants_SSB_nofishing[["Yr"]] <- substring(quants_SSB_nofishing[["Label"]],
first = nchar("Dyn_Bzero_") + 1
)
# get unfished equilibrium uncertainty intervals (same with/without fishing)
Dynamic_Bzero[sub_equil, c("SSB_lo", "SSB_hi")] <-
qnorm(
p = c(0.025, 0.975),
mean = quants["SSB_Virgin", "Value"],
sd = quants["SSB_Virgin", "StdDev"]
)
# loop over years to get uncertainty intervals
for (y in Dynamic_Bzero[["Yr"]]) {
# fill in SSB uncertainty
if (y %in% quants_SSB[["Yr"]]) {
Dynamic_Bzero[Dynamic_Bzero[["Yr"]] == y, c("SSB_lo", "SSB_hi")] <-
qnorm(
p = c(0.025, 0.975),
mean = quants_SSB[paste0("SSB_", y), "Value"],
sd = quants_SSB[paste0("SSB_", y), "StdDev"]
)
}
# fill in dynamic B0 uncertainty
if (y %in% quants_SSB_nofishing[["Yr"]]) {
Dynamic_Bzero[Dynamic_Bzero[["Yr"]] == y, c("SSB_nofishing_lo", "SSB_nofishing_hi")] <-
qnorm(
p = c(0.025, 0.975),
mean = quants_SSB_nofishing[paste0("Dyn_Bzero_", y), "Value"],
sd = quants_SSB_nofishing[paste0("Dyn_Bzero_", y), "StdDev"]
)
}
}
# update ymax
ymax <- max(ymax,
ifelse(equilibrium,
Dynamic_Bzero[["SSB_hi"]][sub_equil],
NA
),
Dynamic_Bzero[["SSB_hi"]][Dynamic_Bzero[["Yr"]] %in% yrs],
Dynamic_Bzero[["SSB_nofishing_hi"]][Dynamic_Bzero[["Yr"]] %in% yrs],
na.rm = TRUE
)
}
plotfun <- function() {
if (!add) {
plot(0,
type = "n", , xlab = "Year", ylab = ylab, xlim = range(yrs),
ylim = c(0, 1.1 * ymax),
yaxs = "i",
cex.lab = 1.0, cex.axis = 1.0, cex = 0.7,
mar = mar,
main = main
)
}
sub <- Dynamic_Bzero[["Yr"]] %in% yrs & Dynamic_Bzero[["Era"]] != "VIRG"
if (uncertainty) {
# add intervals to plot
if (equilibrium) {
# add equilibrium intervals
arrows(
x0 = equil_yr,
y0 = Dynamic_Bzero[sub_equil, "SSB_lo"],
x1 = equil_yr,
y1 = Dynamic_Bzero[sub_equil, "SSB_hi"],
length = 0.01, angle = 90, code = 3, col = col[1],
lwd = 1
)
}
# function to add shaded uncertainty intervals behind line
addpoly <- function(yrvec, lower, upper, col) {
lower[lower < 0] <- 0 # max of value or 0
good <- !is.na(lower) & !is.na(upper)
polygon(
x = c(yrvec[good], rev(yrvec[good])),
y = c(lower[good], rev(upper[good])),
border = NA, col = col
)
}
# add polygons around time series
addpoly(
yrvec = Dynamic_Bzero[["Yr"]][sub],
lower = Dynamic_Bzero[["SSB_lo"]][sub],
upper = Dynamic_Bzero[["SSB_hi"]][sub],
col = adjustcolor(col[2], alpha.f = 0.1)
)
addpoly(
yrvec = Dynamic_Bzero[["Yr"]][sub],
lower = Dynamic_Bzero[["SSB_nofishing_lo"]][sub],
upper = Dynamic_Bzero[["SSB_nofishing_hi"]][sub],
col = adjustcolor(col[1], alpha.f = 0.1)
)
} # end check for uncertainty
lines(
x = Dynamic_Bzero[["Yr"]][sub],
y = Dynamic_Bzero[["SSB"]][sub],
lwd = lwd[2], lty = lty[2], col = col[2]
)
lines(
x = Dynamic_Bzero[["Yr"]][sub],
y = Dynamic_Bzero[["SSB_nofishing"]][sub],
lwd = lwd[1], lty = lty[1], col = col[1]
)
if (equilibrium) {
points(
x = equil_yr,
y = Dynamic_Bzero[["SSB"]][sub_equil],
cex = 1.5,
pch = 16,
col = col[1]
)
# add to vectors used in the legend
lty <- c(NA, lty)
lwd <- c(NA, lwd)
col <- c(col[1], col)
pch <- c(16, NA, NA)
} else {
pch <- NA
legendlabels <- legendlabels[2:3]
}
if (legend) {
legend(
x = legendloc,
legend = legendlabels,
col = col,
lwd = lwd,
lty = lty,
pch = pch,
pt.cex = 1.5,
bty = "n",
ncol = ifelse(equilibrium, 3, 2)
)
}
}
if (plot) {
plotfun()
}
if (print) {
caption <- paste0(
"Dynamic B0 plot. The lower line shows the time series ",
"of estimated ", ylab, " in the presence of fishing ",
"mortality. The upper line shows the time series that ",
"could occur under the same dynamics (including ",
"deviations in recruitment), but without fishing."
)
if (equilibrium) {
caption <- paste0(
caption,
" The point at the left represents the unfished equilibrium."
)
}
file <- "ts_DynamicB0.png"
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
plotfun()
dev.off()
if (!is.null(plotinfo)) {
plotinfo[["category"]] <- "Timeseries"
}
}
if (verbose) message("Plotting Dynamic B0")
return(invisible(plotinfo))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.