R/exploitationrates.r

In jae0/snowcrab: Snow crab assessment suite using aegis* and associated projects and data streams.

  exploitationrates = function( p, areas, labels=NULL, CFA4X.exclude.year.assessment=FALSE  ){

    if ( is.null(labels)) labels = areas

    # regions = list( region=c("cfanorth", "cfasouth", "cfa4x"))
    regions = list(subarea=c("cfanorth", "cfa23", "cfa24", "cfa4x"))

    la = fishery_data( toget="summary_annual", regions=regions )

    la$yr = NULL
    la$landings.kt = as.numeric(as.character(la$landings))/1000/1000
 
    # -----------  biomass data
    p$vars.to.model ="R0.mass"

    p = parallel_run( p=p, runindex=list(y=p$yrs, v=p$vars.to.model ) )

    td = interpolation.db( DS="interpolation.simulation", p=p )
    td = td[ which(td$region %in% areas ), ]
    td$id = paste(td$yr, td$region, sep=".")
    td$biomass = td$total / 1000
    td$lbound =  td$total.ub / 1000
    td$ubound =  td$total.lb / 1000

    td$total = NA

    X = merge( td, la, by="id", all.x=T, all.y=F, sort=F)
    td = X[ which(is.finite( X$biomass)) ,]

    years= sort(unique(td$yr))

    for (rs in areas) {
    for (ys in years) {

      i = which(td$yr == ys & as.character(td$region) == rs )
      j = which(td$yr == (ys-1) & as.character(td$region) == rs )

      # post-2002 surveys have surveys in autumn
        if (ys >= 2002) {
          td$total[i] = td$landings.kt[i] / (td$biomass[i] + td$landings.kt[i])
          td$lb[i] = td$landings.kt[i] / (td$ubound[i] + td$landings.kt[i])
          td$ub[i] = td$landings.kt[i] / (td$lbound[i] + td$landings.kt[i])
        } else if (ys <= 2001 ) {
          td$total[i] = td$landings.kt[i] / td$biomass[i]
          td$lb[i] = td$landings.kt[i] / td$ubound[i]
          td$ub[i] = td$landings.kt[i] / td$lbound[i]
        }
    }}


    td$region = factor(as.character(td$region), levels=areas, labels=labels)

    eps = 1e-6

    varstocheck = c("total", "ubound", "lbound")

    for (vs in varstocheck) {
      td[,vs] = td[,vs]
      kk = which(td[,vs] <= eps)
      if (length(kk)>0) td[kk,vs] = 0
    }

    td = td[is.finite(td$total) ,]
    td = td[order(td$region,td$yr),]

    if (CFA4X.exclude.year.assessment) {
        td$total[ which( td$region=="4X" & td$yr==p$year.assessment) ] = NA  # the fishery is ongoing
        td$ub[ which(td$region=="4X" & td$yr==p$year.assessment) ] = NA  # the fishery is ongoing
        td$lb[ which(td$region=="4X" & td$yr==p$year.assessment) ] = NA  # the fishery is ongoing
    }

    ii = which(td$region=="4X" & td$yr<2004)
    td$total[ ii ] = NA
    td$ub[ ii ] = NA
    td$lb[ ii ] = NA

    return(td)

  }