R/aspic-plot.R

In laurieKell/mpb: Management Procedures

utils::globalVariables(c("ggplot","geom_line","aes","yield",
                         "geom_point","cast","xlab","ylab"))
utils::globalVariables(c("scale_alpha_manual"))
utils::globalVariables('data')


setMethod('plot', signature(x='aspics',y='missing'),
          function(x, y, probs=c(0.95,0.75,0.50,0.25,0.05),type=7,na.rm=FALSE,
                   facet=facet_wrap(~qname,scales='free',ncol=1),
                   fn=list('Stock'  =function(x) stock(x), 
                           'Harvest'=function(x) harvest(x),
                           'Yield'  =function(x) catch(x)),...)
          {  
            res=ldply(x,function(x) plot(x,probs=probs,type=type,na.rm=na.rm)$data)
            mdn=ldply(x,function(x) whooow(x,fn,probs=.5))
            
            if (length(unique(res$CI))>=1){
              p=ggplot(res)+
                geom_ribbon(aes(year,ymax=ymax,ymin=ymin,alpha=CI,group=paste(.id,CI),fill=.id))+
                geom_line(  aes(year,data,col=.id,group=.id),data=mdn)+
                scale_alpha_manual(values=seq(.1,.75,length.out=length(unique(res$CI))))
              
            }else{
              p=ggplot(res)+
                geom_line(aes(year,data,col=.id,group=.id),data=mdn)
              
            }
            
            p = p + facet + 
              expand_limits(y = 0) +
              xlab('Year') + ylab('')
            
            p})

#' plot aspics,missing
#' plotIndex aspic
#' plotIndex aspics

# @param  \code{fn}, a list of functions that estimate the quantities for plotting
# @param  \code{probs}, a vector specifying the percentiles for plotting, these are c(0.95,0.50,0.05) by default.
# @param  \code{size}, thinkness of percentile lines
# @param  \code{lty}, line type for percentiles
# @param \code{facet}, a layer that determines the facetting of the plot


plotProduction=function(object,biomass=FLQuant(seq(0,max(params(object)["k"]),length.out=101))) {
  object=as(object,"biodyn")
  if ((dims(object)$iter>1 | dims(params(object))$iter>1) & dims(biomass)$iter==1) 
    biomass=propagate(biomass,max(dims(object)$iter,dims(params(object))$iter))
  
  p=model.frame(FLQuants(stock=biomass, yield=FLQuant(production(object,biomass))))
  p=ggplot(p)+
    geom_line(aes(stock, yield, group=iter, col=iter)) +
    geom_point(aes(bmsy,msy,col=iter),size=4,
               data=model.frame(refpts(object))) +
    geom_path(aes(stock, yield, group=iter, col=iter),
        data=model.frame(mcf(FLQuants(
                 stock=stock(object),
                 yield=catch(object))))) +
    xlab("Stock") + ylab("Surplus Production")

  p}