Nothing
SSplotTimeseries <-
function(replist,subplot,add=FALSE,areas="all",
areacols="default",areanames="default",
forecastplot=TRUE,uncertainty=TRUE,bioscale="default",
minyr=NULL,maxyr=NULL,
plot=TRUE,print=FALSE,plotdir="default",verbose=TRUE,
btarg="default",minbthresh="default",xlab="Year",
labels=NULL,
pwidth=7,pheight=7,punits="in",res=300,ptsize=12,cex.main=1)
{
# individual function for plotting time series of total or summary biomass
# subplot1 = total biomass total all areas
# subplot2 = total biomass by area
# subplot3 = total biomass in all areas in spawning season
# subplot4 = summary biomass total all areas
# subplot5 = summary biomass by area
# subplot6 = summary biomass in all areas in spawning season
# subplot7 = spawning biomass total (with or without uncertainty)
# subplot8 = spawning biomass by area
# subplot9 = spawning depletion total (with or without uncertainty)
# subplot10 = spawning depletion by area
# subplot11 = recruitment total (with or without uncertainty)
# subplot12 = recruitment by area
# subplot13 = fraction of recruitment by area
# subplot14 = recruitment by birth season
# subplot15 = fraction of recruitment by birth season
if(missing(subplot)) stop("'subplot' input required")
if(length(subplot)>1) stop("function can only do 1 subplot at a time")
pngfun <- function(file,caption=NA){
png(filename=file,width=pwidth,height=pheight,
units=punits,res=res,pointsize=ptsize)
plotinfo <- rbind(plotinfo,data.frame(file=file,caption=caption))
return(plotinfo)
}
plotinfo <- NULL
# get values from replist
SS_versionshort <- replist$SS_versionshort
timeseries <- replist$timeseries
nseasons <- replist$nseasons
spawnseas <- replist$spawnseas
birthseas <- replist$birthseas
startyr <- replist$startyr
endyr <- replist$endyr
nsexes <- replist$nsexes
nareas <- replist$nareas
derived_quants <- replist$derived_quants
FecPar2 <- replist$FecPar2
B_ratio_denominator <- replist$B_ratio_denominator
seasfracs <- replist$seasfracs
recruitment_dist <- replist$recruitment_dist
if(btarg=="default") btarg <- replist$btarg
if(minbthresh=="default") minbthresh <- replist$minbthresh
if(areacols[1]=="default"){
areacols <- rich.colors.short(nareas)
if(nareas > 2) areacols <- rich.colors.short(nareas+1)[-1]
}
if(!is.null(birthseas)){
nbirthseas <- length(birthseas)
seascols <- rich.colors.short(nbirthseas)
if(nbirthseas > 2) seascols <- rich.colors.short(nbirthseas+1)[-1]
}
# temporary fix for SS_output versions prior to 9/20/2010
if(is.null(B_ratio_denominator)) B_ratio_denominator <- 1
# check if spawning output rather than spawning biomass is plotted
if(plotdir=="default") plotdir <- replist$inputs$dir
if(FecPar2!=0) labels[5] <- labels[7]
# check area subsets
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(nareas>1 & areanames[1]=="default"){
areanames <- paste("area",1:nareas)
}
#scaling factor for single sex models
if(bioscale=="default"){
if(nsexes==1) bioscale <- 0.5 else bioscale <- 1
}
# modifying data to subset for a single season
ts <- timeseries
if(nseasons>1){
if(SS_versionshort=="SS-V3.11"){
ts$YrSeas <- ts$Yr + (ts$Seas-1)/nseasons
}else{
ts$YrSeas <- ts$Yr + seasfracs
}
}else{
ts$YrSeas <- ts$Yr
}
# warn about spawning season--seems to no longer be necessary now that title
# is update for to reflect spawning season
## if(spawnseas>1 & subplot %in% c(3,6,7,8,9,10) ){
## cat("Note: spawning seems to be in season ",spawnseas,". Some plots will show only this season.\n",sep="")
## }
# define which years are forecast or not
ts$period <- "time"
ts$period[ts$Yr < startyr] <- "equilibria"
ts$period[ts$Yr > endyr+1] <- "fore"
if(!forecastplot) ts$period[ts$Yr > endyr + 1] <- "exclude"
# a function to make the plot
biofunc <- function(subplot){
# make the logical vector of which time-series entries to use
plot1 <- ts$Area==1 & ts$Era=="VIRG" # T/F for in area & is virgin value
plot2 <- ts$Area==1 & ts$period=="time" & ts$Era!="VIRG" # T/F for in area & not virgin value
plot3 <- ts$Area==1 & ts$period=="fore" & ts$Era!="VIRG" # T/F for in area & not virgin value
if(subplot %in% c(3,6,7,9)){
plot1 <- ts$Area==1 & ts$Era=="VIRG" & ts$Seas == spawnseas # T/F for in area & is virgin value
plot2 <- ts$Area==1 & ts$period=="time" & ts$Era!="VIRG" & ts$Seas == spawnseas # T/F for in area & not virgin value
plot3 <- ts$Area==1 & ts$period=="fore" & ts$Era!="VIRG" & ts$Seas == spawnseas # T/F for in area & not virgin value
}
# switch for which column of the TIME_SERIES table is being plotted
# subplot1,2,3 = total biomass
if(subplot %in% 1:3){
yvals <- ts$Bio_all
ylab <- labels[1]
if(subplot==3){ylab <- paste(labels[2],spawnseas)}
}
# subplot4,5,6 = summary biomass
if(subplot %in% 4:6){
yvals <- ts$Bio_smry
ylab <- labels[3]
if(subplot==6){ylab <- paste(labels[4],spawnseas)}
}
# subplot7&8 = spawning biomass
if(subplot %in% 7:8){
yvals <- bioscale*ts$SpawnBio
ylab <- labels[5]
}
# subplot9&10 = spawning depletion
if(subplot %in% 9:10){
# yvals for spatial models are corrected later within loop over areas
yvals <- ts$SpawnBio/ts$SpawnBio[!is.na(ts$SpawnBio)][1]
ylab <- labels[6]
}
# subplot11&12 = recruitment
if(subplot %in% 11:15){
yvals <- ts$Recruit_0
ylab <- labels[8]
}
# change ylab to represent fractions for those plots
if(subplot %in% c(13,15)) ylab <- labels[9]
# title initially set equal to y-label
main=ylab
# birth season-related calculations
yrshift <- 0 # years of shift for fish spawning to next birth season
if(!is.null(birthseas) && max(birthseas) < spawnseas){
# case where fish are born in the year after spawning
yrshift <- 1
}
if(!is.null(birthseas) && nbirthseas > 1){
if(subplot==11){
# sum total recruitment across birth seasons
for(y in ts$Yr){
yvals[ts$Yr==y & ts$Seas==1] <- sum(yvals[ts$Yr==y],na.rm=TRUE)
yvals[ts$Yr==y & ts$Seas >1] <- 0
}
}
if(subplot==15){
# sum total recruitment across birth seasons
for(y in ts$Yr){
yvals[ts$Yr==y] <- yvals[ts$Yr==y]/sum(yvals[ts$Yr==y],na.rm=TRUE)
}
}
if(subplot %in% c(14,15)) main=paste(main,"by birth season")
}
# sum up total across areas if needed
if(nareas>1){
if(subplot %in% c(2,3,5,6,8,10,12,13)){
# these plots have separate lines for each area
main=paste(main,"by area")
}
if(subplot %in% c(1,4,7,11,13)){
# these plots have sum across areas
yvals2 <- rep(NA,length(ts$YrSeas))
for(iyr in 1:length(yvals)){
y <- ts$YrSeas[iyr]
yvals2[iyr] <- sum(yvals[ts$YrSeas==y])
}
if(subplot==13){
yvals <- yvals/yvals2
}else{
yvals <- yvals2
}
}
if(subplot==9){
# sum up total across areas differently for spawning depletion
yvals2 <- rep(NA,length(ts$YrSeas))
for(iyr in 1:length(yvals)){
y <- ts$YrSeas[iyr]
yvals[iyr] <- sum(ts$SpawnBio[ts$YrSeas==y])
}
yvals <- yvals/yvals[!is.na(yvals)][1] # total depletion
}
ymax <- max(yvals,1,na.rm=TRUE)
# correct ymax value for plot 10 (other plots may need it too)
if(subplot==10){
for(iarea in 1:nareas){
yvals <- ts$SpawnBio[ts$Area==iarea]/(ts$SpawnBio[ts$Area==iarea & ts$Seas == spawnseas][1])
ymax <- max(ymax,yvals,na.rm=TRUE)
}
}
}
if(subplot==10){
yvals[1] <- NA
}
if(forecastplot) main <- paste(main,"with forecast")
# calculating intervals around spawning biomass, depletion, or recruitment
# area specific confidence intervals?
if(uncertainty & subplot %in% c(7,9,11)){
main <- paste(main,"with ~95% asymptotic intervals")
if(!"SPB_Virgin" %in% derived_quants$LABEL){
cat("Skipping spawning biomass with uncertainty plot because 'SPB_Virgin' not in derived quantites.\n",
" Try changing 'min yr for Spbio_sdreport' in starter file to -1.\n")
}else{
# get subset of DERIVED_QUANTITIES
if(subplot==7){ # spawning biomass
stdtable <- derived_quants[grep("SPB_Virgin",derived_quants[,1]):(grep("Recr_Virgin",derived_quants[,1])-1),1:3]
# year as part of the LABEL string starting with 5th character
stdtable$Yr <- substring(stdtable$LABEL,5)
# filling in Virgin and Initial years as 2 and 1 years prior to following years
stdtable$Yr[1:2] <- as.numeric(stdtable$Yr[3])-(2:1) - yrshift
stdtable$Yr <- as.numeric(stdtable$Yr)
}
if(subplot==9){ # spawning depletion
stdtable <- derived_quants[substring(derived_quants$LABEL,1,6)=="Bratio",]
stdtable$Yr <- as.numeric(substring(stdtable$LABEL,8))
### these temporary fixes now replaced using "B_ratio_denominator"
## if(abs(stdtable$Value[1] - 4)<.1) bioscale <- 1/4 # temporary fix
## if(abs(stdtable$Value[1] - 2.5)<.1) bioscale <- 1/2.5 # temporary fix
bioscale <- B_ratio_denominator
}
if(subplot==11){ # recruitment
stdtable <- derived_quants[substring(derived_quants$LABEL,1,5)=="Recr_",]
stdtable <- stdtable[stdtable$LABEL!="Recr_Unfished",]
# year as the part of the LABEL string starting with 6th character
stdtable$Yr <- substring(stdtable$LABEL,6)
# filling in Virgin and Initial years as 2 and 1 years prior to following years
stdtable$Yr[1:2] <- as.numeric(stdtable$Yr[3])-(2:1)
stdtable$Yr <- as.numeric(stdtable$Yr) + yrshift
bioscale <- 1
}
# scaling and calculation of confidence intervals
v <- stdtable$Value * bioscale
std <- stdtable$StdDev * bioscale
if(subplot==11){
# assume recruitments have log-normal distribution
# from first principals (multiplicative survival probabilities)
stdtable$logint <- sqrt(log(1+(std/v)^2))
stdtable$lower <- exp(log(v) - 1.96*stdtable$logint)
stdtable$upper <- exp(log(v) + 1.96*stdtable$logint)
}else{ # assume normal distribution matching internal assumptions of ADMB
stdtable$upper <- v + 1.96*std
stdtable$lower <- pmax(v - 1.96*std, 0) # max of value or 0
}
if(max(stdtable$Yr) < max(floor(ts$YrSeas))){
cat(" !warning:\n",
" ",max(stdtable$Yr),"is last year with uncertainty in Report file, but",max(ts$YrSeas),"is last year of time series.\n",
" Consider changing starter file input for 'max yr for sdreport outputs' to -2\n")
}
}
}
# improved y-range for plot (possibly excluding time periods that aren't plotted)
# only works on single area models
if(nareas==1){
ymax <- max(yvals[plot1 | plot2 | plot3], na.rm=TRUE)
}
if(subplot%in%c(13,15)) ymax <- 1 # these plots show fractions
if(uncertainty & subplot %in% c(7,9,11)) ymax <- max(ymax,stdtable$upper, na.rm=TRUE)
if(print){ # if printing to a file
# adjust file names
filename <- main
filename <- gsub(",","",filename,fixed=TRUE)
filename <- gsub("~","",filename,fixed=TRUE)
filename <- gsub("%","",filename,fixed=TRUE)
if(forecastplot) filename <- paste(filename,"forecast")
if(uncertainty & subplot %in% c(5,7,9)) filename <- paste(filename,"intervals")
filename <- paste("ts",subplot," ",filename,".png",sep="")
filename <- paste(plotdir,filename,sep="")
# if(verbose) cat("printing plot to file:",filename,"\n")
plotinfo <- pngfun(file=filename,caption=main)
}
# move VIRG value from startyr-2 to startyr-1 to show closer to plot
# this one didn't work: if(exists("stdtable")) stdtable$Yr[stdtable$Yr %in% ts$Yr[plot1]] <- stdtable$Yr[stdtable$Yr %in% ts$Yr[plot1]]+1
ts$Yr[ts$Era=="VIRG"] <- ts$Yr[ts$Era=="VIRG"]+1
ts$YrSeas[ts$Era=="VIRG"] <- ts$YrSeas[ts$Era=="VIRG"]+1
# create an empty plot (if not adding to existing plot)
if(!add){
yrvals <- ts$YrSeas[ plot1 | plot2 | plot3]
# axis limits
if(is.null(minyr)) minyr <- min(yrvals)
if(is.null(maxyr)) maxyr <- max(yrvals)
xlim <- c(minyr,maxyr)
plot(yrvals,yvals[plot1 | plot2 | plot3],
type='n', xlab=xlab, ylim=c(0,ymax), ylab=ylab,
main=main, cex.main=cex.main,xlim=xlim)
abline(h=0,col="grey")
}
# add references points to plot of depletion
if(subplot %in% c(9,10))
{
addtarg <- function(){
if(btarg>0 & btarg<1){
abline(h=btarg,col="red")
text(max(startyr,minyr)+4,btarg+0.03,labels[10],adj=0)
}
if(minbthresh>0 & minbthresh<1){
abline(h=minbthresh,col="red")
text(max(startyr,minyr)+4,minbthresh+0.03,labels[11],adj=0)
}
}
addtarg()
}
# add references points to plot of abundance
if(subplot %in% 7:9)
{
addtarg <- function(){
if(btarg>1){
abline(h=btarg,col="red")
text(max(startyr,minyr)+4,btarg+0.02*diff(par()$usr[3:4]),
labels[10],adj=0)
}
if(minbthresh>1){
abline(h=minbthresh,col="red")
text(max(startyr,minyr)+4,minbthresh+0.02*diff(par()$usr[3:4]),
labels[11],adj=0)
}
}
addtarg()
}
if(subplot %in% 14:15){
# these plots show lines for each birth season,
# but probably won't work if there are multiple birth seasons and multiple areas
for(iseas in 1:nbirthseas){
s <- birthseas[iseas]
mycol <- seascols[iseas]
mytype <- "o" # overplotting points on lines for most time series
plot1 <- ts$Seas==s & ts$Era=="VIRG" # T/F for in seas & is virgin value
plot2 <- ts$Seas==s & ts$period=="time" & ts$Era!="VIRG" # T/F for in seas & not virgin value
plot3 <- ts$Seas==s & ts$period=="fore" & ts$Era!="VIRG" # T/F for in seas & is forecast
points(ts$Yr[plot1],yvals[plot1],pch=19, col=mycol) # filled points for virgin conditions
lines(ts$Yr[plot2],yvals[plot2],type=mytype,col=mycol) # open points and lines in middle
points(ts$Yr[plot3],yvals[plot3],pch=19, col=mycol) # filled points for forecast
}
legend("topright",legend=paste("Season",birthseas),lty=1,pch=1,col=seascols,bty="n")
}else{
# always loop over areas, but for plots with only one line,
# change vector of areas to equal 1.
if(subplot %in% c(1,4,7,9,11,14,15)) myareas <- 1 else myareas <- areas
for(iarea in myareas){ # loop over chosen areas
###
# subset for time period to allow different colors in plot
# plot1 = subset for equilibrium (virgin) values
# plot2 = subset for main timeseries
# plot3 = subset for forecast
###
if(subplot==10){
yvals <- ts$SpawnBio/(ts$SpawnBio[ts$Area==iarea & ts$Seas == spawnseas][1])
}
if(subplot %in% c(3,6,7,8,9,10)){
plot1 <- ts$Area==iarea & ts$Era=="VIRG" & ts$Seas == spawnseas # T/F for in area & is virgin value
plot2 <- ts$Area==iarea & ts$period=="time" & ts$Era!="VIRG" & ts$Seas == spawnseas # T/F for in area & not virgin value
plot3 <- ts$Area==iarea & ts$period=="fore" & ts$Era!="VIRG" & ts$Seas == spawnseas # T/F for in area & not virgin value
}else{
plot1 <- yvals>0 & ts$Area==iarea & ts$Era=="VIRG" # T/F for in area & is virgin value
plot2 <- yvals>0 & ts$Area==iarea & ts$period=="time" & ts$Era!="VIRG" # T/F for in area & not virgin value
plot3 <- yvals>0 & ts$Area==iarea & ts$period=="fore" & ts$Era!="VIRG" # T/F for in area & not virgin value
}
if(subplot %in% 9:10){
plot1 <- NULL
plot2[3] <- FALSE
}
mycol <- areacols[iarea]
mytype <- "o" # overplotting points on lines for most time series
if(subplot==11 & uncertainty) mytype <- "p" # just points without connecting lines if plotting recruitment with confidence intervals
if(!uncertainty){
points(ts$YrSeas[plot1],yvals[plot1],pch=19, col=mycol) # filled points for virgin conditions
lines( ts$YrSeas[plot2],yvals[plot2],type=mytype,col=mycol) # open points and lines in middle
points(ts$YrSeas[plot3],yvals[plot3],pch=19, col=mycol) # filled points for forecast
}else{
# add lines for confidence intervals areas if requested
# lines and points on integer years
points(ts$Yr[plot1],yvals[plot1],pch=19, col=mycol) # filled points for virgin conditions
lines( ts$Yr[plot2],yvals[plot2],type=mytype,col=mycol) # open points and lines in middle
points(ts$Yr[plot3],yvals[plot3],pch=19, col=mycol) # filled points for forecast
if(subplot %in% c(7,9,11)){
# subset years for confidence intervals
if(subplot==7){
plot1 <- stdtable$LABEL=="SPB_Virgin"
stdtable$Yr[plot1] <- stdtable$Yr[plot1]+yrshift
}
if(subplot==9){
plot1 <- stdtable$LABEL=="Bratio_Virgin" # note: this doesn't exist
}
if(subplot==11){
plot1 <- stdtable$LABEL=="Recr_Virgin"
stdtable$Yr[plot1] <- stdtable$Yr[plot1]+1 # shifting as in other cases to make Virgin year adjacent to first year of timeseries
}
plot2 <- stdtable$Yr %in% ts$Yr[plot2]
plot3 <- stdtable$Yr %in% ts$Yr[plot3]
plotall <- plot1 | plot2 | plot3
## stdtable$plot1 <- plot1
## stdtable$plot2 <- plot2
## stdtable$plot3 <- plot3
}
if(subplot %in% c(7,9)){
# add lines for main period
lines(stdtable$Yr[plot2], stdtable$upper[plot2], lty=2, col=mycol)
lines(stdtable$Yr[plot2], stdtable$lower[plot2], lty=2, col=mycol)
# add dashes for early period
points(stdtable$Yr[plot1]+1, stdtable$upper[plot1], pch="-", col=mycol) # +1 is because VIRG was shifted right 1 year
points(stdtable$Yr[plot1]+1, stdtable$lower[plot1], pch="-", col=mycol) # +1 is because VIRG was shifted right 1 year
# add dashes for forecast period
points(stdtable$Yr[plot3], stdtable$upper[plot3], pch="-", col=mycol)
points(stdtable$Yr[plot3], stdtable$lower[plot3], pch="-", col=mycol)
}
if(subplot==11){ # confidence intervals as error bars because recruitment is more variable
old_warn <- options()$warn # previous setting
options(warn=-1) # turn off "zero-length arrow" warning
arrows(x0=stdtable$Yr[plotall], y0=stdtable$lower[plotall], y1=stdtable$upper[plotall],
length=0.01, angle=90, code=3, col=mycol)
options(warn=old_warn) #returning to old value
}
} # end if uncertainty
} # end loop over areas
if(nareas>1 & subplot%in%c(2,3,5,6,8,10,12)) legend("topright",legend=areanames[areas],lty=1,pch=1,col=areacols[areas],bty="n")
} # end test for birthseason plots or not
if(verbose) cat(" finished time series subplot ",subplot,": ",main,"\n",sep="")
if(print) dev.off()
return(plotinfo)
} # end biofunc
# make plots
# for(iplot in subplot){ # doesn't work for more than one subplota at a time
skip <- FALSE
# plots 2, 5, 8, 10, and 12 are redundant for 1-area models
if(nareas==1 & subplot %in% c(2,5,8,10,12:13)) skip <- TRUE
# plots 3 and 6 are redundant for 1-season models
if(nseasons==1 & subplot %in% c(3,6)) skip <- TRUE
if(subplot %in% c(14:15) & (is.null(birthseas) || nbirthseas==1)) skip <- TRUE
if(!skip){
plotinfo <- biofunc(subplot=subplot)
if(!is.null(plotinfo)) plotinfo$category <- "Timeseries"
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.