R/plotExploitationRatesGG.R

In wStockhausen/rTCSAM2015: A set of tools to set up input for, and interact with output from, a TCSAM2015 model run, as well as to compare TCSAM2015 and rsimTCSAM runs.

#'
#'@title Plot time series of exploitation rates based on abundance (natural, discards, retained mortality) by fishery and sex.
#'
#'@description Plot time series of exploitation rates based on abundance (natural, discards, retained mortality) by fishery and sex.
#'
#'@param res - results list from a TCSAM2015 model run
#'@param showPlots - flag (T/F) to show plots
#'
#'@return list of ggplot objects
#'
#'@import ggplot2 
#'
#'@export
#'
plotExploitationRatesGG<-function(res,
                                  showPlot=TRUE){
    n<-reshape2::melt(res$pop.quants$n.xmsyz,value.name='N');
    n<-n[n$year<=res$mc$mxYr,];
    n<-reshape2::dcast(n,sex+year~.,value.var='N',fun.aggregate=sum);
    names(n)[3]<-'N';
    
    tm<-reshape2::melt(res$pop.quants$tmN.xmsyz,value.name='mortality');
    tm<-reshape2::dcast(tm,sex+year~.,value.var='mortality',fun.aggregate=sum);
    names(tm)[3]<-'mortality';
    tm$explrate<-tm$mortality/n$N;
    tm$type<-'total mortality';
    
    nm<-reshape2::melt(res$pop.quants$nmN.xmsyz,value.name='mortality');
    nm<-reshape2::dcast(nm,sex+year~.,value.var='mortality',fun.aggregate=sum);
    names(nm)[3]<-'mortality';
    nm$explrate<-nm$mortality/n$N;
#    nm$fishery<-'';
    nm$type<-'natural mortality';
    
    dfr<-NULL;    
    nFsh<-res$mc$nFsh;
    for (f in 1:nFsh){
        cat('f =',f,'\n')
        dm<-reshape2::melt(res$fisheries[[f]]$dm$n.xy,value.name='mortality')
        dm<-reshape2::dcast(dm,sex+year~.,value.var='mortality',fun.aggregate=sum);
        names(dm)[3]<-'mortality';
        dm$explrate<-dm$mortality/n$N;
        dm$fishery<-res$mc$lbls.fsh[f];
        dm$type<-'discards mortality';
        dfr<-rbind(dfr,dm);
        if (!is.null(res$fisheries[[f]]$rm)){
            rm<-reshape2::melt(res$fisheries[[f]]$rm$n.xy,value.name='mortality')
            rm<-reshape2::dcast(rm,sex+year~.,value.var='mortality',fun.aggregate=sum);
            names(rm)[3]<-'mortality';
            rm$explrate<-rm$mortality/n$N;
            rm$fishery<-res$mc$lbls.fsh[f];
            rm$type<-'retained mortality';
            dfr<-rbind(dfr,rm);
        }
    }
    dfr$year<-as.numeric(dfr$year);
    
    ps<-list();
    
    #plot exploitation rates for each sex
    usx<-tolower(unique(tm$sex))
    for (sx in usx){
        cat("plotting exploitation rates for",sx,'\n')
        tmp<-tm[tolower(tm$sex)==sx,]
        nmp<-nm[tolower(nm$sex)==sx,]
        dfp<-dfr[tolower(dfr$sex)==sx,]
        p <- ggplot(data=dfp)
        p <- p + geom_line(aes(x=year,y=explrate,colour=type),data=tmp,stat="identity",position='identity',alpha=1.0)
        p <- p + geom_line(aes(x=year,y=explrate,colour=type),data=nmp,stat="identity",position='identity',alpha=1.0)
        p <- p + geom_line(aes(x=year,y=explrate,colour=type),stat="identity",position='identity',alpha=1.0)
        p <- p + scale_x_continuous(breaks=pretty(dfp$year)) 
        p <- p + geom_hline(yintercept=0,colour='black',size=0.5)
        p <- p + labs(x="year",y="Exploitation Rate")
        p <- p + ggtitle(paste(sx,'s',sep=''))
        p <- p + guides(color=guide_legend(''))
        p <- p + facet_wrap(~fishery,ncol=2) 
        if (showPlot) print(p)
        ps[[sx]]<-p;
    }
    
    return(ps)
}

#ps<-plotExploitationRatesGG(res)