R/exploitationrates.r

In jae0/bio.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

    # ------------ fisheries data
    landings = bio.snowcrab::snowcrab_landings_db()
    # year is year of capture
    # yr is "fishing year" relative to the assessment cycle

    landings.north = landings[ which(landings$cfa =="cfanorth"), ]
    ln = aggregate(landings.north$landings, list(yr=landings.north$yr), function(x) sum(x, na.rm=T))
    names(ln) = c("yr", "landings")
    ln = factor2number(ln, c("yr", "landings"))
    ln$id = paste(ln$yr, "cfanorth", sep=".")

    landings.south = landings[ which(landings$cfa =="cfasouth") , ]
    ls = aggregate(landings.south$landings, list(yr=landings.south$yr), function(x) sum(x, na.rm=T))
    names(ls) = c("yr", "landings")
    ls = factor2number(ls, c("yr", "landings"))
    ls$id = paste(ls$yr, "cfasouth", sep=".")

    landings.4x = landings[ which(landings$cfa =="cfa4x" ) , ]
    lx = aggregate(landings.4x$landings, list(yr=landings.4x$yr), function(x) sum(x, na.rm=T))
    names(lx) = c("yr", "landings")
    lx = factor2number(lx, c("yr", "landings"))
    lx$id = paste(lx$yr, "cfa4x", sep=".")

    # lx$landings[2:(nrow(lx))] = lx$landings[1:(nrow(lx)-1)]  # offset due to the fishing season being after the Survey

    la = rbind(ln, ls, lx)
    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)

  }