Nothing
#' Plot spawner-recruit curve.
#'
#' Plot spawner-recruit curve based on output from Stock Synthesis model.
#'
#'
#' @template replist
#' @param subplot vector of which subplots to show. 1=plot without labels,
#' 2=plot with year labels.
#' @param add add to existing plot?
#' @param plot plot to active plot device?
#' @param print print to PNG files?
#' @param xlim optional control of x range
#' @param ylim optional control of y range
#' @param labels vector containing x-axis label for models with spawning biomass
#' in metric tons, y-axis label, and alternative x-axis for models with a fecundity
#' relationship making spawning output not equal to spawning biomass.
#' @param bioscale multiplier on spawning biomass, set to 0.5 for single-sex
#' models
#' @param plotdir directory where PNG files will be written. by default it will
#' be the directory where the model was run.
#' @template pwidth
#' @template pheight
#' @template punits
#' @template ptsize
#' @template res
#' @param verbose report progress to R GUI?
#' @param colvec vector of length 4 with colors for 3 lines and 1 set of points
#' (where the 4th value for the points is the color of the circle around the
#' background color provided by `ptcol`
#' @param ltyvec vector of length 4 with line types for the 3 lines and 1 set
#' of points, where the points are disconnected (lty=NA) by default
#' @param ptcol vector or single value for the color of the points, "default"
#' will by replaced by a vector of colors of length equal to
#' `nrow(replist[["recruit"]])`
#' @param legend add a legend to the figure?
#' @param legendloc location of legend. By default it is chosen as the first
#' value in the set of "topleft", "topright", "bottomright" that results in no
#' overlap with the points in the plot, but the user can override this with their
#' choice of location. See ?legend for more info on the options.
#' @param minyr minimum year of recruitment deviation to show in plot
#' @param textmindev minimum recruitment deviation for label to be added so
#' only extreme devs are labeled (labels are added to first and last years as
#' well). Default=0.7.
#' @param relative scale both axes so that B0 and R0 are at 1
#' to show spawning output and recruitment relative to the equilibrium
#' @param expected show line for expected recruitment (stock-recruit curve)
#' @param estimated show points for estimated recruitment values
#' (including deviations)
#' @param bias_adjusted show lines for bias adjusted expected recruitment
#' @param show_env add line for expected recruitment with environmental
#' variability
#' @param virg add point for equilibrium conditions (x=B0,y=R0)
#' @param init add point for initial conditions (x=B1,y=R1), only appears
#' if this point differs from virgin values
#' @param forecast include forecast years in the curve?
#' @author Ian Stewart, Ian Taylor
#' @export
#' @seealso [SS_plots()], [SS_output()]
SSplotSpawnrecruit <-
function(replist, subplot = 1:3, add = FALSE, plot = TRUE, print = FALSE,
xlim = NULL, ylim = NULL,
labels = c(
"Spawning biomass (mt)",
"Recruitment (1,000s)",
"Spawning output",
expression(paste("Spawning output (relative to ", italic(B)[0], ")")),
expression(paste("Recruitment (relative to ", italic(R)[0], ")")),
"Log recruitment deviation"
),
bioscale = "default",
plotdir = "default",
pwidth = 6.5, pheight = 6.5, punits = "in", res = 300, ptsize = 10,
verbose = TRUE,
colvec = c("blue", "black", "black", gray(0, 0.7)),
ltyvec = c(1, 2, 1, NA),
ptcol = "default",
legend = TRUE, legendloc = NULL,
# line1="blue",line2="green3",line3="black",ptcol="red",
minyr = "default", textmindev = 0.5, relative = FALSE,
expected = TRUE, estimated = TRUE, bias_adjusted = TRUE,
show_env = TRUE, virg = TRUE, init = TRUE, forecast = FALSE) {
# plot of spawner recruit curve
# table to store information on each plot
plotinfo <- NULL
recruit <- replist[["recruit"]]
if (is.null(recruit)) {
message("Skipping stock-recruit plots: no recruitment information available")
return()
} else {
if (3 %in% subplot && max(abs(recruit[["dev"]]), na.rm = TRUE) < 1e-6) {
subplot <- setdiff(subplot, 3)
}
}
nsexes <- replist[["nsexes"]]
# set axis labels
xlab <- labels[1]
ylab <- labels[2]
# check if spawning output rather than spawning biomass is plotted
if (is.na(replist[["SpawnOutputUnits"]]) ||
replist[["SpawnOutputUnits"]] == "numbers") { # quantity from test in SS_output
xlab <- labels[3]
}
if (relative) {
xlab <- labels[4]
ylab <- labels[5]
}
# scaling factor for single sex models
if (bioscale == "default") {
if (nsexes == 1) bioscale <- 0.5 else bioscale <- 1
}
if (plotdir == "default") plotdir <- replist[["inputs"]][["dir"]]
if (minyr == "default") minyr <- min(recruit[["Yr"]])
recruit <- recruit[recruit[["era"]] %in% c(
"Early", "Main", "Fixed", "Late",
ifelse(forecast, "Forecast", NA)
) &
recruit[["Yr"]] >= minyr, ]
timeseries <- replist[["timeseries"]]
recruit[["spawn_bio"]] <- bioscale * recruit[["SpawnBio"]]
timeseries[["SpawnBio"]] <- bioscale * timeseries[["SpawnBio"]]
# x and y limits
if (is.null(ylim)) {
ylim <- c(0, 1.1 * max(recruit[["pred_recr"]], recruit[["exp_recr"]], recruit[["bias_adjusted"]]))
}
x <- recruit[["spawn_bio"]]
if (is.null(xlim)) {
xlim <- c(0, 1.1 * max(x))
}
# only add lines for environmentally dependent recruitment if it differs
# from expected recruitment without environmental link
show_env <- show_env & any(recruit[["with_env"]] != recruit[["exp_recr"]])
# store virgin and initial values
B0 <- sum(timeseries[["SpawnBio"]][timeseries[["Era"]] == "VIRG"], na.rm = TRUE)
B1 <- sum(timeseries[["SpawnBio"]][timeseries[["Era"]] == "INIT"], na.rm = TRUE)
R0 <- sum(timeseries[["Recruit_0"]][timeseries[["Era"]] == "VIRG"], na.rm = TRUE)
R1 <- sum(timeseries[["Recruit_0"]][timeseries[["Era"]] == "INIT"], na.rm = TRUE)
# work around for issue with Shepherd function producing 0 values in equilibrium
# use first non-zero value for each
if (B0 == 0) {
B0 <- head(recruit[["spawn_bio"]][recruit[["spawn_bio"]] != 0], 1)
}
if (R0 == 0) {
R0 <- head(recruit[["exp_recr"]][recruit[["exp_recr"]] != 0], 1)
}
if (B0 == B1 & R0 == R1) {
init <- FALSE
}
# scaling factor for axes (relative to B0/R0 or absolute)
if (relative) {
x.mult <- 1 / B0
y.mult <- 1 / R0
} else {
x.mult <- 1
y.mult <- 1
}
# color for points
if (ptcol[1] == "default") {
ptcol <- rev(rich.colors.short(nrow(recruit) + 10, alpha = 0.8))[-(1:10)]
color.caption <- paste(
" Point colors indicate year, with warmer",
"colors indicating earlier years and cooler colors",
"in showing later years."
)
} else {
color.caption <- ""
}
# prepare for legend
if (legend) {
legend_entries <- c(
expected, # expected
show_env, # with environmental link
bias_adjusted, # bias adjusted
estimated, # estimated points
virg, # virgin
init
) # initial equilibrium
legend_col <- colvec[c(3, 1, 2, 4, 3, 3)][legend_entries]
legend_bg <- c(NA, NA, NA, tail(ptcol)[1], NA, NA)[legend_entries]
legend_lwd <- c(2, 1, 1, NA, NA, NA)[legend_entries]
legend_lty <- ltyvec[c(3, 1, 2, 4, 3, 3)][legend_entries]
legend_pch <- c(NA, NA, NA, 21, 3, 4)[legend_entries]
legend_cex <- c(1, 1, 1, 1, 1.5, 1.5)[legend_entries]
legend_lab <- c(
"Exp. recruitment",
"Exp. recruitment with env. link",
"Exp. recruitment after bias adj.",
"Estimated recruitments",
"Unfished equilibrium",
"Initial equilibrium"
)[legend_entries]
}
StockRecruitCurve.fn <- function(text = FALSE) {
### a function to make the plots
if (!add) {
par(mar = c(4.5, 4.5, 1, 1))
# make empty plot (if not adding to existing plot)
plot(0,
type = "n", xlim = xlim * x.mult, ylim = ylim * y.mult,
xaxs = "i", yaxs = "i", xlab = xlab, ylab = ylab
)
}
if (show_env) {
# add line for expected recruitment with environmental variability
lines(x[order(x)] * x.mult, recruit[["with_env"]][order(x)] * y.mult,
lwd = 1, lty = ltyvec[1], col = colvec[1]
)
}
if (expected) {
# add line for expected recruitment
lines(x[order(x)] * x.mult, recruit[["exp_recr"]][order(x)] * y.mult,
lwd = 2, lty = ltyvec[3], col = colvec[3]
)
}
if (bias_adjusted) {
# add line for adjusted recruitment
lines(x * x.mult, recruit[["bias_adjusted"]] * y.mult,
lwd = 1, lty = ltyvec[2], col = colvec[2]
)
}
if (estimated) {
# add points for individual estimates
points(x * x.mult, recruit[["pred_recr"]] * y.mult, pch = 21, col = colvec[4], bg = ptcol)
}
if (text) {
# add text, but only label values with larger devs (in abs value)
show <- abs(recruit[["dev"]]) > textmindev
show[1] <- show[length(show)] <- TRUE # also include first & last years
text(x[show] * x.mult, recruit[["pred_recr"]][show] * y.mult,
labels = recruit[["Yr"]][show], pos = 2, cex = .7
)
}
# add point for virgin biomass/recruitment (if requested)
if (virg) {
points(B0 * x.mult, R0 * y.mult, pch = 3, cex = 1.5)
}
# add point for initial biomass/recruitment (if requested)
if (init) {
points(B1 * x.mult, R1 * y.mult, pch = 4, cex = 1.5)
}
# add legend
if (legend) {
# add sub-function:
legend.overlap <- function(x, y, ...) {
# function to figure out if a legend with inputs given by ...
# overlaps any of the points with coordinates x and y
# run legend without plotting
legend.out <- legend(..., plot = FALSE)
# get coordinates of legend boundaries
leg.left <- legend.out[["rect"]][["left"]]
leg.right <- legend.out[["rect"]][["left"]] + legend.out[["rect"]][["w"]]
leg.top <- legend.out[["rect"]][["top"]]
leg.bottom <- legend.out[["rect"]][["top"]] - legend.out[["rect"]][["h"]]
# test for overlap
if (any(x >= leg.left & x <= leg.right &
y >= leg.bottom & y <= leg.top)) {
return(TRUE)
} else {
return(FALSE)
}
}
# indicator if legend has been added successfully
legend_added <- FALSE
# if not ready to plot, look for best location
if (is.null(legendloc)) {
for (legendloc in c("topleft", "topright", "bottomright")) {
has_overlap <- legend.overlap(
x = x * x.mult,
y = recruit[["pred_recr"]] * y.mult,
legendloc, legend = legend_lab,
col = legend_col, pt.bg = legend_bg, lwd = legend_lwd,
lty = legend_lty, pch = legend_pch, bg = rgb(1, 1, 1, .6)
)
if (!has_overlap & !legend_added) {
legend(legendloc,
legend = legend_lab,
col = legend_col, pt.bg = legend_bg, lwd = legend_lwd,
lty = legend_lty, pch = legend_pch, bg = rgb(1, 1, 1, .6)
)
legend_added <- TRUE
}
} # end loop over legendloc values
if (!legend_added) {
warning(
"Legend in spawner-recruit plot overlaps at least 1 point\n",
"in the plot. Try running SSplotSpawnrecruit() with\n",
"adjustments to 'ylim' and/or 'legendloc' inputs."
)
legendloc <- "topleft"
}
}
# add legend at user-requested location or topleft with warning
if (!legend_added) {
legend(legendloc,
legend = legend_lab,
col = legend_col, pt.bg = legend_bg, lwd = legend_lwd,
lty = legend_lty, pch = legend_pch, bg = rgb(1, 1, 1, .5)
)
}
} # end commands to add legend
} # end StockRecruitCurve.fn
stock_vs_devs.fn <- function(text = FALSE) {
### a function to make the plots
if (!add) {
par(mar = c(4.5, 4.5, 1, 1))
# maximum spawning output relative to unfished equilibrium (usually 1)
xmax <- 1.05 * max(x / B0)
# make empty plot (if not adding to existing plot)
plot(0,
type = "n", xlim = c(0, xmax),
ylim = c(-1.1, 1.1) * max(abs(recruit[["dev"]]), na.rm = TRUE),
las = 1, xaxs = "i", yaxs = "i", xlab = labels[4], ylab = labels[6]
)
}
abline(h = 0, col = "grey")
points(x / B0, recruit[["dev"]], pch = 21, bg = ptcol, col = colvec[4], cex = 1.5)
if (text) {
# add text, but only label values with larger devs (in abs value)
show <- abs(recruit[["dev"]]) > textmindev
show[1] <- show[length(show)] <- TRUE # also include first & last years
text(x[show] / B0, recruit[["dev"]][show],
labels = recruit[["Yr"]][show], pos = 2, cex = .7
)
}
# add point for virgin biomass/recruitment (if requested)
if (virg) {
points(1, 0, pch = 3, cex = 2, lwd = 2)
}
# add point for initial biomass/recruitment (if requested)
if (init) {
points(B1 / B0, 0, pch = 4, cex = 2)
}
## # add legend
## if(legend){
## legend(legendloc, legend=legend_lab, col=legend_col, lwd=legend_lwd,
## pch=legend_pch, bg=rgb(1,1,1,.9))
## }
}
if (plot) {
if (1 %in% subplot) {
StockRecruitCurve.fn()
}
if (2 %in% subplot) {
StockRecruitCurve.fn(text = TRUE)
}
if (3 %in% subplot) {
stock_vs_devs.fn(text = TRUE)
}
}
if (print) {
if (1 %in% subplot) {
file <- "SR_curve.png"
caption <- paste("Stock-recruit curve.", color.caption)
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
StockRecruitCurve.fn()
dev.off()
}
if (2 %in% subplot) {
file <- "SR_curve2.png"
caption <- paste0(
"Stock-recruit curve with labels on first, last, and ",
"years with (log) deviations > ", textmindev, ".",
color.caption
)
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
StockRecruitCurve.fn(text = TRUE)
dev.off()
}
if (3 %in% subplot) {
file <- "SR_resids.png"
caption <- paste0(
"Deviations around the stock-recruit curve. ",
"Labels are on first, last, and ",
"years with (log) deviations > ", textmindev, ".",
color.caption
)
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
stock_vs_devs.fn(text = TRUE)
dev.off()
}
}
if (!is.null(plotinfo)) plotinfo[["category"]] <- "S-R"
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.