#' Plot numbers-at-age related data and fits.
#'
#' Plot numbers-at-age related data and fits from Stock Synthesis output.
#' Plots include bubble plots, mean age, equilibrium age composition,
#' sex-ratio, and ageing imprecision patterns.
#'
#'
#' @template replist
#' @param subplots vector controlling which subplots to create
#' Numbering of subplots is as follows,
#' \itemize{
#' \item 1: Expected numbers at age
#' \item 2: Mean age in the population
#' \item 3: Fraction female in numbers at age
#' \item 4: Equilibrium age distribution
#' \item 5: Ageing imprecision: SD of observed age (plot using image() formerly
#' included in this group but now replaced by better distribution plots)
#' \item 6: Expected numbers at length
#' \item 7: Mean length in the population
#' \item 8: Fraction female in numbers at length
#' \item 9: no plot yet
#' \item 10: Distribution of observed age at true age by ageing error type
#' }
#' @template plot
#' @template print
#' @param numbers.unit Units for numbers. Default (based on typical Stock Synthesis
#' setup) is thousands (numbers.unit=1000).
#' @param areas optional subset of areas to plot for spatial models
#' @param areanames names for areas. Default is to use Area1, Area2,...
#' @template areacols
#' @param pntscalar maximum bubble size for bubble plots; each plot scaled
#' independently based on this maximum size and the values plotted. Often some
#' plots look better with one value and others with a larger or smaller value.
#' Default=2.6
#' @param bub.bg background color for bubbles
#' (no control over black border at this time)
#' @param bublegend Add legend with example bubble sizes?
#' @param period indicator of whether to make plots using numbers at age just
#' from the beginning ("B") or middle of the year ("M") (new option starting
#' with SSv3.11)
#' @param meanlines add lines for mean age or length on top of bubble plots
#' @param add add to existing plot? (not yet implemented)
#' @template labels
#' @template pwidth
#' @template pheight
#' @template punits
#' @template ptsize
#' @template res
#' @template cex.main
#' @template plotdir
#' @template mainTitle
#' @template verbose
#' @author Ian Stewart, Ian Taylor
#' @export
#' @seealso [SS_output()], [SS_plots()]
SSplotNumbers <-
function(replist, subplots = c(1:10),
plot = TRUE, print = FALSE,
numbers.unit = 1000,
areas = "all",
areanames = "default",
areacols = NULL,
pntscalar = 2.6,
bub.bg = gray(0.5, alpha = 0.5),
bublegend = TRUE,
period = c("B", "M"),
meanlines = TRUE,
add = FALSE,
labels = c(
"Year", # 1
"Age", # 2
"True age (yr)", # 3
"SD of observed age (yr)", # 4
"Mean observed age (yr)", # 5
"Mean age (yr)", # 6
"mean age in the population", # 7
"Ageing imprecision", # 8
"Numbers at age at equilibrium", # 9
"Equilibrium age distribution", # 10
"Fraction female in numbers at age", # 11
"Length", # 12
"Mean length (cm)", # 13
"mean length (cm) in the population", # 14
"expected numbers at age", # 15
"Beginning of year", # 16
"Middle of year", # 17
"expected numbers at length", # 18
# 19 below is out of order, runumbering others would be tedious
"Fraction female in numbers at length"
), # 19
pwidth = 6.5, pheight = 6.5, punits = "in", res = 300, ptsize = 10,
cex.main = 1,
plotdir = "default",
mainTitle = FALSE,
verbose = TRUE) {
# table to store information on each plot
plotinfo <- NULL
natage <- replist[["natage"]]
if (is.null(natage)) {
message(
"Skipped some plots because NUMBERS_AT_AGE unavailable\n",
" in the report file. The starter file may be set to produce\n",
" limited report detail."
)
} else {
# get stuff from replist
ngpatterns <- max(natage[["Bio_Pattern"]])
natlen <- replist[["natlen"]]
nsexes <- replist[["nsexes"]]
nareas <- replist[["nareas"]]
nseasons <- replist[["nseasons"]]
spawnseas <- replist[["spawnseas"]]
morph_indexing <- replist[["morph_indexing"]]
accuage <- replist[["accuage"]]
agebins <- replist[["agebins"]]
endyr <- replist[["endyr"]]
N_ageerror_defs <- replist[["N_ageerror_defs"]]
AAK <- replist[["AAK"]]
age_error_mean <- replist[["age_error_mean"]]
age_error_sd <- replist[["age_error_sd"]]
lbinspop <- replist[["lbinspop"]]
nlbinspop <- replist[["nlbinspop"]]
mainmorphs <- replist[["mainmorphs"]]
SS_versionNumeric <- replist[["SS_versionNumeric"]]
if (plotdir == "default") {
plotdir <- replist[["inputs"]][["dir"]]
}
if (areas[1] == "all") {
areas <- 1:nareas
} else {
if (length(intersect(areas, 1:nareas)) != length(areas)) {
stop("Input 'areas' should be 'all' or a vector of values between 1 and nareas.")
}
}
if (areanames[1] == "default") {
areanames <- paste0("area", 1:nareas)
}
# set default area-specific colors if not specified
areacols <- get_areacols(areacols, nareas)
if (SS_versionNumeric <= 3.1) {
stop("numbers at age plots no longer supported for SS version 3.10 and earlier")
}
###########
# Numbers at age plots
# combining across platoons/submorphs and growth patterns
column1 <- 12
if (SS_versionNumeric >= 3.30) {
column1 <- 13
}
remove <- -(1:(column1 - 1)) # removes first group of columns
if ("Settlement" %in% names(natage)) {
settlement <- unique(natage[["Settlement"]])
if (length(settlement) > 1) {
message("Numbers at age plots only show first settlement event")
}
}
birth_seas_name <- grep("^BirthSeas", colnames(natage), value = TRUE)
bseas <- unique(natage[[birth_seas_name]])
if (length(bseas) > 1) {
message("Numbers at age plots are for only the first birth season")
}
if (ngpatterns > 1) {
message(
"Numbers at age plots may not deal correctly with growth patterns:",
"no guarantees!"
)
}
if (nseasons > 1) {
message("Numbers at age plots are for season 1 only")
# change plot labels for seasonal models
labels[16] <- gsub(pattern = "of year", replacement = "of season 1", x = labels[16])
labels[17] <- gsub(pattern = "of year", replacement = "of season 1", x = labels[17])
}
for (iarea in areas) {
for (iperiod in seq_along(period)) {
for (m in 1:nsexes) {
# warning! implementation of birthseasons may not be correct in this section
# data frame to combine values across factors
natagetemp_all <- natage[natage[["Area"]] == iarea &
natage[["Sex"]] == m &
natage[["Seas"]] == 1 &
natage[["Era"]] != "VIRG" &
!is.na(natage$"0") &
natage[["Yr"]] < (endyr + 2) &
natage[[birth_seas_name]] == min(bseas), ]
# natage[["Bio_Pattern"]]==1,] # formerly filtered
natagetemp_all <- natagetemp_all[natagetemp_all$"Beg/Mid" == period[iperiod], ]
# create data frame with 0 values to fill across platoons/submorphs
morphlist <- unique(natagetemp_all[["Platoon"]])
natagetemp0 <- natagetemp_all[natagetemp_all[["Platoon"]] == morphlist[1] &
natagetemp_all[["Bio_Pattern"]] == 1, ]
for (iage in 0:accuage) {
# matrix of zeros for upcoming calculations
natagetemp0[, column1 + iage] <- 0
}
for (imorph in seq_along(morphlist)) {
for (igp in 1:ngpatterns) {
natagetemp_imorph_igp <-
natagetemp_all[natagetemp_all[["Platoon"]] == morphlist[imorph] &
natagetemp_all[["Bio_Pattern"]] == igp, ]
natagetemp0[, column1 + 0:accuage] <-
natagetemp0[, column1 + 0:accuage] +
natagetemp_imorph_igp[, column1 + 0:accuage]
} # end growth pattern loop
} # end morph loop
if (ngpatterns > 0) {
natagetemp0[["Bio_Pattern"]] == 999
}
nyrsplot <- nrow(natagetemp0)
resx <- rep(natagetemp0[["Yr"]], accuage + 1)
resy <- NULL
for (i in 0:accuage) {
resy <- c(resy, rep(i, nyrsplot))
}
resz <- NULL
for (i in column1 + 0:accuage) {
# numbers here are scaled by units
resz <- c(resz, numbers.unit * natagetemp0[, i])
}
# clean up big numbers
units <- ""
if (max(resz, na.rm = TRUE) > 1e9) {
resz <- resz / 1e9
units <- "billion"
}
if (max(resz, na.rm = TRUE) > 1e6 & units == "") {
resz <- resz / 1e6
units <- "million"
}
if (max(resz, na.rm = TRUE) > 1e3 & units == "") {
resz <- resz / 1e3
units <- "thousand"
}
# assign unique name to data frame for area, sex (beginning of year only)
if (iperiod == 1) {
assign(paste0("natagetemp0area", iarea, "sex", m), natagetemp0)
}
if (m == 1 & nsexes == 1) sextitle <- ""
if (m == 1 & nsexes == 2) sextitle <- " of females"
if (m == 2) sextitle <- " of males"
if (nareas > 1) sextitle <- paste0(sextitle, " in ", areanames[iarea])
if (!period[iperiod] %in% c("B", "M")) {
stop("'period' input to SSplotNumbers should include only 'B' or 'M'")
}
if (period[iperiod] == "B") {
periodtitle <- labels[16]
fileperiod <- "_beg"
}
if (period[iperiod] == "M") {
periodtitle <- labels[17]
fileperiod <- "_mid"
}
# title for use in title or caption
plottitle1 <- paste0(
periodtitle, " ", labels[15], sextitle,
" in (max ~ ", format(round(max(resz, na.rm = TRUE), 1), nsmall = 1),
" ", units, ")"
)
# title if requested
main <- ""
if (mainTitle) {
main <- plottitle1
}
### calculations related to mean age
# removing the first columns to get just numbers
natagetemp1 <- as.matrix(natagetemp0[, remove])
ages <- 0:accuage
natagetemp2 <- as.data.frame(natagetemp1)
natagetemp2[["sum"]] <- as.vector(apply(natagetemp1, 1, sum))
# remove rows with 0 fish (i.e. no growth pattern in this area)
natagetemp0 <- natagetemp0[natagetemp2[["sum"]] > 0, ]
natagetemp1 <- natagetemp1[natagetemp2[["sum"]] > 0, ]
natagetemp2 <- natagetemp2[natagetemp2[["sum"]] > 0, ]
prodmat <- t(natagetemp1) * ages
prodsum <- as.vector(apply(prodmat, 2, sum))
natagetemp2[["sumprod"]] <- prodsum
natagetemp2[["meanage"]] <-
natagetemp2[["sumprod"]] / natagetemp2[["sum"]] - (natagetemp0[[birth_seas_name]] - 1) / nseasons
natageyrs <- sort(unique(natagetemp0[["Yr"]]))
if (iperiod == 1) {
natageyrsB <- natageyrs # unique name for beginning of year
}
meanage <- 0 * natageyrs
for (i in seq_along(natageyrs)) {
# averaging over values within a year (depending on birth season)
meanage[i] <- sum(natagetemp2[["meanage"]][natagetemp0[["Yr"]] == natageyrs[i]] *
natagetemp2[["sum"]][natagetemp0[["Yr"]] == natageyrs[i]]) /
sum(natagetemp2[["sum"]][natagetemp0[["Yr"]] == natageyrs[i]])
}
if (m == 1 & nsexes == 2) {
meanagef <- meanage # save value for females in 2 sex models
}
ylab <- labels[6]
main <- ""
plottitle2 <- paste(periodtitle, labels[7])
if (nareas > 1) plottitle2 <- paste(plottitle2, "in", areanames[iarea])
if (mainTitle) {
main <- plottitle2
}
ageBubble.fn <- function() {
# bubble plot with line
bubble3(
x = resx, y = resy, z = resz,
xlab = labels[1], ylab = labels[2],
legend = bublegend, bg.open = bub.bg,
main = main, maxsize = (pntscalar + 1.0),
las = 1, cex.main = cex.main, allopen = TRUE
)
# add line for mean age
if (meanlines & all(!is.nan(meanage))) {
lines(natageyrs, meanage, col = "red", lwd = 3)
}
}
meanAge.fn <- function() {
# mean age for males and females
ylim <- c(0, max(meanage, meanagef, na.rm = TRUE))
plot(natageyrs, meanage,
col = "blue", lty = 1, pch = 4, xlab = labels[1],
ylim = ylim, type = "o", ylab = ylab, main = main, cex.main = cex.main
)
points(natageyrs, meanagef, col = "red", lty = 2, pch = 1, type = "o")
legend("bottomleft",
bty = "n", c("Females", "Males"), lty = c(2, 1),
pch = c(1, 4), col = c("red", "blue")
)
}
if (plot) {
if (1 %in% subplots) ageBubble.fn()
if (2 %in% subplots & m == 2 & nsexes == 2) {
meanAge.fn()
}
}
if (print) {
filepartsex <- paste0("_sex", m)
filepartarea <- ""
if (nareas > 1) {
filepartarea <- paste0("_", areanames[iarea])
}
if (1 %in% subplots) {
file <- paste0(
"numbers1", filepartarea, filepartsex,
fileperiod, ".png"
)
caption <- plottitle1
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
ageBubble.fn()
dev.off()
}
# make 2-sex plot after looping over both sexes
if (2 %in% subplots & m == 2 & nsexes == 2) {
file <- paste0("numbers2_meanage", filepartarea,
fileperiod, ".png",
sep = ""
)
caption <- plottitle2
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
meanAge.fn()
dev.off()
}
} # end printing to PNG file
} # end sex loop
} # end period loop
} # end area loop
if (nsexes > 1) {
for (iarea in areas) {
plottitle3 <- paste0(labels[11], sep = "")
if (nareas > 1) {
plottitle3 <- paste0(plottitle3, " for ", areanames[iarea])
}
main <- ""
if (mainTitle) {
main <- plottitle3
}
# get objects that were assigned earlier
natagef <- get(paste0("natagetemp0area", iarea, "sex", 1))
natagem <- get(paste0("natagetemp0area", iarea, "sex", 2))
natagefyrs <- natagef[["Yr"]]
natageratio <- as.matrix(natagef[, remove] /
(natagem[, remove] + natagef[, remove]))
natageratio[is.nan(natageratio)] <- NA
if (diff(range(natageratio, finite = TRUE)) != 0) {
numbersRatioAge.fn <- function(...) {
contour(natagefyrs, 0:accuage, natageratio,
xaxs = "i", yaxs = "i",
xlab = labels[1], ylab = labels[2], main = main,
cex.main = cex.main, ...
)
}
if (plot & 3 %in% subplots) {
numbersRatioAge.fn(labcex = 1)
}
if (print & 3 %in% subplots) {
filepart <- ""
if (nareas > 1) filepart <- paste0("_", areanames[iarea], filepart)
file <- paste0("numbers3_frac_female_age", filepart, ".png")
caption <- plottitle3
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
numbersRatioAge.fn(labcex = 0.6)
dev.off()
}
} else {
if (verbose) {
message(
"skipped sex ratio contour plot because females=males ",
"for all ages and years"
)
}
}
} # end area loop
} # end if nsexes>1
##########
# repeat code above for numbers at length (subplots 6, 7 and 8)
if (length(intersect(6:8, subplots)) > 0 & !is.null(natlen) & all(!is.na(lbinspop))) {
column1 <- column1 - 1 # because index of lengths starts at 1, not 0 as in ages
for (iarea in areas) {
for (iperiod in seq_along(period)) {
for (m in 1:nsexes) {
# warning! implementation of birthseasons may not be correct in this section
# data frame to combine values across factors
natlentemp_all <- natlen[natlen[["Area"]] == iarea &
natlen[["Sex"]] == m &
natlen[["Seas"]] == 1 &
natlen[["Era"]] != "VIRG" &
# natlen[["Era"]]!="FORE" &
natlen[["Yr"]] < (endyr + 2) &
natlen[[birth_seas_name]] == min(bseas), ]
# natlen[["Bio_Pattern"]]==1,] # formerly filtered
natlentemp_all <- natlentemp_all[natlentemp_all$"Beg/Mid" == period[iperiod], ]
# create data frame with 0 values to fill across platoons/submorphs
morphlist <- unique(natlentemp_all[["Platoon"]])
natlentemp0 <- natlentemp_all[natlentemp_all[["Platoon"]] == morphlist[1] &
natlentemp_all[["Bio_Pattern"]] == 1, ]
for (ilen in 1:nlbinspop) {
# matrix of zeros for upcoming calculations
natlentemp0[, column1 + ilen] <- 0
}
for (imorph in seq_along(morphlist)) {
for (igp in 1:ngpatterns) {
natlentemp_imorph_igp <-
natlentemp_all[natlentemp_all[["Platoon"]] == morphlist[imorph] &
natlentemp_all[["Bio_Pattern"]] == igp, ]
natlentemp0[, column1 + 1:nlbinspop] <-
natlentemp0[, column1 + 1:nlbinspop] +
natlentemp_imorph_igp[, column1 + 1:nlbinspop]
} # end growth pattern loop
} # end morph loop
if (ngpatterns > 0) natlentemp0[["Bio_Pattern"]] == 999
nyrsplot <- nrow(natlentemp0)
resx <- rep(natlentemp0[["Yr"]], nlbinspop)
resy <- NULL
for (ilen in 1:nlbinspop) resy <- c(resy, rep(lbinspop[ilen], nyrsplot))
resz <- NULL
for (ilen in column1 + 1:nlbinspop) {
# numbers here are scaled by units
resz <- c(resz, numbers.unit * natlentemp0[, ilen])
}
# clean up big numbers
units <- ""
if (max(resz, na.rm = TRUE) > 1e9) {
resz <- resz / 1e9
units <- "billion"
}
if (max(resz, na.rm = TRUE) > 1e6 & units == "") {
resz <- resz / 1e6
units <- "million"
}
if (max(resz, na.rm = TRUE) > 1e3 & units == "") {
resz <- resz / 1e3
units <- "thousand"
}
# assign unique name to data frame for area, sex
assign(paste0("natlentemp0area", iarea, "sex", m), natlentemp0)
if (m == 1 & nsexes == 1) sextitle <- ""
if (m == 1 & nsexes == 2) sextitle <- " of females"
if (m == 2) sextitle <- " of males"
if (nareas > 1) sextitle <- paste0(sextitle, " in ", areanames[iarea])
if (period[iperiod] == "B") {
periodtitle <- labels[16]
} else if (period[iperiod] == "M") {
periodtitle <- labels[17]
} else {
stop("'period' input to SSplotNumbers should include only 'B' or 'M'")
}
plottitle1 <- paste0(
periodtitle, " ", labels[18], sextitle,
" in (max ~ ", format(round(max(resz, na.rm = TRUE), 1), nsmall = 1),
" ", units, ")"
)
main <- ""
if (mainTitle) {
main <- plottitle1
}
### calculations related to mean len
# removing the first columns to get just numbers
natlentemp1 <- as.matrix(natlentemp0[, remove])
natlentemp2 <- as.data.frame(natlentemp1)
natlentemp2[["sum"]] <- as.vector(apply(natlentemp1, 1, sum))
# remove rows with 0 fish (i.e. no growth pattern in this area)
natlentemp0 <- natlentemp0[natlentemp2[["sum"]] > 0, ]
natlentemp1 <- natlentemp1[natlentemp2[["sum"]] > 0, ]
natlentemp2 <- natlentemp2[natlentemp2[["sum"]] > 0, ]
prodmat <- t(natlentemp1) * lbinspop
prodsum <- as.vector(apply(prodmat, 2, sum))
natlentemp2[["sumprod"]] <- prodsum
natlentemp2[["meanlen"]] <-
natlentemp2[["sumprod"]] / natlentemp2[["sum"]] - (natlentemp0[[birth_seas_name]] - 1) / nseasons
natlenyrs <- sort(unique(natlentemp0[["Yr"]]))
if (iperiod == 1) {
natlenyrsB <- natlenyrs # unique name for years associated with period=="B"
}
meanlen <- 0 * natlenyrs
for (i in seq_along(natlenyrs)) {
# averaging over values within a year (depending on birth season)
meanlen[i] <-
sum(natlentemp2[["meanlen"]][natlentemp0[["Yr"]] == natlenyrs[i]] *
natlentemp2[["sum"]][natlentemp0[["Yr"]] == natlenyrs[i]]) /
sum(natlentemp2[["sum"]][natlentemp0[["Yr"]] == natlenyrs[i]])
}
if (m == 1 & nsexes == 2) {
meanlenf <- meanlenf <- meanlen # save value for females in 2 sex models
}
ylab <- labels[13]
plottitle2 <- paste(periodtitle, labels[14])
if (nareas > 1) plottitle2 <- paste(plottitle2, "in", areanames[iarea])
main <- ""
if (mainTitle) {
main <- plottitle2
}
lenBubble.fn <- function() {
# bubble plot with line
bubble3(
x = resx, y = resy, z = resz,
xlab = labels[1], ylab = labels[12],
legend = bublegend, bg.open = bub.bg,
main = main, maxsize = (pntscalar + 1.0),
las = 1, cex.main = cex.main, allopen = TRUE
)
# add line for mean length
if (meanlines & all(!is.nan(meanlen))) {
lines(natlenyrs, meanlen, col = "red", lwd = 3)
}
}
meanLen.fn <- function() {
# mean length for males and females
ylim <- c(0, max(meanlen, meanlenf))
plot(natlenyrs, meanlen,
col = "blue", lty = 1, pch = 4, xlab = labels[1], ylim = ylim,
type = "o", ylab = ylab, main = main, cex.main = cex.main
)
points(natlenyrs, meanlenf, col = "red", lty = 2, pch = 1, type = "o")
legend("bottomleft",
bty = "n", c("Females", "Males"),
lty = c(2, 1), pch = c(1, 4), col = c("red", "blue")
)
}
if (plot) {
if (6 %in% subplots) lenBubble.fn()
if (7 %in% subplots & m == 2 & nsexes == 2) meanLen.fn()
}
if (print) {
filepartsex <- paste0("_sex", m)
filepartarea <- ""
if (nareas > 1) {
filepartarea <- paste0("_", areanames[iarea])
}
if (6 %in% subplots) {
file <- paste0(
"numbers6_len", filepartarea,
filepartsex, ".png"
)
caption <- plottitle1
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
lenBubble.fn()
dev.off()
}
# make 2-sex plot after looping over both sexes
if (7 %in% subplots & m == 2 & nsexes == 2) {
file <- paste0("numbers7_meanlen", filepartarea, ".png")
caption <- plottitle2
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
meanLen.fn()
dev.off()
}
} # end printing of plot 14
} # end sex loop
} # end period loop
} # end area loop
if (nsexes > 1) {
for (iarea in areas) {
# get objects that were assigned earlier
natlenf <- get(paste0("natlentemp0area", iarea, "sex", 1))
natlenm <- get(paste0("natlentemp0area", iarea, "sex", 2))
natlenratio <- as.matrix(natlenf[, remove] /
(natlenm[, remove] + natlenf[, remove]))
if (diff(range(natlenratio, finite = TRUE)) != 0) {
main <- ""
caption <- labels[19]
if (nareas > 1) {
caption <- paste0(main, " for ", areanames[iarea])
}
if (mainTitle) {
main <- caption
}
numbersRatioLen.fn <- function(...) {
z <- natlenratio
# check for mismatch that caused crash in one particular model
# note from Ian (11/1/2018): taking too long to sort this out so
# just skipping the plot for the rare case
# that has the error
if (length(natlenyrsB) == nrow(z)) {
contour(natlenyrsB, lbinspop, z,
xaxs = "i", yaxs = "i", xlab = labels[1], ylab = labels[12],
main = main, cex.main = cex.main, ...
)
} else {
warning(
"Skipping plot of sex ratio by length and year\n ",
"due to mismatch in length of table and vector of years.\n ",
"This may be due to 0 values in the table."
)
}
}
if (plot & 8 %in% subplots) {
numbersRatioLen.fn(labcex = 1)
}
if (print & 8 %in% subplots) {
filepart <- ""
if (nareas > 1) {
filepart <- paste0("_", areanames[iarea], filepart)
}
file <- paste0("numbers8_frac_female_len", filepart, ".png")
caption <- caption
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
numbersRatioLen.fn(labcex = 0.6)
dev.off()
}
} else {
if (verbose) {
message(
"skipped sex ratio contour plot because females=males",
" for all lengths and years"
)
}
}
} # end area loop
} # end if nsexes>1
} # end numbers at length plots
##########
# plot of equilibrium age composition by sex and area
equilibfun <- function() {
# subset to unfished equilibrium
equilage <- natage[natage[["Era"]] == "VIRG", ]
# remove rows with all zeros
equilage <- equilage[as.vector(apply(equilage[, remove], 1, sum)) > 0, ]
# figure out birth season / spawning seaso
BirthSeas <- spawnseas
if (!(spawnseas %in% bseas)) BirthSeas <- min(bseas)
if (length(bseas) > 1) {
message("showing equilibrium age for first birth season", BirthSeas)
}
# sort out plot title
main <- ""
if (mainTitle) {
main <- labels[10]
}
ymax <- max(equilage[
equilage[[birth_seas_name]] == BirthSeas &
equilage[["Seas"]] == BirthSeas &
equilage[["Beg/Mid"]] == "B",
"0"
])
plot(0,
type = "n", xlim = c(0, accuage),
ylim = c(0, 1.05 * ymax),
xaxs = "i", yaxs = "i", xlab = "Age", ylab = labels[9],
main = main, cex.main = cex.main
)
# now fill in legend
legendlty <- NULL
legendcol <- NULL
legendlegend <- NULL
for (iarea in areas) {
for (m in 1:nsexes) {
# subset to beginning of the year and birth season within each area and sex
equilagetemp <- equilage[equilage[["Area"]] == iarea &
equilage[["Sex"]] == m &
equilage[[birth_seas_name]] == BirthSeas &
equilage[["Seas"]] == BirthSeas &
equilage[["Beg/Mid"]] == "B", ]
# subset to columns with numbers at age values only
equilagetemp <- equilagetemp[, remove]
if (nrow(equilagetemp) == 1) {
lines(0:accuage, equilagetemp, lty = m, lwd = 3, col = areacols[iarea])
} else {
matplot(
x = 0:accuage, y = t(equilagetemp), type = "l",
lty = m, lwd = 3, col = areacols[iarea], add = TRUE
)
message(
"Multiple matching lines in equilibrium plot indicate multiple ",
"morphs or platoons within each area/sex combination."
)
}
legendlty <- c(legendlty, m)
legendcol <- c(legendcol, areacols[iarea])
if (m == 1 & nsexes == 1) sextitle <- ""
if (m == 1 & nsexes == 2) sextitle <- "Females"
if (m == 2) sextitle <- "Males"
if (nareas > 1) {
sextitle <- paste0(sextitle, " in ", areanames[iarea])
}
legendlegend <- c(legendlegend, sextitle)
}
}
if (length(legendlegend) > 1) {
legend("topright", legend = legendlegend, col = legendcol, lty = legendlty, lwd = 3)
}
}
if (plot & 4 %in% subplots) {
equilibfun()
}
if (print & 4 %in% subplots) {
file <- "numbers4_equilagecomp.png"
caption <- labels[10]
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
equilibfun()
dev.off()
} # end print to PNG
# plot the ageing imprecision for all age methods
if (!is.null(N_ageerror_defs) && N_ageerror_defs > 0) {
xvals <- age_error_sd[["age"]] + 0.5
yvals <- age_error_sd[, -1]
ylim <- c(0, max(yvals))
if (N_ageerror_defs == 1) {
colvec <- "black"
} else {
colvec <- rich.colors.short(N_ageerror_defs)
}
main <- ""
if (mainTitle) {
main <- labels[8]
}
ageingfun <- function() {
matplot(xvals, yvals,
ylim = ylim, type = "o", pch = 1, lty = 1, col = colvec,
xlab = labels[3], ylab = labels[4], main = main, cex.main = cex.main
)
abline(h = 0, col = "grey") # grey line at 0
legend("topleft",
bty = "n", pch = 1, lty = 1, col = colvec,
# more columns for crazy models like hake with many ageing matrices
ncol = ifelse(N_ageerror_defs < 20, 1, 2),
legend = paste("Ageing method", 1:N_ageerror_defs)
)
}
# check for bias in ageing error pattern
ageingbias <- age_error_mean[, -1] - (age_error_mean[["age"]] + 0.5)
if (mean(ageingbias == 0) != 1) {
ageingfun2 <- function() {
yvals <- age_error_mean[, -1]
ylim <- c(0, max(yvals))
matplot(xvals, yvals,
ylim = ylim, type = "o", pch = 1, lty = 1, col = colvec,
xlab = labels[3], ylab = labels[5], main = main
)
abline(h = 0, col = "grey") # grey line at 0
abline(0, 1, col = "grey") # grey line with slope = 1
legend("topleft",
bty = "n", pch = 1, lty = 1, col = colvec,
# more columns for crazy models like hake with many ageing matrices
ncol = ifelse(N_ageerror_defs < 20, 1, 2),
legend = paste("Ageing method", 1:N_ageerror_defs)
)
}
}
# run functions to make requested plots
if (plot & 5 %in% subplots) {
# make plots of age error standard deviations
ageingfun()
# make plots of age error means
if (mean(ageingbias == 0) != 1) ageingfun2()
}
if (print & 5 %in% subplots) {
# make files with plots of age error standard deviations
file <- "numbers5_ageerrorSD.png"
caption <- paste0(labels[8], ": ", labels[4])
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
ageingfun()
dev.off()
# make files with plots of age error means
if (mean(ageingbias == 0) != 1) {
file <- "numbers5_ageerrorMeans.png"
caption <- paste0(labels[8], ": ", labels[5])
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
ageingfun2()
dev.off()
}
} # end print to PNG
if (10 %in% subplots) {
# newer plot of ageing matrix as histogram-style figure
plotinfo.tmp <- SSplotAgeMatrix(
replist = replist, option = 2,
plot = plot, print = print,
plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits,
res = res, ptsize = ptsize,
cex.main = cex.main,
mainTitle = mainTitle
)
plotinfo <- rbind(plotinfo, plotinfo.tmp)
}
} # end if AAK
} # end if data available
if (!is.null(plotinfo)) plotinfo[["category"]] <- "Numbers"
return(invisible(plotinfo))
} # end function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.