-
Notifications
You must be signed in to change notification settings - Fork 35
/
SS_plots.R
1715 lines (1654 loc) · 71.4 KB
/
SS_plots.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#' plot many quantities related to output from Stock Synthesis
#'
#' Creates a user-chosen set of plots, including biological quantities, time
#' series, and fits to data. Plots are sent to R GUI, single PDF file, or
#' multiple PNG files. This is now just a wrapper which calls on separate
#' functions to make all the plots.
#'
#' @template replist
#' @param plot Plot sets to be created, see list of plots below. Use to
#' specify only those plot sets of interest, e.g., c(1,2,5,10). Plots for data
#' not available in the model run will automatically be skipped, whether called
#' or not.
#' Current grouping of plots is as follows:
#' \enumerate{
#' \item Biology
#' \item Selectivity and retention
#' \item Timeseries
#' \item Recruitment deviations
#' \item Recruitment bias adjustment
#' \item Spawner-recruit
#' \item Catch
#' \item SPR
#' \item Discards
#' \item Mean weight
#' \item Indices
#' \item Numbers at age
#' \item Length comp data (and generalized size comp data)
#' \item Age comp data
#' \item Conditional age-at-length data
#' \item Length comp fits (and generalized size comp fits)
#' \item Age comp fits
#' \item Conditional age-at-length fits
#' \item Francis and Punt conditional age-at-length comp fits
#' \item Mean length-at-age and mean weight-at-age
#' \item Tags
#' \item Yield
#' \item Movement
#' \item Data range
#' \item Parameter distributions
#' \item Diagnostic tables
#' }
#'
#' @param pdf Send plots to PDF file instead of R GUI?
#' @param png Send plots to PNG files instead of R GUI?
#' @param html Run [SS_html()] on completion? By default has same
#' value as `png`.
#' @param printfolder The sub-directory under 'dir' (see below) in which the
#' PNG files will be located. The default sub-directory is "plots".
#' The directory will be created if it doesn\'t exist.
#' If 'printfolder' is set to "", it is ignored and the PNG files will be located
#' in the directory specified by 'dir'.
#' @param dir The directory in which a PDF file (if requested) will be created
#' and within which the printfolder sub-directory (see above) will be created
#' if png=TRUE. By default it will be the same directory that the report file
#' was read from by the `SS_output` function. Alternatives to the default
#' can be either relative (to the working directory) or absolute paths.
#' The function will attempt to create the directory it doesn't exist, but it
#' does not do so recursively.
#' @template fleets
#' @param areas Either the string "all", or a vector of numerical values, like
#' c(1,3), listing areas for which plots should be made in a multi-area model.
#' By default, plots will be made for all areas (excepting cases where the
#' function has not yet been updated for multi-area models). Default="all".
#' @template fleetnames
#' @param fleetcols Either the string "default", or a vector of colors to use
#' for each fleet. Default="default".
#' @param fleetlty Vector of line types used for each fleet in some plots.
#' Default=1.
#' @param fleetpch Vector of point types used for each fleet in some plots.
#' Default=1.
#' @template lwd
#' @template areacols
#' @param areanames Optional vector of names for each area used in titles.
#' Default="default".
#' @template verbose
#' @param uncertainty Include values in plots showing estimates of uncertainty
#' (requires positive definite hessian in model? Default=TRUE.
#' @param forecastplot Include forecast years in the timeseries plots and
#' plots of time-varying quantities?
#' @param datplot Plot the data by itself? This is useful in document
#' preparation, but doesn't change across alternative model runs with the same
#' data, so can be committed to save time once the plots have been created once.
#' Setting datplot=FALSE is equivalent to leaving off plots 15 and 16.
#' Default=TRUE.
#' @param Natageplot Plot the expected numbers at age bubble plots and mean-age
#' time series? Default=TRUE.
#' @param samplesizeplots Show sample size plots? Default=TRUE.
#' @param compresidplots Show residuals for composition plots?
#' @param comp.yupper Upper limit on ymax for polygon/histogram composition
#' plots. This avoids scaling all plots to have max=1 if there is a vector
#' with only a single observed fish in it. Default=0.4.
#' @param sprtarg Specify the F/SPR proxy target. Default=0.4.
#' @param btarg Target %unfished to be used in plots showing %unfished. May be
#' omitted by setting to NA.
#' @param minbthresh Threshold depletion to be used in plots showing depletion.
#' May be omitted by setting to NA.
#' @param pntscalar This scalar defines the maximum bubble size for bubble
#' plots. This option is still available but a better choice is to use
#' bub.scale.pearson and bub.scale.dat, which are allow the same scaling
#' throughout all plots.
#' @param pntscalar.nums This scalar defines the maximum bubble size for
#' numbers-at-age and numbers-at-length plots.
#' @param pntscalar.tags This scalar defines the maximum bubble size for
#' tagging plots.
#' @param bub.scale.pearson Character expansion (cex) value for a proportion of
#' 1.0 in bubble plot of Pearson residuals. Default=1.5.
#' @param bub.scale.dat Character expansion (cex) value for a proportion of 1.0
#' in bubble plot of composition data. Default=3.
#' @param minnbubble This defines the minimum number of years below which blank
#' years will be added to bubble plots to avoid cropping. Default=8.
#' @param aalyear Years to plot multi-panel conditional age-at-length fits for
#' all length bins; must be in a "c(YYYY,YYYY)" format. Useful for checking the
#' fit of a dominant year class, critical time period, etc. Default=-1.
#' @param aalbin The length bin for which multi-panel plots of the fit to
#' conditional age-at-length data will be produced for all years. Useful to
#' see if growth curves are ok, or to see the information on year classes move
#' through the conditional data. Default=-1.
#' @param aalresids Plot the full set of conditional age-at-length Pearson
#' residuals? Turn to FALSE if plots are taking too long and you don't want
#' them.
#' @param maxneff The maximum value to include on plots of input and effective
#' sample size. Occasionally a calculation of effective N blows up to very
#' large numbers, rendering it impossible to observe the relationship for other
#' data. Default=5000.
#' @param cohortlines Optional vector of birth years for cohorts for which to
#' add growth curves to numbers at length bubble plots. Default=c().
#' @param smooth Add loess smoother to observed vs. expected index plots and
#' input vs. effective sample size? Default=TRUE.
#' @param showsampsize Display sample sizes on composition plots? Default=TRUE.
#' @param showeffN Display effective sample sizes on composition plots?
#' Default=TRUE.
#' @param sampsizeline show line for input sample sizes on top of conditional
#' age-at-length plots (TRUE/FALSE, still in development)
#' @param effNline show line for effective sample sizes on top of conditional
#' age-at-length plots (TRUE/FALSE, still in development)
#' @param showlegend Display legends in various plots?
#' @template pwidth
#' @template pheight
#' @template pheight_tall
#' @template punits
#' @template ptsize
#' @template res
#' @template mainTitle
#' @template cex.main
#' @param selexlines Vector controlling which lines should be shown on
#' selectivity plots if the model includes retention. Default=1:5.
#' @param rows Number of rows to use for single panel plots. Default=1.
#' @param cols Number of columns to use for single panel plots. Default=1.
#' @param maxrows Maximum number of rows to for multi-panel plots.
#' @param maxcols Maximum number of columns for multi-panel plots.
#' @param maxrows2 Maximum number of rows for conditional age-at-length
#' multi-panel plots.
#' @param maxcols2 Maximum number of rows for conditional age-at-length
#' multi-panel plots.
#' @param andrerows Number of rows of Andre's conditional age-at-length plots
#' within each page.
#' @param tagrows Number of rows for tagging-related plots.
#' @param tagcols Number of columns for tagging-related plots.
#' @param parrows Number of rows for parameter distribution plots.
#' @param parcols Number of columns for parameter distribution plots.
#' @param fixdims Control whether multi-panel plots all have dimensions equal
#' to maxrows by maxcols, or resized within those limits to fit number of
#' plots. Default=TRUE.
#' @param new Open a new window or add to existing plot windows. Default=TRUE.
#' @param SSplotDatMargin Size of right-hand margin in data plot (may be too
#' small if fleet names are long)
#' @param filenotes Optional vector of character strings to be added to intro
#' HTML page (if created) with notes about the model.
#' @param catchasnumbers Is catch input in numbers instead of biomass?
#' Default=F.
#' @param catchbars show catch by fleet as barplot instead of stacked polygons
#' (default=TRUE)
#' @template legendloc
#' @param minyr First year to show in time-series and time-varying plots
#' @param maxyr Last year to show in time-series and time-varying plots. This
#' can either be an alternative to, or redundant with, the forecastplot input.
#' @param sexes Which sexes to show in composition plots. Default="all".
#' @param scalebins Rescale expected and observed proportions in composition
#' plots by dividing by bin width for models where bins have different widths?
#' Caution!: May not work correctly in all cases.
#' @param scalebubbles scale data-only bubbles by sample size, not just
#' proportion within sample? Default=FALSE.
#' @param tslabels Either NULL to have default labels for timeseries plots or
#' a vector of appropriate length (currently 11) with labels for each figure
#' @param catlabels Either NULL to have default labels for catch plots or
#' a vector of appropriate length (currently 10) with labels for each figure
#' @param maxsize The size of the largest bubble in the datasize
#' plot. Default is 1.0.
#' @param showmle Show MLE estimate and asymptotic variance estimate with blue
#' lines in the parameter distribution plots?
#' @param showprior Show prior distribution as black line in the parameter
#' distribution plots?
#' @param showpost Show posterior distribution as bar graph in parameter
#' distribution plots (requires MCMC results to be available in `replist`)?
#' @param showinit Show initial value as red triangle in the parameter
#' distribution plots?
#' @param showdev Include devs in the parameter distribution plots?
#' @param fitrange Fit range in parameter distribution plots tightly around MLE
#' and posterior distributions instead of full parameter range?
#' @param \dots Additional arguments that will be passed to some subfunctions.
#' @author Ian Stewart, Ian Taylor
#' @export
#' @seealso [SS_output()], [SSplotBiology()],
#' [SSplotCatch()], [SSplotComps()],
#' [SSplotDiscard()], [SSplotIndices()],
#' [SSplotMnwt()], [SSplotNumbers()],
#' [SSplotRecdevs()], [SSplotSelex()],
#' [SSplotSpawnrecruit()], [SSplotSPR()],
#' [SSplotTags()], [SSplotTimeseries()],
#' [SSplotYield()]
#' @references Walters, Hilborn, and Christensen, 2008, Surplus production
#' dynamics in declining and recovering fish populations. Can. J. Fish. Aquat.
#' Sci. 65: 2536-2551.
SS_plots <-
function(replist = NULL, plot = 1:26, pdf = FALSE, png = TRUE, html = png,
printfolder = "plots", dir = "default", fleets = "all", areas = "all",
fleetnames = "default", fleetcols = "default", fleetlty = 1, fleetpch = 1,
lwd = 1, areacols = NULL, areanames = "default",
verbose = TRUE, uncertainty = TRUE, forecastplot = FALSE,
datplot = TRUE, Natageplot = TRUE, samplesizeplots = TRUE, compresidplots = TRUE,
comp.yupper = 0.4,
sprtarg = "default", btarg = "default", minbthresh = "default", pntscalar = NULL,
bub.scale.pearson = 1.5, bub.scale.dat = 3, pntscalar.nums = 2.6,
pntscalar.tags = 2.6, minnbubble = 8, aalyear = -1, aalbin = -1, aalresids = TRUE,
maxneff = 5000, cohortlines = c(), smooth = TRUE, showsampsize = TRUE,
showeffN = TRUE, sampsizeline = FALSE, effNline = FALSE,
showlegend = TRUE,
pwidth = 6.5, pheight = 4.0, pheight_tall = 6.5,
punits = "in", ptsize = 10, res = 300,
mainTitle = FALSE, cex.main = 1, selexlines = 1:6, rows = 1, cols = 1,
maxrows = 6, maxcols = 4, maxrows2 = 4, maxcols2 = 4, andrerows = 4,
tagrows = 3, tagcols = 3, parrows = 4, parcols = 2, fixdims = TRUE, new = TRUE,
SSplotDatMargin = 8, filenotes = NULL, catchasnumbers = NULL, catchbars = TRUE,
legendloc = "topleft", minyr = -Inf, maxyr = Inf, sexes = "all", scalebins = FALSE,
scalebubbles = FALSE, tslabels = NULL, catlabels = NULL, maxsize = 1.0,
showmle = TRUE, showpost = TRUE, showprior = TRUE, showinit = TRUE, showdev = FALSE,
fitrange = FALSE, ...) {
flush.console()
# label table is a step toward internationalization of the code
# in the future, this could be read from a file, or we could have multiple columns
# in the table to choose from
if (is.null(replist) || !is.list(replist) || !"nfleets" %in% names(replist)) {
stop(
"The input 'replist' should refer to an R object created by",
" the function 'SS_output'."
)
}
# get quantities from the big list
nfleets <- replist[["nfleets"]]
nfishfleets <- replist[["nfishfleets"]]
nareas <- replist[["nareas"]]
nseasons <- replist[["nseasons"]]
timeseries <- replist[["timeseries"]]
lbins <- replist[["lbins"]]
inputs <- replist[["inputs"]]
endyr <- replist[["endyr"]]
SS_version <- replist[["SS_version"]]
SS_versionNumeric <- replist[["SS_versionNumeric"]]
StartTime <- replist[["StartTime"]]
Files_used <- replist[["Files_used"]]
FleetNames <- replist[["FleetNames"]]
rmse_table <- replist[["rmse_table"]]
comp_data_exists <- replist[["comp_data_exists"]]
# check for internal consistency
if (pdf & png) {
stop("Inputs 'pdf' and 'png' are mututally exclusive. You need to set one of them to FALSE")
}
if (html & !png) {
stop("You can't set 'html=TRUE' without also setting 'png=TRUE'")
}
if (uncertainty & !inputs[["covar"]]) {
message("covar information unavailable, changing 'uncertainty' to FALSE")
uncertainty <- FALSE
}
if (forecastplot & max(timeseries[["Yr"]] > endyr + 1) == 0) {
message("Changing 'forecastplot' input to FALSE because all years up to endyr+1 are included by default")
forecastplot <- FALSE
}
# derived quantities
if (fleets[1] == "all") {
fleets <- 1:nfleets
} else {
if (length(intersect(fleets, 1:nfleets)) != length(fleets)) {
return("Input 'fleets' should be 'all' or a vector of values between 1 and nfleets.")
}
}
if (areas[1] == "all") {
areas <- 1:nareas
} else {
if (length(intersect(areas, 1:nareas)) != length(areas)) {
return("Input 'areas' should be 'all' or a vector of values between 1 and nareas.")
}
}
if (verbose) {
message("Finished defining objects")
}
# set fleet-specific names, and plotting parameters
if (fleetnames[1] == "default") {
fleetnames <- FleetNames
}
if (fleetcols[1] == "default") {
fleetcols <- rich.colors.short(nfishfleets)
if (nfishfleets > 2) fleetcols <- rich.colors.short(nfishfleets + 1)[-1]
}
if (length(fleetlty) < nfishfleets) {
fleetlty <- rep(fleetlty, nfishfleets)
}
if (length(fleetpch) < nfishfleets) {
fleetpch <- rep(fleetpch, nfishfleets)
}
# set default area-specific colors if not specified
areacols <- get_areacols(areacols, nareas)
#### prepare for plotting
# number of plot groups
nplots <- length(intersect(1:50, plot))
# make plot window (hopefully no-longer operating system specific)
if (nplots > 0 & !png & !pdf & new) {
### Note: the following line has been commented out because it was identified
### by Brian Ripley as "against CRAN policies".
# if(exists(".SavedPlots",where=1)) rm(.SavedPlots,pos=1)
dev.new(width = pwidth, height = pheight, pointsize = ptsize, record = TRUE)
}
if (nplots > 0 & !new) {
if (verbose) {
message("Adding plots to existing plot window. Plot history not erased.")
}
}
### deal with directories in which to create PNG or PDF files
if (dir == "default") {
# directory within which printfolder will be created
# by default it is assumed to be the location of the model files
dir <- inputs[["dir"]]
}
if (png | pdf) {
# get info on directory where subfolder will go
# (typically folder with model output files)
dir.isdir <- file.info(dir)$isdir
# create directory
if (is.na(dir.isdir) | !dir.isdir) {
message("Directory doesn't exist, attempting to create:\n", dir)
dir.create(dir)
}
# test again (even though failure to create dir should have already caused error)
# get info on directory where subfolder will go
# (typically folder with model output files)
dir.isdir <- file.info(dir)$isdir
# create
if (is.na(dir.isdir) | !dir.isdir) {
stop("Not able to create directory:\n", dir, "\n")
}
}
plotdir <- "default" # dummy value passed to functions that ignore it if png=FALSE
if (png) {
# add subdirectory for PNG and HTML files if that option is chosen
# close any in-process plots still open
graphics.off()
# figure out path to where PNG files will go
plotdir <- file.path(dir, printfolder)
plotdir.isdir <- file.info(plotdir)$isdir
if (is.na(plotdir.isdir) | !plotdir.isdir) {
dir.create(plotdir)
}
if (verbose) {
message(
"Plots will be written to PNG files in the directory:\n ",
plotdir
)
}
# get info on any older plots inside the plotdir directory
csv.files <- grep("plotInfoTable.+csv", dir(plotdir), value = TRUE)
if (length(csv.files) > 0) {
StartTimes.old <- NULL
for (ifile in seq_along(csv.files)) {
plotInfo.old <- read.csv(file.path(plotdir, csv.files[ifile]),
stringsAsFactors = FALSE
)
StartTimes.old <- c(StartTimes.old, unique(plotInfo.old[["StartTime"]]))
}
if (any(StartTimes.old != StartTime)) {
# if there are plots that are older than those from the current model,
# rename the directory to something containing the older model start time
StartTimeName <- gsub(":", ".", StartTimes.old[1], fixed = TRUE)
StartTimeName <- gsub(" ", "_", StartTimeName, fixed = TRUE)
StartTimeName <- gsub("._", "_", StartTimeName, fixed = TRUE)
plotdir.old <- file.path(dir, paste0("plots_", StartTimeName))
message(
"NOTE: the directory\n ",
plotdir,
"\n contains plots from a previous model run, renaming to\n ",
plotdir.old
)
file.rename(plotdir, plotdir.old)
# create a new, empty directory for the new plots
dir.create(plotdir)
}
}
}
# create PDF file if requested
if (pdf) {
pdffile <- file.path(dir, paste0(
"SS_plots_",
format(Sys.time(), "%d-%m-%Y_%H.%M"),
".pdf"
))
pdf(file = pdffile, width = pwidth, height = pheight)
if (verbose) {
message("PDF file with plots will be:", pdffile)
}
}
# blank table to store plot info
plotInfoTable <- NULL
if (new & !png) {
# make multi-panel plot if requested (not available for PNG files)
par(mfcol = c(rows, cols))
}
if (pdf) {
mar0 <- par()$mar # current margins
par(mar = rep(0, 4))
plot(0, type = "n", xlab = "", ylab = "", axes = FALSE, xlim = c(0, 1), ylim = c(0, 1))
y <- 0.9
ystep <- -.05
text(0, y, "Plots created using the 'r4ss' package in R", pos = 4)
y <- y + ystep
text(0, y, paste("Stock Synthesis version:", substr(SS_version, 1, 9)), pos = 4)
y <- y + ystep
text(0, y, StartTime, pos = 4)
y <- y + ystep
Files2 <- strsplit(Files_used, " ")[[1]]
text(0, y, paste(Files2[[1]], Files2[2]), pos = 4)
y <- y + ystep
text(0, y, paste(Files2[[3]], Files2[4]), pos = 4)
if (!is.null(filenotes)) {
y <- y + ystep
text(0, y, "Notes:", pos = 4)
for (i in seq_along(filenotes)) {
y <- y + ystep
text(0, y, filenotes[i], pos = 4)
}
}
par(mar = mar0) # replace margins
}
mar0 <- par()$mar # current inner margins
oma0 <- par()$oma # current outer margins
##########################################
# Biology plots (mean weight, maturity, fecundity, spawning output)
# and Time-varying growth
#
igroup <- 1
if (igroup %in% plot | length(cohortlines) > 0) {
if (verbose) {
message("Starting biology plots (group ", igroup, ")")
}
plotinfo <- SSplotBiology(
replist = replist,
areacols = areacols,
forecast = forecastplot, minyr = minyr, maxyr = maxyr,
plot = !png, print = png,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res, mainTitle = mainTitle,
cex.main = cex.main, plotdir = plotdir
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
}
##########################################
# Selectivity and retention plots
#
igroup <- 2
if (igroup %in% plot) {
if (verbose) {
message("Starting selectivity and retention plots (group ", igroup, ")")
}
selexinfo <-
SSplotSelex(
replist = replist, selexlines = selexlines,
fleets = fleets, fleetnames = fleetnames,
minyr = minyr, maxyr = maxyr,
plot = !png, print = png,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res, cex.main = cex.main,
plotdir = plotdir,
mainTitle = mainTitle,
)
plotinfo <- selexinfo[["plotinfo"]]
if (!is.null(plotinfo)) {
plotInfoTable <- rbind(plotInfoTable, plotinfo)
}
# add plots of unavailable (cryptic) spawning output
if (FALSE) {
# needs revision to work in SS 3.30 and models with multiple birth seasons
plotinfo <-
SSunavailableSpawningOutput(
replist = replist,
plot = !png, print = png,
plotdir = plotdir,
pwidth = pwidth, pheight = pheight,
punits = punits, res = res,
ptsize = ptsize, cex.main = cex.main
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
}
} # end if igroup in plot or print
##########################################
# Basic time series
#
igroup <- 3
if (igroup %in% plot) {
if (verbose) {
message("Starting timeseries plots (group ", igroup, ")")
}
# which subplots to make (those for spawn bio first)
subplot_list <- c(7:10, 1:6, 11:15)
# loop over subplots
for (isubplot in subplot_list) {
# loop over add forecast or not
for (doforecast in unique(c(FALSE, forecastplot))) {
# subset of plots for which uncertainty might be added
if (isubplot %in% c(7, 9, 11)) {
for (douncertainty in unique(c(FALSE, uncertainty))) { # add uncertainty or not
plotinfo <-
SSplotTimeseries(
replist = replist,
subplot = isubplot,
areas = areas,
areacols = areacols,
areanames = areanames,
forecastplot = doforecast,
uncertainty = douncertainty,
plot = !png, print = png,
verbose = verbose,
btarg = btarg,
minbthresh = minbthresh,
minyr = minyr, maxyr = maxyr,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res,
mainTitle = mainTitle,
cex.main = cex.main,
labels = tslabels,
plotdir = plotdir
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
} # end loop over uncertainty or not
} else { # these plots don't have the option for uncertainty
plotinfo <-
SSplotTimeseries(
replist = replist,
subplot = isubplot,
areas = areas,
areacols = areacols,
areanames = areanames,
forecastplot = doforecast,
uncertainty = FALSE,
plot = !png, print = png,
verbose = verbose,
btarg = btarg,
minbthresh = minbthresh,
minyr = minyr, maxyr = maxyr,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res,
mainTitle = mainTitle,
cex.main = cex.main,
labels = tslabels,
plotdir = plotdir
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
}
}
} # end loop over timeseries subplots
### add plot of Summary F
# first get vector of years
yrs <- replist[["startyr"]]:replist[["endyr"]]
yrs <- yrs[yrs >= minyr & yrs <= maxyr]
# now run plot function
plotinfo <- SSplotSummaryF(
replist = replist,
yrs = yrs,
uncertainty = uncertainty,
plot = !png, print = png,
verbose = verbose,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res,
plotdir = plotdir
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
### add plot of Dynamic B0
if (is.null(replist[["Dynamic_Bzero"]])) {
message("Skipping dynamic B0 plot because output not available")
} else {
# first get vector of years
yrs <- replist[["startyr"]]:(replist[["endyr"]] + 1)
yrs <- yrs[yrs >= minyr & yrs <= maxyr]
# check for uncertainty in Dynamic B0 estimates
Dyn_Bzero_uncertainty <-
any(grepl("Dyn_Bzero", replist[["derived_quants"]][["Label"]]))
# now run plot function
plotinfo <- SSplotDynamicB0(
replist = replist,
yrs = yrs,
uncertainty = (uncertainty & Dyn_Bzero_uncertainty),
plot = !png, print = png,
verbose = verbose,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res,
plotdir = plotdir
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
}
} # end if igroup in plot or print
##########################################
# Recruitment deviation plots
#
igroup <- 4
if (igroup %in% plot) {
if (verbose) {
message("Starting recruitment deviation plots (group ", igroup, ")")
}
plotinfo <-
SSplotRecdevs(
replist = replist,
plot = !png, print = png,
forecastplot = forecastplot,
uncertainty = uncertainty,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res, cex.main = cex.main,
plotdir = plotdir
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
if (nareas > 1 & nseasons > 1) {
plotinfo <-
SSplotRecdist(
replist = replist,
plot = !png, print = png,
verbose = verbose,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res, cex.main = cex.main,
plotdir = plotdir
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
}
} # end if igroup in plot or print
##########################################
# Estimating recruitment bias adjustment plots
#
igroup <- 5
if (igroup %in% plot) {
if (uncertainty) {
if (verbose) {
message("Starting estimation of recruitment bias adjustment and associated plots (group ", igroup, ")")
}
if (is.numeric(rmse_table[["RMSE"]])) {
if (max(rmse_table[["RMSE"]]) > 0) {
temp <-
SS_fitbiasramp(
replist = replist,
plot = !png, print = png,
twoplots = FALSE,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res, cex.main = cex.main,
plotdir = plotdir
)
plotinfo <- temp[["plotinfo"]]
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
} else {
message("Skipping bias adjustment fit because root mean squared error of recruit devs is 0.")
}
} else {
message(
"skipping bias adjustment fit because\n",
"input list element 'rmse_table' has non-numeric 'RMSE' column"
)
}
} else {
if (verbose) {
message("Skipping estimation of recruitment bias adjustment (group ", igroup, ") because uncertainty=FALSE")
}
}
} # end if igroup in plot or print
##########################################
# spawner-recruit curve
#
igroup <- 6
if (igroup %in% plot) {
if (verbose) {
message("Starting spawner-recruit curve plot (group ", igroup, ")")
}
plotinfo <-
SSplotSpawnrecruit(
replist = replist,
plot = !png, print = png,
virg = TRUE, # add point on curve at equilibrium values (B0,R0)
init = FALSE, # add point on curve at initial values (B1,R1)
pwidth = pwidth, pheight = pheight_tall, punits = punits,
ptsize = ptsize, res = res,
plotdir = plotdir
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
} # end if igroup in plot or print
##########################################
# time series of catch
#
igroup <- 7
if (igroup %in% plot) {
if (verbose) {
message("Starting catch plots (group ", igroup, ")")
}
temp <-
SSplotCatch(
replist = replist,
plot = !png, print = png,
fleetnames = fleetnames,
fleetlty = fleetlty,
fleetpch = fleetpch,
fleetcols = fleetcols,
minyr = minyr, maxyr = maxyr,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res,
cex.main = cex.main,
catchasnumbers = catchasnumbers,
order = "default",
catchbars = catchbars,
labels = catlabels,
legendloc = legendloc,
plotdir = plotdir,
verbose = verbose
)
plotinfo <- temp[["plotinfo"]]
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
} # end if igroup in plot or print
##########################################
# SPR and fishing intensity plots
#
igroup <- 8
if (igroup %in% plot) {
if (verbose) {
message("Starting SPR plots (group ", igroup, ")")
}
plotinfo <-
SSplotSPR(
replist = replist,
plot = !png, print = png,
uncertainty = uncertainty,
sprtarg = sprtarg, btarg = btarg,
pwidth = pwidth, pheight = pheight, pheight_tall = pheight_tall,
punits = punits,
ptsize = ptsize, res = res, cex.main = cex.main,
plotdir = plotdir
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
} # end if igroup in plot or print
##########################################
# discard fit plots (if present)
#
igroup <- 9
if (igroup %in% plot) {
if (!is.null(replist[["discard"]]) &&
!is.na(replist[["discard"]][[1]][1]) &&
nrow(replist[["discard"]]) > 0) {
if (verbose) {
message("Starting discard plot (group ", igroup, ")")
}
plotinfo <-
SSplotDiscard(
replist = replist,
plot = !png, print = png,
fleets = fleets,
fleetnames = fleetnames,
datplot = datplot,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res, cex.main = cex.main,
plotdir = plotdir
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
} else {
if (verbose) {
message("Skipping discard plot (group ", igroup, ") because no discard data")
}
}
} # end if igroup in plot or print
##########################################
# mean body weight (if present)
#
igroup <- 10
if (igroup %in% plot) {
if (!is.null(replist[["mnwgt"]]) &&
!is.na(replist[["mnwgt"]][[1]][1]) &&
nrow(replist[["mnwgt"]]) > 0) {
if (verbose) {
message("Starting mean body weight plot (group ", igroup, ")")
}
plotinfo <-
SSplotMnwt(
replist = replist,
plot = !png, print = png,
fleets = fleets,
fleetnames = fleetnames,
datplot = datplot,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res, cex.main = cex.main,
plotdir = plotdir
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
} else {
if (verbose) {
message("Skipping mean weight plot (group ", igroup, ") because no mean weight data")
}
}
} # end if igroup in plot or print
##########################################
# Index plots
#
igroup <- 11
if (igroup %in% plot) {
if (!is.null(dim(replist[["cpue"]]))) {
if (verbose) {
message("Starting index plots (group ", igroup, ")")
}
plotinfo <- SSplotIndices(
replist = replist,
fleets = fleets,
fleetnames = fleetnames,
fleetcols = fleetcols,
plot = !png, print = png,
datplot = datplot,
pwidth = pwidth, pheight = pheight, punits = punits,
ptsize = ptsize, res = res,
mainTitle = mainTitle, cex.main = cex.main,
plotdir = plotdir,
minyr = minyr,
maxyr = maxyr
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
} else {
if (verbose) {
message(
"Skipping index plots (group ", igroup,
") because no indices in model (or are not reported)"
)
}
}
} # end if igroup in plot or print
##########################################
# Numbers at age plots
#
igroup <- 12
if (igroup %in% plot) {
if (!is.null(replist[["natage"]])) {
if (verbose) {
message("Starting numbers at age plots (group ", igroup, ")")
}
plotinfo <-
SSplotNumbers(
replist = replist,
areas = areas,
areanames = areanames,
areacols = areacols,
pntscalar = pntscalar.nums,
bublegend = showlegend,
plot = !png, print = png,
pwidth = pwidth, pheight = pheight_tall, punits = punits,
ptsize = ptsize, res = res,
mainTitle = mainTitle, cex.main = cex.main,
plotdir = plotdir
)
if (!is.null(plotinfo)) {
plotInfoTable <- rbind(plotInfoTable, plotinfo)
}
} else {
message(
"Skipping numbers plots (group ", igroup,
") because numbers-at-age table not included in output"
)
# end check for numbers-at-age table available
}
} # end if igroup in plot or print
##########################################
# Composition data plots
#
# use of SSplotcomps function to make composition plots
if (is.null(comp_data_exists) || !comp_data_exists) {
message("No composition data, skipping all composition plots")
} else {
lenCompDatGroup <- 13
ageCompDatGroup <- 14
condCompDatGroup <- 15
if (!datplot) {
if (length(intersect(
c(lenCompDatGroup, ageCompDatGroup, condCompDatGroup),
plot
)) > 0) {
message(
"Skipping plot groups ",
lenCompDatGroup,
"-",
condCompDatGroup,
" (comp data without fit) because input 'datplot=FALSE'"
)
}
} else {
if (lenCompDatGroup %in% plot) # data only aspects
{
if (verbose) {
message("Starting length comp data plots (group ", lenCompDatGroup, ")")
}
# length comp polygon and bubble plots
plotinfo <-
SSplotComps(
replist = replist, datonly = TRUE, kind = "LEN", bub = TRUE, verbose = verbose, fleets = fleets,
fleetnames = fleetnames,
samplesizeplots = samplesizeplots, showsampsize = showsampsize, showeffN = FALSE,
minnbubble = minnbubble, pntscalar = pntscalar, cexZ1 = bub.scale.dat,
bublegend = showlegend,
maxrows = maxrows, maxcols = maxcols, fixdims = fixdims, rows = rows, cols = cols,
plot = !png, print = png,
plotdir = plotdir, mainTitle = mainTitle, cex.main = cex.main,
sexes = sexes, yupper = comp.yupper,
scalebins = scalebins, scalebubbles = scalebubbles,
pwidth = pwidth, pheight = pheight_tall, punits = punits,
ptsize = ptsize, res = res,
...
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
# length comp sex ratios (data only, for 2-sex models only)
if (replist[["nsexes"]] == 2) {
plotinfo <-
SSplotSexRatio(
replist = replist,
datonly = TRUE,
kind = "LEN", bub = TRUE, verbose = verbose, fleets = fleets,
fleetnames = fleetnames,
linescol = 0, # turn off line showing expected value
samplesizeplots = samplesizeplots, showsampsize = showsampsize, showeffN = showeffN,
minnbubble = minnbubble, pntscalar = pntscalar, cexZ1 = bub.scale.pearson,
bublegend = showlegend,
maxrows = maxrows, maxcols = maxcols, fixdims = fixdims, rows = rows, cols = cols,
plot = !png, print = png, smooth = smooth, plotdir = plotdir,
maxneff = maxneff, mainTitle = mainTitle, cex.main = cex.main,
cohortlines = cohortlines,
# sexes=sexes,
# yupper=comp.yupper,
scalebins = scalebins,
pwidth = pwidth, pheight = pheight_tall, punits = punits,
ptsize = ptsize, res = res,
...
)
if (!is.null(plotinfo)) plotInfoTable <- rbind(plotInfoTable, plotinfo)
}
# size comp polygon and bubble plots (data only)
for (sizemethod in sort(unique(replist[["sizedbase"]][["method"]]))) {
plotinfo <-
SSplotComps(
replist = replist, datonly = TRUE, kind = "SIZE", sizemethod = sizemethod,
bub = TRUE, verbose = verbose, fleets = fleets, fleetnames = fleetnames,
samplesizeplots = samplesizeplots, showsampsize = showsampsize, showeffN = FALSE,
minnbubble = minnbubble, pntscalar = pntscalar, cexZ1 = bub.scale.dat,
bublegend = showlegend,
maxrows = maxrows, maxcols = maxcols, fixdims = fixdims, rows = rows, cols = cols,
plot = !png, print = png,
plotdir = plotdir, mainTitle = mainTitle, cex.main = cex.main,
sexes = sexes, yupper = comp.yupper,
scalebins = scalebins, scalebubbles = scalebubbles,
pwidth = pwidth, pheight = pheight_tall, punits = punits,
ptsize = ptsize, res = res,