In iantaylor-NOAA/Lingcod_2021: Estimate the Status of Lingcod off the US West Coast in 2021

###############################################################################
# Working directory: in doc, e.g., lingcod_2021\doc
# rmarkdown::render("catch.Rmd")
# Results are saved in ..\data-raw\catch.RData
#
# TO DO
# 1. Add source to the commercial catch stream
# 2. standardize catch versus landing
# 3. Assign fleet numbers and fleet names
# 4. write a function to format data_catch into SS.dat format
# 5. create a figure showing similarities of composition data from OR and
# northern CA recreational landings
# 7. Think about splitting catch_comm_CA_CALCOM ERK into area
# where right now it is all going to the North
# 9. for catch_rec_2000 maybe assign only a portion of Redwood
# from early years to the north model
# 9. Add variable at top for scrollable tables to either be included or omitted
# and create a table for every data set. these can then be included if we want
# or omitted when making the final .md file that we will source.
#
# PIE IN THE SKY
# 1. Why are WDFW 1995 to 1981 landings different than those in PacFIN
# MH (2017) noted that no one at WDFW knew why and that Phil was going to investigate?
# 2. Reconstructions using Sette and Fiedler for other species,
# JF noted data-based reconstructions may be better than a linear ramp
# 3. Recreational catch: MODE, how party boat vs. private boat landings differ
# 4. inquire to WDFW if reconstruction includes tribal catches via
# https://www.pcouncil.org/documents/2017/03/i2_att1_catch_reconstruction_workshop_report_mar2017bb.pdf/
# 5. Create a file name syntax with attributes to store who provided the file
# and perhaps other metadata that would be useful for latex text in the document,
# e.g., citation, date provided.
###############################################################################

###############################################################################
# Set up objects, file paths,
###############################################################################

# Input required by user to change these things
calcom_fleet_TW <- c("TWL", "NET")
erddapmaxyrmeanbyarea <- 1933
catch_comm_WA_yearstartpacfin <- 1995

# patterns for dir(file.path("..", "data-raw"), pattern = grep_...)
grep_pacfin <- "PacFIN.+FT.+RData"
grep_recweightfile <- "rec_bio_lingcod"
grep_comm_or <- "FINAL COMMERCIAL LANDINGS"
grep_comm_ca_or <- "CAlandingsCaughtORWA"
# grep_comm_ca_ralston provided by E.J. Dick
grep_comm_ca_ralston <- "reconstruction_Ralston_et_al_2010"
grep_comm_ca_calcom <- "1969-1980_CALCOM"
grep_rec_501 <- "CTE501"
grep_rec_mrfss <- "MRFSS-CATCH-ESTIMATES"
grep_rec_or <- "FINAL RECREATIONAL LANDINGS"
grep_rec_wa <- "Lingcod_RecCatch"

# paths for file.path("..", "data-raw", file_...)
file_comm_wa <- "WA_historical_Lingcod_forMelissa_V2.csv"
file_albin <- "Albin_et_al_1993_Lingcod_rows.csv"
file_erddap <- dir(
  path = file.path("..", "data-raw"),
  pattern = "earlyCAlingcod",
  full.names = TRUE,
  recursive = TRUE
)
# file_splitralston provided by MM
file_splitralston <- "Historical commercial catch for Kelli.xlsx"
file_password <- file.path(
    dirname(system.file(package = "PacFIN.Utilities")),
    "password.txt"
  )

###############################################################################
# Commercial data
###############################################################################

# Read in PacFIN file
if (file.exists(file_password)) {
  password <- readLines(file_password)
  # Download the data
  catch.pacfin <- PullCatch.PacFIN("LCOD",
    password = password, verbose = FALSE,
    addnominal = TRUE)
  rm(password)
}
load(dir_recent(file.path("..", "data-raw"), pattern = grep_pacfin))
catch.pacfin <- PacFIN.Utilities::getGearGroup(catch.pacfin, verbose = FALSE)
catch.pacfin <- catch.pacfin %>%
  dplyr::rename(Year = LANDING_YEAR) %>%
  dplyr::filter(
    !(AGENCY_CODE == "W" & Year %in% 0:(catch_comm_WA_yearstartpacfin - 1)),
    !(AGENCY_CODE == "O" & Year <
      # Read in this file from OR to determine which years they provide
      # that should be removed from the PacFIN data
      readxl::read_excel(
        path = dir_recent(
          dir = file.path("..", "data-raw"),
          pattern = grep_comm_or,
        ),
        sheet = 2
      ) %>% dplyr::filter(SOURCE != "PacFIN") %>% dplyr::pull(YEAR) %>% max
    )
  ) %>%
  dplyr::mutate(
    PACFIN_GEAR_CODE = dplyr::case_when(
      GEAR_NAME == "SCALLOP DREDGE" ~ "SCD",
      TRUE ~ PACFIN_GEAR_CODE
    ),
    fleet = use_fleetabb(geargroup),
    eez = dplyr::case_when(
      INPFC_AREA_TYPE_CODE == "XX" ~ "outside",
      ORIG_PACFIN_CATCH_AREA_CODE %in% c("02", "4A") ~ "outside",
      TRUE ~ "inside"
    )
  ) %>%
  dplyr::mutate(eezcatch = sum(ifelse(
    eez == "inside", 0, ROUND_WEIGHT_MTONS),
    na.rm = TRUE)
  ) %>%
  dplyr::filter(eez == "inside")
removenonEEZcatch <- catch.pacfin[1,"eezcatch"]

# Washington commercial data
catch_comm_WA <- utils::read.csv(file = file.path("..", "data-raw", file_comm_wa)) %>%
  dplyr::filter(Year < catch_comm_WA_yearstartpacfin) %>%
  dplyr::mutate(
    LINE = ifelse(LINE == 0, NA_real_, LINE),
    LINE = stats::approx(Year, LINE, n = length(LINE))$y,
    FG = LINE, TW = TRAWL + NET
  ) %>%
  dplyr::select(-Total, -LINE, -NET, -TRAWL) %>%
  tidyr::gather(key = "fleet", value = "mt", -Year) %>%
  dplyr::mutate(area = "North", state = "WA")

# Oregon commercial data
# pacfin.psmfc.org/pacfin_pub/data_rpts_pub/code_lists/agency_gears.txt
# pacfin.psmfc.org/pacfin_pub/data_rpts_pub/code_lists/gr_tree.txt
# 140 == Gill net; 210 == Gill net; 440 == Bait shrimp pump; 490 == Other H&L
catch_comm_OR <- readxl::read_excel(
  path = dir(path = file.path("..", "data-raw"), pattern = grep_comm_or, full.names = TRUE),
  sheet = 2
) %>%
  dplyr::rename(Year = YEAR) %>%
  dplyr::filter(SOURCE != "PacFIN") %>%
  dplyr::mutate(
    area = "North",
    state = "OR",
    mt = Total,
    PACFIN_GEAR_CODE = catch.pacfin[
      match(GEAR_CODE, catch.pacfin[, "GEAR_CODE"]),
      "PACFIN_GEAR_CODE"],
    PACFIN_GEAR_CODE = dplyr::case_when(
      GEAR_CODE %in% c(330) ~ "ONT", # squid nets
      GEAR_CODE %in% c(440, 490) ~ "OHL",
      GEAR_CODE %in% c(470) ~ "OTH", # Hand device
      GEAR_CODE %in% c(140, 210) ~ "GLN",
      GEAR_CODE %in% c(999) ~ "BTT",
      TRUE ~ PACFIN_GEAR_CODE
    )
  ) %>%
  data.frame
catch_comm_OR <- PacFIN.Utilities::getGearGroup(catch_comm_OR, verbose = FALSE)
catch_comm_OR[, "fleet"] <- use_fleetabb(catch_comm_OR[["geargroup"]])
catch_comm_OR_summary <- catch_comm_OR %>%
  dplyr::group_by(Year) %>%
  dplyr::summarize(sum = sum(mt), .groups = "keep") %>%
  dplyr::ungroup() %>%
  dplyr::mutate(mean = mean(sum)) %>%
  dplyr::slice_max(order_by = sum, n = 1)

# California commercial data
# erddap
catch_erddap <- readxl::read_excel(
  path = file_erddap,
  sheet = 3
) %>%
  dplyr::filter(!is.na(fish), port != "All") %>%
  dplyr::rename(Year = year) %>%
  dplyr::rename_with(~ gsub("landings|\\.|\\(|\\)", "", .x)) %>%
  dplyr::group_by(Year) %>%
  dplyr::mutate(mt_y = sum(mt)) %>%
  dplyr::group_by(port, Year) %>%
  dplyr::mutate(mt_yp = sum(mt)) %>%
  dplyr::ungroup()
erddapmeanpropEureka <- catch_erddap %>%
  dplyr::mutate(prop_p = mt_yp/mt_y) %>%
  dplyr::filter(port == "Eureka", Year < erddapmaxyrmeanbyarea) %>%
  dplyr::pull(prop_p) %>% mean
erddapmeanpropNEureka <- readxl::read_excel(
  path = file_erddap,
  sheet = 1, n_max = 1,
  skip = 30, .name_repair = "minimal",
) %>% dplyr::pull(8)

# sette
catch_sette_1926fa <- data.frame(
  district = rep(
    c("Northern", "San Francisco", "Monterey", "Southern"),
    times = c(3, 4, 4, 1)
  ),
  gear = c(
    "lines", "gill nets", "paranzella nets",
    "lines", "gill nets", "lampara nets", "paranzella nets",
    "lines", "paranzella nets", "gill nets", "lampara nets",
    "trammel nets"
  ),
  lbs = c(
    35107, 0, 13449,
    309151, 0, 0, 140363,
    98383, 47265, 0, 34,
    1248
  )
) %>% dplyr::mutate(
  fleet = dplyr::case_when(
    gear %in% c("gill nets", "paranzella nets", "lampara nets", "trammel nets") ~ "TW",
    TRUE ~ "FG"
  ),
  mt = lbs * mult_lbs2mt()
) %>%
  dplyr::mutate(
    area = ifelse(district == "Northern", "North", "South")
  ) %>%
  dplyr::group_by(fleet, area) %>%
  dplyr::mutate(mt_fa = sum(mt)) %>%
  dplyr::group_by(area) %>%
  dplyr::mutate(prop_fa = mt_fa / sum(mt))
# Calculate the proportion for fleet and area from 1926 to use for
# multiple data sets
catch_sette_prop_fa <- catch_sette_1926fa %>% dplyr::distinct(fleet, area, prop_fa)
# Sette and Fiedler: transcribed by KFJ
catch_sette1928 <- data.frame(
  Year = c(
    1888, 1892, 1895, 1899, 1904, 1908, 1915, 1918, # pg 528
    1919:1926 # pg 529
    ),
  lbs = c(
    0, 231.0, 139.0, 148.0, 293.0, 167.0, 578.0, 916.0, # pg 528
    1063.0, 688.0, 426.0, 568.0, 467.0, 400.0, 683.0, 645.0 # pg 529
    ) * 1000
) %>%
  tidyr::crossing(., catch_sette_prop_fa) %>%
  dplyr::mutate(
    prop_a = dplyr::case_when(
      area == "North" ~ (erddapmeanpropNEureka * erddapmeanpropEureka),
      TRUE ~ 1 - (erddapmeanpropNEureka * erddapmeanpropEureka)
    ),
    mt = lbs * mult_lbs2mt() * prop_a * prop_fa,
    source = "sette"
  )

# CA landings in OR waters
# 1948 - 1968 additional catch in OR waters landed in CA from John F.
catch_comm_CA_ORwaters <- readxl::read_excel(
  path = dir(path = file.path("..", "data-raw"), pattern = grep_comm_ca_or, full.names = TRUE),
  sheet = 3,
  skip = 5,
  .name_repair = "minimal"
) %>% dplyr::select(1, 8, 9) %>%
  dplyr::rename(Year = 1) %>%
  dplyr::mutate(
    area = "North",
    fleet = "TW",
    mt = (Oregon + Washington),
    source = "Field"
  )

# Ralston
catch_comm_CA_ralston <- utils::read.csv(
  file = dir(
    file.path("..", "data-raw"),
    pattern = grep_comm_ca_ralston,
    full.names = TRUE
  )
) %>%
  dplyr::rename(Year = 1) %>%
  dplyr::mutate(species = gsub("\\s", "", species)) %>%
  dplyr::mutate(
    area = dplyr::case_when(
      region_caught == 2 ~ "North_South",
      region_caught %in% c(1, 11, 12) ~ "North",
      TRUE ~ "South"
    )
  ) %>% tidyr::separate_rows(area, sep = "_") %>%
  dplyr::mutate(
    mt = dplyr::case_when(
      region_caught == 2 & area == "North" ~ pounds * mult_lbs2mt() * erddapmeanpropNEureka,
      region_caught == 2 & area == "South" ~ pounds * mult_lbs2mt() * (1-erddapmeanpropNEureka),
      TRUE ~ pounds * mult_lbs2mt()
    ),
    state = dplyr::case_when(
      region_caught %in% 11 ~ "OR",
      region_caught %in% 12 ~ "WA",  #maybe, but I know it is N,
      TRUE ~ "CA"
    ),
    source = "ralston"
  ) %>%
  dplyr::filter(state == "CA")

# CALCOM
catch_comm_CA_calcom <- utils::read.csv(
  file = dir(
    file.path("..", "data-raw"),
    pattern = grep_comm_ca_calcom,
    full.names = TRUE
  )
) %>% dplyr::rename(Year = 1)
catch_comm_CA_MM <- readxl::read_excel(
  path = file.path("..", "data-raw", file_splitralston),
  .name_repair = "minimal",
  skip = 1
) %>%
  dplyr::rename(Year = 1, FG = 2, TW = 3) %>%
  dplyr::mutate(
    area = ifelse(grepl("[NS]", Year), Year, NA_character_)
  ) %>%
  tidyr::fill(area, .direction = "down") %>%
  dplyr::filter(
    !grepl("cape|North|South|Year", Year),
    !is.na(Year)
  ) %>%
  dplyr::mutate_at(c("Year","FG","TW"), as.numeric) %>%
  tidyr::gather(key = "type", value = "prop_fa", -Year, -area) %>%
  dplyr::mutate(
    area = dplyr::case_when(
      is.na(area) & type == "FG" ~ "North_all",
      is.na(area) & type == "TW" ~ "South_all",
      TRUE ~ area
    ),
    type = ifelse(grepl("all", area), "ALL", type),
    area = gsub("_all", "", area)
  ) %>%
  data.frame

# Find landings in Eureka
# Sort commercial catches out to area
catch_comm_CA_ERA <- catch.pacfin %>%
  dplyr::filter(PACFIN_GROUP_PORT_CODE %in% c("ERA")) %>%
  dplyr::group_by(PORT_NAME) %>%
  dplyr::summarize(.groups = "keep",
    n = length(unique(VESSEL_ID)),
    mt = sum(ROUND_WEIGHT_MTONS)
  ) %>%
  dplyr::arrange(desc(mt)) %>%
  dplyr::filter(PORT_NAME != "SHELTER COVE") %>%
  dplyr::mutate(
    out = ifelse(n < 4, "", paste0(" (", round(mt, 2), " mt)"))
  ) %>%
  dplyr::mutate(PORT_NAME = paste0(stringr::str_to_title(PORT_NAME, locale = "en"), out)) %>%
  dplyr::pull(PORT_NAME)
#
# Find vessels in CA with unassigned port
# Determine vessels fished by visually inspecting
# all ports with lingcod landings by vessel_id
catch_comm_CA_noport <- catch.pacfin[
  catch.pacfin[["AGENCY_CODE"]] == "C" &
  catch.pacfin[["PORT_CODE"]] %in% c(0, 896),
  ] %>% dplyr::group_by(VESSEL_ID, Year) %>%
  dplyr::summarize(sumMT = sum(ROUND_WEIGHT_MTONS), .groups = "keep")
catch_comm_CA_noport_vesselid <- catch.pacfin %>%
  dplyr::filter(VESSEL_ID %in% catch_comm_CA_noport[["VESSEL_ID"]]) %>%
  dplyr::group_by(VESSEL_ID, PACFIN_PORT_NAME) %>%
  dplyr::summarize(mt = sum(ROUND_WEIGHT_MTONS), .groups = "keep") %>%
  dplyr::mutate(north = dplyr::case_when(
      PACFIN_PORT_NAME %in% c("CRESCENT", "EUREKA") ~ sum(mt),
      TRUE ~ 0
      )
  ) %>%
  dplyr::ungroup() %>% dplyr::group_by(VESSEL_ID) %>%
  dplyr::summarize(
    ratio = sum(north) / sum(mt)
  ) %>% dplyr::arrange(ratio, VESSEL_ID)

# Summarize commercial pacfin catches
catch_comm_pacfin <- catch.pacfin %>%
  dplyr::mutate(
    area = dplyr::case_when(
      AGENCY_CODE %in% c("W", "O") ~ "North",
      PORT_NAME %in% "SHELTER COVE" ~ "South",
      PORT_NAME %in% "BAYSIDE" ~ "North",
      PACFIN_GROUP_PORT_CODE %in% "ERA" ~ "North",
      PACFIN_GROUP_PORT_CODE %in% "CCA" ~ "North",
      PORT_CODE %in% c(0, 896) &
      VESSEL_ID %in% catch_comm_CA_noport_vesselid[catch_comm_CA_noport_vesselid[["ratio"]] > 0.5, "VESSEL_ID", drop = TRUE] ~ "North",
      TRUE ~ "South"
    ),
  mt = ROUND_WEIGHT_MTONS,
  state = PacFIN.Utilities:::getState(catch.pacfin)[["state"]]
  )
catch_comm_pacfin_summary <- catch_comm_pacfin %>%
dplyr::filter(!(state == "OR" & Year <= max(catch_comm_OR[["Year"]]))) %>%
dplyr::group_by(Year, state, fleet) %>%
dplyr::summarize(sum = sum(mt), .groups = "keep") %>%
dplyr::group_by(Year, state) %>%
dplyr::mutate(
  prop = sum / sum(sum),
  yearsum = sum(sum)
) %>%
dplyr::group_by(state) %>%
dplyr::mutate(
  max = sum[which.max(sum)],
  min = min(sum),
  maxYear = Year[which.max(sum)],
  mean = mean(prop)
)

# Fill in missing years for California commercial catches
catch_comm_CA_interpolate <- dplyr::full_join(
    by = c("Year", "area", "fleet", "source", "mt"),
    # Build up to 1930
    x = dplyr::full_join(
      by = c("Year", "source", "fleet", "area", "mt", "prop_fa"),
      x = catch_sette1928,
      y = dplyr::left_join(
          by = "area",
          x = catch_erddap %>%
            dplyr::group_by(Year, port) %>%
            dplyr::slice(1L) %>%
            dplyr::mutate(
              North = dplyr::case_when(
                port == "Eureka" ~ mt_yp * erddapmeanpropNEureka,
                TRUE ~ 0
                ),
              South = dplyr::case_when(
                port == "Eureka" ~ mt_yp * (1-erddapmeanpropNEureka),
                TRUE ~ mt_yp
              )
            ) %>%
            dplyr::select(-time, -fish, -(lbs:mt_yp)) %>%
            tidyr::gather(key = area, value = mt, -Year, -port) %>%
            dplyr::filter(
              mt != 0,
              Year < min(catch_comm_CA_ralston[["Year"]])
            ) %>%
            dplyr::mutate(source = "erddap"),
          y = catch_sette_prop_fa
        ) %>%
        dplyr::mutate(mt = mt * prop_fa) %>%
        dplyr::group_by_at(vars(-mt, -port)) %>%
        dplyr::summarize(mt = sum(mt), .groups = "keep") %>%
        dplyr::ungroup()
    ),
    # Build since 1931
    y = dplyr::full_join(
      by = c("Year", "area", "fleet", "mt", "source"),
      x = catch_comm_CA_ralston %>%
        dplyr::group_by(Year, area) %>%
        dplyr::summarize(mt = sum(mt), .groups = "keep") %>%
        dplyr::ungroup()  %>%
        tidyr::crossing(data.frame(fleet = c("FG", "TW"))) %>%
        dplyr::left_join(
          by = c("Year", "area", "fleet"),
          x = .,
          y = dplyr::full_join(
            # Create time series of prop with MM and erddap
            by = c("Year", "fleet", "area", "prop_fa"),
            x = catch_sette_prop_fa %>%
              dplyr::mutate(Year = min(catch_comm_CA_ralston[["Year"]]) - 1),
            y = catch_comm_CA_MM %>%
              dplyr::filter(type != "ALL") %>%
              dplyr::rename(fleet = "type")
          ) %>%
          dplyr::group_by(fleet, area) %>%
          tidyr::complete(fleet, area, Year = min(Year):max(Year)) %>%
          dplyr::mutate(prop_fa = stats::approx(Year, prop_fa, n = length(Year))$y) %>%
          dplyr::ungroup()
        ) %>%
        dplyr::mutate(mt = mt * prop_fa, source = "ralston"),
      y = catch_comm_CA_calcom %>%
        dplyr::mutate(
          source = "CALCOM",
          fleet = ifelse(GEAR_GRP %in% calcom_fleet_TW, "TW", "FG"),
          area = ifelse(PORT_COMPLEX %in% c("CRS", "ERK"), "North", "South") # might need to split ERK
        ) %>%
        dplyr::group_by(Year, fleet, area, source) %>%
        dplyr::summarize(mt = sum(POUNDS) * mult_lbs2mt(), .groups = "keep") %>%
        dplyr::ungroup()
      ) %>% dplyr::full_join(
          by = c("Year", "source", "mt", "fleet", "area"),
          x = .,
          y = catch_comm_CA_ORwaters
        )
  ) %>%
  # Work with the combined early and late data set
  tidyr::complete(area, fleet,
    Year = min(Year):max(Year),
    fill = list(source = "interpolate")
  ) %>%
  dplyr::group_by(Year, area, fleet) %>%
  dplyr::summarize(.groups = "keep", mt = sum(mt)) %>%
  dplyr::ungroup(Year) %>%
  dplyr::mutate(
    mt_orig = mt,
    mt = stats::approx(Year, mt, n = length(Year))$y,
    state = "CA"
  ) %>%
  data.frame

myby <- c("Year", "fleet", "area", "mt", "state")
catch_comm_reconstruction_notsummed <- dplyr::full_join(
  by = myby,
  x = catch_comm_WA,
  y = catch_comm_OR) %>%
  dplyr::full_join(
    by = myby,
    x = .,
    y = catch_comm_CA_interpolate
  )
catch_comm_reconstruction <- catch_comm_reconstruction_notsummed %>%
  dplyr::group_by(Year, area, fleet) %>%
  dplyr::summarize(mt = sum(mt), .groups = "keep") %>%
  dplyr::ungroup()

###############################################################################
# Summarize the commercial catches three ways
# this information is helpful for future runs and for building the document
# 1. For the models: year, area, fleet
# 2. For the document: year, area, state, fleet
# 3. For comparison with the 2017 assessment boundaries at state lines: year, state, fleet
###############################################################################

#### # 1
catch_comm <- dplyr::full_join(
  by = c("Year", "area", "fleet", "mt"),
  x = catch_comm_reconstruction,
  y = catch_comm_pacfin
) %>% dplyr::group_by(Year, area, fleet) %>%
  dplyr::summarize(mt = sum(mt), .groups = "keep") %>%
  dplyr::ungroup() %>%
  data.frame

#### # 2 Table 1. Catches by state x fleet
catch_comm_fordoc <- dplyr::full_join(
    by = c("Year", "area", "state", "fleet", "mt"),
    x = catch_comm_pacfin,
    y = catch_comm_reconstruction_notsummed
  ) %>%
  dplyr::group_by(Year, area, state, fleet) %>%
  dplyr::summarize(mt = sum(mt), .groups = "keep") %>%
  dplyr::ungroup()

#### # 3 Run the model with the previous stratification at state lines
####   combine with the old data file catches
catch_comm_state <- dplyr::full_join(
    by = c("Year", "state", "fleet", "mt"),
    x = catch_comm_pacfin,
    y = catch_comm_reconstruction_notsummed
  ) %>%
  dplyr::mutate(
    state = ifelse(state != "CA", "N", state),
    source = "2021"
  ) %>%
  dplyr::group_by(Year, state, fleet, source) %>%
  dplyr::summarize(mt = sum(mt), .groups = "keep") %>%
  dplyr::ungroup() %>%
  dplyr::full_join(
    by = c("Year", "fleet", "state", "source"),
    x = .,
    y = dplyr::bind_rows(.id = "id",
  lapply(list(last_ssdat[[1]], last_ssdat[[2]]),
    "[[",
    "catch"
  )
) %>%
  dplyr::rename(Year = "year") %>%
  dplyr::mutate(
    source = "2019",
    fleet = dplyr::case_when(
      id == 1 ~ last_ssdat[[1]][["fleetnames"]][fleet],
      id == 2 ~ last_ssdat[[2]][["fleetnames"]][fleet],
    )
  ) %>%
  tidyr::separate(
    col = "fleet",
    into = c("fleetn", "state", "fleet"),
    sep = "_") %>%
  dplyr::filter(Year > 0, !grepl("REC", fleet)) %>%
  dplyr::mutate(fleet = ifelse(fleet == "FIX", "FG", "TW"))
  ) %>%
  dplyr::mutate(
    state = factor(state, levels = c("N", "CA"), labels = c("WA and OR", "CA"))
  )
# Perform a test to make sure that sum by state is same as sum by gear
testthat::expect_equal(dplyr::full_join(by = "Year",
  catch_comm_state %>%
    dplyr::group_by(Year) %>%
    dplyr::summarize(sum = sum(mt, na.rm = TRUE)),
  catch_comm %>%
    dplyr::group_by(Year) %>%
    dplyr::summarize(sum = sum(mt, na.rm = TRUE))
  ) %>% 
    dplyr::mutate(diff = sum.x-sum.y) %>%
      dplyr::pull(diff) %>% sum,
  0
)
# Save a file to re-run the last model with updated catches
catch_writetofile <- catch_comm_state %>%
  dplyr::filter(source == 2021) %>%
  dplyr::select(Year, state, fleet, mt) %>%
  dplyr::rename(year = Year, catch = "mt") %>%
  dplyr::mutate(
    seas = 1,
    fleet = dplyr::case_when(
      state == "CA" & fleet == "TW" ~ 1,
      state == "CA" & fleet == "FG" ~ 2,
      state == "WA and OR" & fleet == "TW" ~ 1,
      state == "WA and OR" & fleet == "FG" ~ 2
    ),
    catch_se = 0.00
  )
utils::write.csv(
  x = catch_writetofile %>% dplyr::filter(state == "CA"),
  file = file.path("..", "data-raw", "catch_comm_south_2019structure.csv")
)
utils::write.csv(
  x = catch_writetofile %>% dplyr::filter(state == "WA and OR"),
  file = file.path("..", "data-raw", "catch_comm_north_2019structure.csv")
)

###############################################################################
# Recreational data
###############################################################################
catch_rec_2000 <- RecFIN::read_cte501(
  file = dir(file.path("..", "data-raw"), pattern = grep_rec_501, full.names = TRUE, recursive = TRUE)
) %>%
  dplyr::mutate(Source = "RecFIN") %>%
  dplyr::filter(!grepl("CAN|MEX|PUG", WATER_AREA_NAME)) %>%
  dplyr::mutate(
    area = dplyr::case_when(
      grepl("REDWOOD", PORT_NAME) ~ "North",
      state == "CA" ~ "South",
      TRUE ~ "North"
      ),
    source = "CRFS"
  )

catch_rec_1980 <- RecFIN::read_mrfss(
  file = dir(file.path("..", "data-raw"), pattern = grep_rec_mrfss, full.names = TRUE, recursive = TRUE)
) %>%
  dplyr::filter(
    SPECIES_NAME == lingcod::utils_name(type = "Common"),
    Year > min(Year)
  ) %>%
  dplyr::mutate(Source = "RecFIN")

bio_rec_recfin <- utils::read.csv(
  file = dir_recent(file.path("..", "data-raw"), pattern = grep_recweightfile),
  header = TRUE
) %>% dplyr::rename(Year = "RECFIN_YEAR") %>%
  dplyr::mutate(
    AGENCY_WEIGHT_KG = dplyr::case_when(
      grepl("^G", AGENCY_WEIGHT_UNITS) ~ AGENCY_WEIGHT / 1000,
      TRUE ~ AGENCY_WEIGHT,
    ),
    sex = dplyr::case_when(
      FISH_SEX %in% c(2, "F") ~ "F",
      FISH_SEX %in% c(1, "M") ~ "M",
      FISH_SEX %in% c(3, 8, "U") ~ "U"
    )
  ) %>%
  # Remove bad sample for lingcod, already emailed Jason and Theresa T.
  dplyr::filter(!(RECFIN_DATE == "07/25/2007" & RECFIN_LENGTH_MM > 7499))

# Mean length to convert numbers to weight
bio_rec_recfin_meanlength <- bio_rec_recfin %>%
  dplyr::group_by(STATE_NAME, IS_RETAINED) %>%
  dplyr::summarize(
    mean_length = mean(RECFIN_LENGTH_MM, na.rm = TRUE),
    n = dplyr::n(),
    .groups = "keep"
  )

# Washington recreational data
catch_rec_WA <- suppressWarnings(suppressMessages(
  dplyr::bind_rows(.id = "Type", mapply(readxl::read_excel,
    MoreArgs = list(
      col_names = FALSE,
      col_types = "numeric",
      .name_repair = "minimal",
      path = dir(
        file.path("..", "data-raw"),
        pattern = grep_rec_wa,
        full.names = TRUE,
        recursive = TRUE
      )
    ),
    sheet = c(OSP = 1, PSSP = 2, historical = 3),
    skip = c(2, 1, 1),
    n_max = c(Inf, 15, Inf),
    SIMPLIFY = FALSE
  )) %>%
  dplyr::rename(Year = 2, LINGCOD = 3, RELEASED = 4)
) %>%
  dplyr::select(1:4)
) %>%
  dplyr::filter(!is.na(Year)) %>%
  dplyr::group_by(Year) %>%
  dplyr::summarize(.groups = "keep",
    num = sum(LINGCOD + ifelse(is.na(RELEASED), 0, RELEASED) * 0.07),
    numdeaddiscard = sum(ifelse(is.na(RELEASED), 0, RELEASED) * 0.07),
    numretained = sum(LINGCOD)
  ) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(interpolate = 0) %>%
  tidyr::complete(Year = tidyr::full_seq(Year, 1), fill = list(interpolate = 999)) %>%
  dplyr::mutate(
    num = stats::approx(Year, num, n = length(Year))$y,
    interpolate = dplyr::case_when(
      interpolate == 999 ~ num,
      TRUE ~ interpolate
    ),
    mean_length = bio_rec_recfin_meanlength %>% dplyr::filter(STATE_NAME == "WASHINGTON") %>% dplyr::pull(mean_length),
    mean_weight = lw.WCGBTS[["all"]][["a"]] * (mean_length / 10)^(lw.WCGBTS[["all"]][["b"]]),
    mt = mean_weight / 1000 * num,
    state = "WA",
    mtdeaddiscard = mean_weight / 1000 * numdeaddiscard,
    mtretained = mean_weight / 1000 * numretained,
  )

gg <- ggplot2::ggplot(catch_rec_WA, ggplot2::aes(Year, mtretained)) +
ggplot2::geom_line() +
ggplot2::geom_line(ggplot2::aes(y = mtdeaddiscard), col = "red") +
ggplot2::geom_line(ggplot2::aes(y = mt), col = "green") +
ggplot2::labs(ylab = "mt") +
ggplot2::theme_bw() + 
ggplot2::xlim(c(2000, 2020)) +
ggplot2::ylab("mt")
ggplot2::ggsave(gg, filename = file.path("..", "figures", "catch_rec_WA_retainedvdead.png"))

# Oregon recreational data
catch_rec_OR <- dplyr::full_join(by = c("Year", "mt"), .id = "dataset",
  # OR data from Ali
  readxl::read_excel(
    path = dir(
      file.path("..", "data-raw"),
      pattern = grep_rec_or,
      full.names = TRUE,
      recursive = TRUE
    ),
    sheet = 2
  ) %>%
    dplyr::rename(Year = YEAR) %>%
    dplyr::group_by(Year) %>%
    dplyr::mutate(mt = sum(`TOTAL MORTALITY_MT`)) %>%
    dplyr::select(!dplyr::matches("RET|REL")) %>%
    tidyr::spread(key = MODE, value = `TOTAL MORTALITY_MT`),
  # 2017 dat file from Northern model
  last_ssdat[[1]][["catch"]] %>%
    dplyr::filter(
      fleet == grep("OR", ignore.case = TRUE, last_ssdat[[1]][["fleetnames"]]),
      catch > 0
    ) %>%
    dplyr::rename(Year = year, mt = catch) %>%
    dplyr::select(Year, mt)
) %>%
  dplyr::distinct(Year, .keep_all = TRUE) %>%
  dplyr::arrange(Year) %>%
  dplyr::mutate(area = "North", state = "OR")

# California recreational data
# albin
# 2020-11-06 E.J. Dick, M. Monk, C. Wetzel, J. Budrick, I. Taylor
# engaged in an email conversation but it pertained to Point Reyes, i.e.,
# split between Mendocino/Sonoma and San Francisco, around 38.25
#
# 2021-04-14 E.J. Dick, M. Monk, I. Taylor, and K.F. Johnson
# engaged in a video chat regarding how to split at 40-10 because
# Redwood includes Humboldt County with Shelter Cove rather than
# Humboldt County without Shelter Cove, which is what we want.
# E.J. does not know of any other way to partition MRFSS data that has
# two groups, Northern and Southern, for SUB_REGION_NAME == CA
#   * 1990-1992 San Luis Obispo County in Southern region
#   * 1993-.... San Luis Obispo County in Northern region
#
# csv is group = "35. LINGCOD" of Table 1 in Albin et al. 1993 which
# contains county-specific catch estimates for 1981-1986
data_albin <- read.csv(file = file.path("..", "data-raw", file_albin),
  skip = 2, header = TRUE, check.names = FALSE
) %>%
  rlang::set_names(paste(sep = "_",
    read.csv(
      file = file.path("..", "data-raw", file_albin),
      skip = 1, header = FALSE, check.names = FALSE, nrows = 1
    ),
    colnames(.))
  ) %>%
  dplyr::rename(Year = "NA_Year") %>%
  tidyr::gather("type", "value", -Year) %>%
  tidyr::separate(type, into = c("Area", "type"), sep = "_") %>%
  tidyr::spread(key = "type", value = "value") %>%
  dplyr::arrange(Area) %>%
  dplyr::mutate(Source = "albinetal1993") %>%
  dplyr::filter(Area != "Total") %>%
  dplyr::group_by(Year) %>%
  dplyr::mutate(sum = sum(Est)) %>%
  dplyr::group_by(Area, Year) %>%
  dplyr::mutate(prop_source = Est / sum) %>%
  dplyr::ungroup()

albinmeanpropN <- data_albin %>%
          dplyr::group_by(Year) %>%
          dplyr::mutate(YearT = sum(Est)) %>%
          dplyr::ungroup() %>% dplyr::group_by(Area) %>%
          dplyr::summarize(wm = stats::weighted.mean(prop_source, w = YearT)) %>%
          dplyr::filter(grepl("Del", Area)) %>% dplyr::pull(wm)
catch_rec_CA <- catch_rec_1980 %>%
    dplyr::filter(state == "CA") %>%
    dplyr::group_by(Year, state) %>%
    dplyr::summarize(
      mt = sum(WGT_AB1, na.rm = TRUE) / 1000,
      .groups = "keep"
    ) %>%
    dplyr::ungroup() %>%
    dplyr::mutate(source = "MRFSS") %>%
  dplyr::full_join(
  by = c("Year", "state", "mt", "source"),
  x = .,
  y = dplyr::anti_join(
    by = c("Year", "mt", "state", "source"),
    y = .,
    x = last_ssdat[[2]][["catch"]] %>%
      dplyr::mutate(state = "CA", source = "old SS") %>%
      dplyr::filter(
        fleet == grep("CA_REC", last_ssdat[[2]][["fleetnames"]]),
        catch > 0,
        year < min(catch_rec_1980$Year)
      ) %>%
      dplyr::select(-seas, -catch_se, -fleet) %>%
      dplyr::rename(Year = year, mt = catch)
  )) %>%
    tidyr::crossing(area = get_groups(info_groups)) %>%
    dplyr::mutate(
      mt = dplyr::case_when(
        state == "CA" & area == "North" ~ mt * albinmeanpropN,
        state == "CA" & area == "South" ~ mt * (1 - albinmeanpropN),
        TRUE ~ mt
      )
    ) %>%
  dplyr::full_join(
    by = c("Year", "state", "mt", "source", "area"),
    x = .,
    y = catch_rec_2000 %>%
    dplyr::group_by(area, state, Year, source) %>%
    dplyr::summarize(mt = sum(TOTAL_MORTALITY_MT), .groups = "keep") %>%
    dplyr::filter(state == "CA") %>%
    dplyr::ungroup()
  ) %>%
  dplyr::group_by(area) %>%
  tidyr::complete(
    Year = tidyr::full_seq(Year, 1),
    fill = list(source = "interpolate", state = "CA")
  ) %>%
  dplyr::mutate(mt = stats::approx(Year, mt, n = length(Year))$y) %>%
  dplyr::ungroup()

###############################################################################
# Summarize recreational catches
# 1 For use in the stock assessment
# 2 For Table 1 in the document
###############################################################################

#### # 1
catch_rec <- dplyr::full_join(
  by = c("Year", "mt", "state"),
  catch_rec_WA,
  catch_rec_OR,
) %>%
  dplyr::select(Year, mt, state) %>%
  dplyr::mutate(area = "North") %>%
  dplyr::full_join(
    y = catch_rec_CA,
    by = c("Year", "area", "state", "mt"),
  )
# Save the file
data_catch <- dplyr::full_join(
  by = c("Year", "area", "mt", fleet = "state"),
  x = catch_comm,
  y = catch_rec %>% dplyr::select(-source)
)

# 2021-07-13 Day 2 STAR requests to fix a mistake in WA rec landings
# gg <- catch_rec %>%
#   dplyr::filter(area=="North") %>%
#   ggplot2::ggplot(ggplot2::aes(Year, mt,group=state,col=state)) +
#   ggplot2::geom_line(size = 1.1) +
#   ggplot2::theme_bw() +
#   ggplot2::theme(legend.position = c(0.1, 0.87)) +
#   ggplot2::scale_color_manual(values = c(rev(get_fleet("Rec", col = "col.n")), "black"))
# ggplot2::ggsave(plot = gg, filename = file.path("..", "figures", "catch_rec_WA_correcttotallandings.png"))

#### # 2
catch_fordoc <- dplyr::full_join(
  x = catch_comm_fordoc %>%
    dplyr::mutate(state = paste0(fleet, "-", state)) %>%
    dplyr::select(-fleet) %>%
    tidyr::spread("state", "mt"),
  y = catch_rec %>%
    # dplyr::mutate(state = paste0("rec-", state)) %>%
    dplyr::select(-source) %>%
    tidyr::spread("state", "mt"),
  by = c("Year", "area")
) %>%
  dplyr::mutate_if(is.numeric, tidyr::replace_na, 0) %>%
  dplyr::mutate(Total = purrr::reduce(dplyr::select(., !dplyr::matches("Year|area")), `+`)) %>%
  dplyr::rename_with(~ gsub("FG-", "  ", .x)) %>%
  dplyr::rename_with(~ gsub("TW-", " ", .x)) %>%
  dplyr::filter(Total != 0, Year < as.numeric(substr(Sys.Date(), 1, 4))) %>%
  dplyr::mutate_at(3:NCOL(.), ~ sprintf("%.2f", .x))

catchcaption <- "Commercial (by fleet and state) and recreational (by state) landings (mt) with yearly totals. Commercial landings were modeled as gear-specific landings and separated here by state for information purposes only."
catchlabel <- "table-catchbystate"
kableExtra::kbl(catch_fordoc %>% dplyr::filter(area == "North") %>% dplyr::select(-area),
  booktabs = TRUE, longtable = TRUE,
  caption = catchcaption, label = catchlabel, format = "latex",
  align = "r",
) %>%
kableExtra::kable_styling(latex_options = c("repeat_header")) %>%
kableExtra::add_header_above(c(" ", "Fixed gear" = 3, "Trawl" = 3, "Recreational" = 3, "")) %>%
kableExtra::save_kable(file = file.path("..", "tables", paste0(catchlabel, "-n-", ".tex")))
kableExtra::kbl(catch_fordoc %>% dplyr::filter(area == "South") %>% dplyr::select(-area),
  booktabs = TRUE, longtable = TRUE,
  caption = catchcaption, label = catchlabel, format = "latex",
  align = "r",
) %>%
kableExtra::add_header_above(c(" ", "Fixed gear" = 3, "Trawl" = 3, "Recreational" = 3, "")) %>%
kableExtra::kable_styling(latex_options = c("repeat_header")) %>%
kableExtra::save_kable(file = file.path("..", "tables", paste0(catchlabel, "-s-", ".tex")))

###############################################################################
# Save the data
###############################################################################
save.image(file = file.path("..", "data-raw", "catch.RData"))
usethis::use_data(data_catch, overwrite = TRUE)
usethis::use_data(catch_info, overwrite = TRUE)

###############################################################################
# Make the figures for the catches
###############################################################################
gg <- ggplot2::ggplot(
  catch_comm_state,
  ggplot2::aes(Year, mt, lty = source, pch = fleet, col = fleet)
) +
  ggplot2::geom_line() +
  # ggplot2::geom_point() +
  ggplot2::geom_line(ggplot2::aes(y = catch)) +
  # ggplot2::geom_point(aes(y = catch)) +
  ggplot2::facet_grid(state ~ .) +
  ggplot2::theme_bw() +
  ggplot2::theme(strip.background = ggplot2::element_rect(fill = NA)) +
  ggplot2::theme(legend.position = c(0.1, 0.8)) +
  ggplot2::ylab("Landings (mt)") +
  ggplot2::scale_color_manual(
    values = unlist(t(get_fleet(value="Comm")[,c("col.n")]),use.names=FALSE)
  )
ggsave2(
  gg,
  filename = file.path("..", "figures", "catch-comm-state.png"),
  caption = paste0(
    "Comparison of current landings (dashed line) and ",
    "those used in the previous assessment model (solid line) by state (panel). ",
    "Fixed gear (FG; dark blue) and trawl (TW; light blue) landings ",
    "are shown for Washington and Oregon (top panel) combined and ",
    "California (bottom panel) because ",
    "the previous assessment model used the Oregon-California border to define the stocks."
  ),
  alttext = paste(
    "Time series of commercial catches are less in the early years for the new southern model than the previous assessment."
  )
)

gg <- ggplot2::ggplot(catch_comm_CA_interpolate,
  ggplot2::aes(
    x = Year,
    y = mt_orig,
    col = interaction(area, fleet, sep = " - "),
  )
) +
  ggplot2::geom_line() +
  ggplot2::geom_line(ggplot2::aes(y = mt), lty = 2) +
  ggplot2::geom_point(cex = 2) +
  ggplot2::ylab("Commercial landings (mt)") +
  ggplot2::guides(
    colour = ggplot2::guide_legend(title = "area x fleet")
  ) +
  ggplot2::scale_colour_manual(
    values = unlist(t(get_fleet(value="Comm")[,c("col.n","col.s")]),use.names=FALSE)
  ) +
  ggplot2::theme_bw() + ggplot2::theme(legend.position = c(0.1, 0.8))
ggsave2(
  gg,
  filename = file.path("..", "figures", "catch-comm-CA-interpolate-ts.png"),
  caption = paste0(
    "Reconstructed commercial landings for the state of California by fleet ",
    "(fixed gear, FG, circles; trawl gear, TW, triangles) and ", 
    "area (northern region, blue; southern region, red). ",
    "Dashed line indicates data were interpolated across years."
  ),
  alttext = paste(
    ""
  )
)

gg <- catch_comm_CA_ralston %>%
  dplyr::group_by(Year, area, state) %>%
  dplyr::summarize(mt = sum(mt), .groups = "keep") %>%
  dplyr::ungroup() %>%
  dplyr::group_by(Year, state) %>%
  dplyr::mutate(prop = mt / sum(mt)) %>%
  dplyr::full_join(
    by = c("Year", "area"),
    x = .,
    y = catch_comm_CA_MM
  ) %>% data.frame %>%
  dplyr::filter(type == "ALL" & area == "North") %>%
ggplot2::ggplot(
  size = 2,
  ggplot2::aes(x = Year, y = prop)
) +
  ggplot2::geom_point(pch = 1) +
  ggplot2::geom_point(ggplot2::aes(y = prop_fa)) +
  ggplot2::ylim(c(0, 1.0)) +
  ggplot2::theme_bw() +
  ggplot2::ylab("Proportion of CA commercial catch in the northern area")
ggsave2(
  gg,
  filename = file.path("..", "figures", "catch-comm-CA-proportionnorth.png"),
  caption = paste0(
    "Yearly proportion of California (CA) commercial catch landed in the ",
    "the northern region (i.e., north of ",
    "forty degrees ten minutes latitude", ") ",
    "of all CA commercial catches from two data sources, ",
    "the Raltson et al. (2010) catch reconstruction effort (open circles) and ",
    "fish ticket data in CALCOM, the database used by the ",
    "California Cooperative Groundfish Survey to store and manage ",
    "commercial market sample data, (filled circles)."
  ),
  alttext = paste0(
    "Proportions largely agree between both data sets."
  ),
  label = "catch-comm-CA-proportionnorth"
)

gg <- ggplot2::ggplot(data = catch_comm_pacfin %>%
  dplyr::filter(Year < as.numeric(format(Sys.Date(), "%Y"))) %>%
  dplyr::group_by(Year, area, fleet) %>%
  dplyr::summarize(MT = sum(ROUND_WEIGHT_MTONS), .groups = "keep")) +
  ggplot2::geom_bar(stat="identity", lwd = 0.2,
    ggplot2::aes(
      x = Year,
      y = MT,
      fill = interaction(area, fleet, sep = " ")
    )
  ) +
  ggplot2::scale_fill_manual(
    values  = unlist(t(get_fleet(value="Comm")[,c("col.n","col.s")]),use.names=FALSE)
  ) +
  ggplot2::facet_grid(area ~ .) +
  ggplot2::xlab("Year") + ggplot2::ylab("Landings (mt)") +
  ggplot2::labs(fill = ggplot2::guide_legend(title = "area x fleet")) +
  ggplot2::theme_bw() + ggplot2::theme(
    strip.background = ggplot2::element_rect(fill = NA),
    text = ggplot2::element_text(size = 16)
    )
ggsave2(
  gg,
  filename = file.path("..", "figures", "catch-comm-ts-pacfin.png"),
  caption = paste0(
    "Yearly commercial landings (mt) from the PacFIN database since 1981 ",
    " by area (panel and color) and fleet (shading). ",
    "Trawl gear (TW) includes all trawl, nets, and dredging. ",
    "Fixed gear (FG) includes all other gear types."
  ),
  alttext = paste(
    ""
  )
)

gg <- catch_comm_pacfin %>%
  dplyr::filter(AGENCY_CODE == "C", Year != as.numeric(format(Sys.Date(), "%Y"))) %>%
  dplyr::mutate(gear = factor(fleet)) %>%
  dplyr::group_by(Year, fleet, area)  %>%
  dplyr::summarize(mt = sum(ROUND_WEIGHT_MTONS), .groups = "keep") %>%
  dplyr::group_by(Year, fleet) %>%
  dplyr::mutate(percentage = mt/sum(mt)) %>%
  ggplot2::ggplot(ggplot2::aes(
    x = Year,
    y = percentage,
    fill = interaction(area, fleet)
  )) +
  ggplot2::scale_fill_manual(
    values = unlist(t(get_fleet(value="Comm")[,c("col.n","col.s")]),use.names=FALSE)
  ) +
  ggplot2::facet_grid(fleet ~ .) +
  ggplot2::geom_bar(position = "stack", stat = "identity") +
  ggplot2::scale_y_continuous(labels = scales::percent) +
  ggplot2::theme_bw() + ggplot2::theme(
    strip.background = ggplot2::element_rect(fill = NA),
    # text = ggplot2::element_text(size = 10),
    legend.position = "none"
  ) +
  ggplot2::xlab("Year") + ggplot2::ylab("California commercial catch by area")
ggsave2(
  gg,
  filename = file.path("..", "figures", "catch-comm-CA-gearprop.png"),
  caption = paste0(
    "Percentage of California commercial catch by area within each fleet ",
    "(fixed gear, FG, upper panel; trawl, TW, lower panel) ",
    "since 1981. ",
    "Darker shades represent the northern area and lighter shades represent the southern area."
  ),
  alttext = paste(
    ""
  )
)

gg <- bio_rec_recfin %>%
  dplyr::group_by(Year, STATE_NAME, AGENCY_WEIGHT_UNITS, IS_RETAINED, sex) %>%
  dplyr::summarize(.groups = "keep",
    mean = mean(RECFIN_LENGTH_MM, na.rm = TRUE)
  ) %>%
  ggplot2::ggplot(ggplot2::aes(Year, mean, col = STATE_NAME, lty = STATE_NAME)) +
  ggplot2::geom_point(cex = 5, alpha = 0.5) +
  ggplot2::geom_smooth(method = "lm", se = FALSE) +
  ggplot2::scale_color_manual(
    values = get_fleet(value="rec")[,c("col.n")]
  ) +
  ggplot2::ylab("Mean length (mm)") +
  ggplot2::labs(col = "", lty = "") +
  ggplot2::theme_bw() +
  ggplot2::theme(
    legend.position = c(0.15, 0.90),
    strip.background = ggplot2::element_rect(fill = NA)
  ) +
  ggplot2::facet_grid(tolower(IS_RETAINED) ~ factor(sex, labels = c("female", "male", "unknown sex")))
ggsave2(
  gg,
  filename = file.path("..", "figures", "catch-rec-tsrecfinmeanlength.png"),
  caption = paste0(
    "Mean length (mm) of released (top panel) and retained (bottom panel) fish from the ",
    "recreational fishery for Washington (long dashes), ",
    "Oregon (short dashes), and California (solid line) ",
    " by sex (columns). ",
    "The average length of retained fish across all years was used to translate ",
    "numbers to weight for Washington."
  ),
  alttext = paste(
    ""
  )
)

gg <- ggplot2::ggplot(
  catch_rec %>%
    dplyr::group_by(Year, state) %>%
    dplyr::summarize(mt = sum(mt), .groups = "keep") %>%
    dplyr::ungroup(),
    ggplot2::aes(Year, mt)
  ) +
  ggplot2::geom_line(ggplot2::aes(col = state), lty = 2) +
  ggplot2::scale_color_manual(
    values = rev(get_fleet(value="rec")[,c("col.n")])
  ) +
  ggplot2::geom_line(
    lty = 1,
    col = get_fleet(value="CA")[,c("col.n")],
    data = last_ssdat[[2]][["catch"]] %>%
      dplyr::filter(fleet == grep("REC", last_ssdat[[2]][["fleetnames"]]), catch > 0),
    ggplot2::aes(year, catch)
  ) +
  ggplot2::geom_line(
    lty = 1,
    col = get_fleet(value="OR")[,c("col.n")],
    data = last_ssdat[[1]][["catch"]] %>%
      dplyr::filter(fleet == grep("OR", last_ssdat[[1]][["fleetnames"]]), catch > 0),
    ggplot2::aes(year, catch)
  ) +
  ggplot2::geom_line(
    lty = 1,
    col = get_fleet(value="WA")[,c("col.n")],
    data = last_ssdat[[1]][["catch"]] %>%
      dplyr::filter(fleet == grep("WA", last_ssdat[[1]][["fleetnames"]]), catch > 0),
    ggplot2::aes(year, catch)
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = c(0.1, 0.8)) +
  ggplot2::ylab("Recreational landings (mt)")
ggsave2(
  gg,
  filename = file.path("..", "figures", "catch-rec-ca-oldts.png"),
  caption = paste0(
    "Comparison of recreational landings from this assessment (dashed line) versus ",
    "the previous assessment (solid line) for each state (colors). ",
    "The previous assessment used numbers for WA recreational landings, whereas this assessment used weight."
  ),
  alttext = paste(
    ""
  )
)

gg <- ggplot2::ggplot(
  catch_rec_CA,
  ggplot2::aes(Year, mt, col = area)
) +
  ggplot2::geom_point(ggplot2::aes(pch = source)) +
  ggplot2::geom_line() +
  ggplot2::scale_colour_manual(
    values = unlist(use.names = FALSE, get_fleet("CA", col = c("col.n", "col.s")))
  ) +
  ggplot2::theme_bw() +
  ggplot2::ylab("Recreational California landings (mt)")
ggsave2(
  gg,
  filename = file.path("..", "figures", "catch-rec-CA-ts.png"),
  caption = paste0(
    "Time series of California recreational landings (mt) for the ",
    "northern (darker color) and southern (lighter color) areas. ",
    "The shape of the points indicates the information source ",
    "(California Recreational Fisheries Survey, CRFS; linear interpolation, interpolate; ",
    "Marine Recreational Fisheries Statistics Survey, MRFSS; ",
    "old Stock Synthesis, SS, model)."
  ),
  alttext = paste(
    ""
  )
)

gg <- ggplot2::ggplot(
  catch_rec,
  ggplot2::aes(Year, mt, col = interaction(area, state, sep = " "))
) +
  ggplot2::geom_line(size=1.2) +
  ggplot2::scale_color_manual(
    values = unlist(use.names = FALSE, get_fleet(value="rec")[c(3,3,2,1), c("col.n", "col.s")])[c(1,5,3,4)]
  ) +
  ggplot2::labs(col = "fleet") +
  ggplot2::ylab("Recreational landings (mt)") +
  ggplot2::theme_bw()
ggsave2(
  gg,
  filename = file.path("..", "figures", "catch-rec-ts.png"),
  caption = paste0(
    "Time series of recreational landings by state and area."
  ),
  alttext = paste(
    ""
  )
)

catch.pacfin %>%
  dplyr::select(matches("Year|^area$|Round.+MT")) %>%
  dplyr::mutate(Year = Year) %>%
  dplyr::group_by(Year, area) %>%
  dplyr::summarize(MT = sum(ROUND_WEIGHT_MTONS), .groups = "keep") %>%
  tidyr::pivot_wider(names_from = area, values_from = MT) %>%
  kableExtra::kbl(
    longtable = TRUE,
    caption = "Yearly commercial landings (mt) from the PacFIN database by area, North and South, since 1981."
  ) %>%
  kableExtra::kable_paper(full_width = TRUE) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "condensed", "responsive")) %>%
  kableExtra::scroll_box(width = "500px", height = "200px")

dplyr::full_join(
  by = c(Year = "year"),
  x = catch_rec %>%
    dplyr::group_by(Year) %>%
    dplyr::filter(state == "OR") %>%
    dplyr::summarize(mt=sum(mt)),
  y = last_ssdat[[1]][["catch"]] %>%
      dplyr::filter(fleet == grep("OR", last_ssdat[[1]][["fleetnames"]]), catch > 0)
  ) %>%
  dplyr::mutate(diff = mt - catch) %>%
  dplyr::select(Year, mt, catch, diff) %>%
  kable(.)
  data.frame()

gg <- ggplot2::ggplot(data = catch_comm_pacfin %>%
  dplyr::group_by(Year, LANDING_MONTH, area, fleet) %>%
  dplyr::summarize(.groups = "keep",
    MT = sum(ROUND_WEIGHT_MTONS),
    n = length(unique(VESSEL_ID))) %>%
  dplyr::mutate(fleet = ifelse(fleet == "FG", "Fixed gear", "Trawl")) %>%
  dplyr::ungroup() %>%
  dplyr::filter(n > 3),
  ggplot2::aes(y = MT, x = LANDING_MONTH, col = Year)) +
  ggplot2::xlab("") + ggplot2::ylab("Commercial landings (mt; log scale)") +
  ggplot2::geom_point(size=4, alpha = 0.3) +
  ggplot2::guides(alpha = "none") +
  ggplot2::theme_bw() + ggplot2::theme(
    legend.position = "bottom",
    legend.key.width = grid::unit(3, "cm"),
    strip.background = ggplot2::element_rect(fill = NA),
    text = ggplot2::element_text(size = 16)
    ) +
  ggplot2::labs(colour = ggplot2::guide_legend(title = "Year")) +
  ggplot2::coord_polar() +
  ggplot2::facet_grid(fleet~area) +
  ggplot2::scale_colour_viridis_c(direction = -1) +
  ggplot2::scale_x_continuous(
    limits = c(0, 12), expand = c(0, 0),
    breaks = seq(1, 12, by = 1),
    labels = c("Jan", "Feb", "Mar", "Apr", "May", "June", "July", "Aug",
      "Sep", "Oct", "Nov", "Dec")) +
  ggplot2::scale_y_continuous(trans = "log", breaks = c(0, 1, 10, 100, 1000))
ggsave2(
  gg,
  filename = file.path("..", "figures", "catch-comm-seasonality.png"),
  caption = paste0(
    "Exploratory figure, not ready for production. ",
    "Seasonality (month; x axis) of landings (mt) by area (column) and fleet (row). ",
    "The y axis is on the log scale, colors represent years, and ",
    "circles are transparent to facilitate visualization of months with many records."
  ),
  alttext = paste(
    ""
  )
)


gg <- ggplot2::ggplot(
  catch_comm,
  ggplot2::aes(Year, mt, col = interaction(area, fleet, sep = " "))
) +
  ggplot2::geom_line(size = 1.2) +
  ggplot2::scale_colour_manual(
    values = unlist(t(get_fleet(value="Comm")[,c("col.n","col.s")]),use.names=FALSE)
  ) +
  ggplot2::labs(ylab = "Landings (mt)", colour = ggplot2::guide_legend(title = "area x fleet")) +
  ggplot2::theme_bw()
ggsave2(
  gg,
  filename = file.path("..", "figures", "catch-comm-ts-full.png"),
  caption = paste0(
    "Commercial landings (mt) by year (x axis), area, and fleet, ",
    "with darker colors representing catch from the northern area and ",
    "lighter colors from the southern area."
  ),
  alttext = paste(
    "Recent landings are much smaller than pre-2000 landings for all commercial fleets."
  )
)

load(file.path("..", "data-raw", "catch.Rdata"))

Commercial landings

Commercial fleet structure

sa4ss::read_child("catch-comm-fleetstructure.Rmd")

cat(text_wcgop)

Reconstruction of commercial landings

```{asis, catch-comm-WA-text, opts.label = 'north', echo = TRUE}

Washington commercial reconstruction

The reconstruction of commercial landings for coastal waters off of Washington was provided by \gls{wdfw}. This reconstruction included landings starting in r min(catch_comm_WA[["Year"]]). Data from the reconstruction was used instead of data in \gls{pacfin} when there was overlap because \gls{wdfw} separates landings from each fish ticket by area. This is important for fish tickets that include landings from Alaskan, Canadian, Puget Sound, and oceanic waters. In \gls{pacfin}, it is more than likely that landings from different areas included on a single fish ticket would be assigned to just one area and partitioning out these landings to area would require accessing logbook information as well as fish ticket information, which is difficult and not part of the current \gls{pacfin} protocols. The reconstruction should largely match what is in \gls{pacfin} because \gls{pacfin} does not currently have landings for years prior to 1980 and Canadian waters have been closed to U.S. fishers targeting groundfish since 1978.

The reconstruction includes data from many sources, but consistently recorded data were largely available starting in 1943 from U.S. Fish Commission reports. Landings from prior to 1941 were converted to round fish weight using a conversion factor of 1.431 because it is assumed that reports were of filleted fish. Linear interpretation between adjacent years was used to fill in years without landings information. In the 2017 assessment, [@haltuch2019lingcod] missing years were filled forward rather than linearly interpolated (Figure \@ref(fig:catch-comm-state)). These differences were minor, and there were no major changes in the time series compared to what was used in the previous assessment.

```{asis, catch-comm-OR-reconstruction, opts.label = 'north', echo = TRUE}
##### Oregon commercial reconstruction

In Oregon, historical commercial landings from
`r min(catch_comm_OR[["Year"]])` to
`r max(catch_comm_OR[["Year"]])`
were provided by \gls{odfw} [@karnowski2014].
Historical landings began with exclusively longline landings,
which was the primary gear type
until the development of the trawl fishery in the 1940s.
Historical landings exhibited an increasing trend until peaking at
`r sprintf("%.0f mt in %.0f", catch_comm_OR_summary[["sum"]], catch_comm_OR_summary[["Year"]])`
and averaged `r sprintf("%.1f mt", catch_comm_OR_summary[["mean"]])` annually.

California commercial reconstruction

@sette1928 provided information from interviews and state records on fishing patterns from 1888 to 1926 for eight regions within U.S. waters. States along the Pacific Coast comprised one region, though state-specific landings were provided for Washington, Oregon, and California by species or species groups. For r utils_name("common"), the first positive record was from r catch_sette1928 %>% dplyr::filter(lbs > 0) %>% dplyr::pull(Year) %>% min and positive landings were documented for r length(catch_sette1928 %>% dplyr::filter(lbs > 0) %>% dplyr::pull(Year) %>% unique) - 1 years. We used linear interpolation to fill in years with missing data, ramping up from zero in r min(catch_sette1928[["Year"]]) to create a time series of r length(min(catch_sette1928[["Year"]]):max(catch_sette1928[["Year"]])) years (Figure \@ref(fig:catch-comm-CA-interpolate-ts)).

Catches by gear type were only available from 1926, and thus, the calculated proportion of of the landings caught by FG and TW in 1926 was applied to all years because it was assumed that ratios were similar across the time series. Proportionally, FG represented r catch_sette_prop_fa %>% dplyr::filter(fleet == "TW") %>% dplyr::pull(prop_fa) %>% sprintf("%.2f",.) %>% knitr::combine_words() of the total catch for the r catch_sette_prop_fa %>% dplyr::filter(fleet == "TW") %>% dplyr::pull(area) %>% tolower() %>% knitr::combine_words(.) areas, respectively. The Northern California district was assumed to represent north of 40 degrees ten minutes and all other districts combined were used to represent the southern area.

Landings from California fish market data, available within the ERDDAP database, were used to estimate the proportion of early landings that occurred in the northern area versus the southern area because fish market data were recorded by region on a yearly basis [@mason2004] within this data set. Whereas, @sette1928 only contained information on area for a single year. California fish market data represent a multi-organizational effort, but most landings are from fish ticket information collected by \gls{cdfw}. First, we calculated the yearly proportion of landings that occurred within the Eureka region north of Point Arena compared to all other regions r sprintf("(%.2f)", erddapmeanpropEureka) from port-specific landings from r min(catch_erddap[["Year"]]) to r sprintf("%s.", erddapmaxyrmeanbyarea) Second, the proportion of landings within Eureka region that occurred north versus south of Cape Mendocino r sprintf("(%.2f)", erddapmeanpropNEureka) was calculated from 100-200 block data [@Miller_2014] starting in 1925 to 1931. The product of the means of these two proportions was used to partition data from @sette1928 to area.

California fish market data, were available from the ERDDAP database over many years, but only those years that were missing between @sette1928 and the California Catch Reconstruction Project [@NOAA-TM-NMFS-SWFSC-461] were used. This resulted in keeping data from r paste0(paste(collapse = " to ", catch_erddap %>% dplyr::filter(Year < catch_comm_CA_ralston$Year %>% min) %>% dplyr::pull(Year) %>% range),".")

Information on region of landing was available and provided a means to assign the landings to the northern and southern areas. Though as previously mentioned, the Eureka region needed to be partitioned to area. We used the mean proportion of fish landed in the northern Eureka region r sprintf("(%.2f)", erddapmeanpropNEureka) from block data [@Miller_2014] to partition the sum of yearly landings within the Eureka region between areas.

@NOAA-TM-NMFS-SWFSC-461 represents the effort led by the \gls{swfsc} to reconstruct groundfish landings for the \gls{pfmc}, which are seen as the best available data for historical commercial landings from California ports. The data includes information on region of landing based on block assignments. Landings within region nine were assumed to be caught off of Mexico and were removed. Landings with a region code of two were partitioned to the northern and southern areas using the same method used above for California fish market landings. To check the validity of this assumption, we compared the proportion of landings assigned to the north versus the south to proportions calculated from confidential fish ticket data available in \gls{calcom} that have information on port of landing for available years between r paste(collapse = " - ", range(catch_comm_CA_MM[catch_comm_CA_MM[["type"]] == "ALL","Year"])) (pers. comm., M. Monk, \gls{swfsc}). The proportions showed similar trends, though the former were consistently higher than the latter for all years (Figure \@ref(fig:catch-comm-CA-proportionnorth)).

The @NOAA-TM-NMFS-SWFSC-461 data also had to be partitioned to fleet given it does not contain information about gear. Fish ticket information in \gls{calcom} was used to calculate the proportion of landings landed by each fleet for the northern and southern areas separately; thus, partitioning landings by year into four groups, northern TW, northern FG, southern TW, and southern FG. Proportions were only available for the following years: r knitr::combine_words(catch_comm_CA_MM[catch_comm_CA_MM[["type"]] == "ALL","Year"]). Years with no information were back filled using adjacent years.

Starting in r sprintf("%.0f,", min(catch_comm_CA_calcom[["Year"]])) commercial landings were available for California by port-group complex and gear group from \gls{calcom}. The following gear groups r knitr::combine_words(unique(catch_comm_CA_calcom %>% dplyr::filter(!GEAR_GRP %in% calcom_fleet_TW) %>% dplyr::pull(GEAR_GRP))) were combined to encompass FG and r knitr::combine_words(calcom_fleet_TW) gear groups were combined to encompass TW.

Unfortunately, the port-group complexes did not exactly align with the north-south split. But, it was assumed that the amount of landings within the Eureka port-group complex that occurred in the south was minor, and thus, all landings within the Crescent City and Eureka port-group complexes were assigned to the northern area and all other ports were assigned to the southern area.

``{asis, eval = params$model == "North", echo = TRUE} Landings caught off of the coast of Oregonr sprintf("(%.2f mt)", sum(catch_comm_CA_ORwaters[["Oregon"]]))and Washingtonr sprintf("(%.2f mt)", sum(catch_comm_CA_ORwaters[["Washington"]]))but landed in California during the periodr format_range(catch_comm_CA_ORwaters[["Year"]], parentheses = FALSE)` were added to the California reconstruction as was done in the previous assessment. These landings were assigned to the northern area TW fleet. In the future, the assignment of species and gear should be investigated more thoroughly for these landings.

<!-- Missing data -->
For combinations of year, area, and fleet that were missing
in the reconstruction of California commercial landings,
landings were interpolated based on a linear approximation
between adjacent years with data
(Figure \@ref(fig:catch-comm-CA-interpolate-ts)).
Thus, the reconstruction ramped up from zero starting in
`r catch_sette1928[["Year"]][1]`
to
`r sprintf("%.2f", catch_sette1928 %>% dplyr::filter(mt > 0) %>% dplyr::pull(mt) %>% dplyr::first())`
mt in
`r catch_sette1928 %>% dplyr::filter(mt > 0) %>% dplyr::pull(Year) %>% dplyr::first()`
and all subsequent missing years of data were filled in based on
linear interpolation between missing years
for a given area and fleet combination.

#### \glsentrylong{pacfin} (\glsentryshort{pacfin})

Commercial data were downloaded from the PacFIN database and provided
landings for Washington, Oregon, and California
(Figure \@ref(fig:catch-comm-ts-pacfin)).
These landings were treated as the best available information for
California since `r min(catch.pacfin[["Year"]])` and
for Washington and Oregon since the beginning of
`r sprintf("%s and %s,", 1 + max(catch_comm_WA[["Year"]]), 1 + max(catch_comm_OR[["Year"]]))`
respectively.

```{asis, catch-comm-OR-recent, opts.label = 'north', echo = TRUE}
Though `r utils_name("common")` are encountered in a wide variety of gear types in Oregon,
most landings are from bottom trawl gear types
`r sprintf("(%.1f percent from %.0f - %.0f).", catch_comm_pacfin_summary %>% dplyr::filter(state == "OR",fleet=="TW") %>% dplyr::pull(prop) %>% mean * 100, max(catch_comm_OR[["Year"]]) + 1, max(catch.pacfin[["Year"]]) - 1)`
Landings from all other gear types are minimal relative to TW and FG and relatively sporadic.
Commercial landings from Oregon waters between
`r sprintf("%.0f - %.0f", max(catch_comm_OR[["Year"]]), max(catch.pacfin[["Year"]]) - 1)`
peaked in
`r catch_comm_pacfin_summary %>% dplyr::filter(state == "OR") %>% dplyr::pull(maxYear) %>% dplyr::first()`
at
`r catch_comm_pacfin_summary %>% dplyr::filter(state == "OR") %>% dplyr::pull(max) %>% dplyr::first() %>% sprintf("%.2f", .)`
before declining and fluctuating between
`r catch_comm_pacfin_summary %>% dplyr::filter(state == "OR" & Year >= 2000) %>% dplyr::pull(yearsum) %>% range %>% sprintf("%.2f",.) %>% knitr::combine_words()`
mt since 2000.

Before splitting the commercial landings to area, all landings that were known to have been caught outside of the U.S. Exclusive Economic Zone r sprintf("(%.2f mt)", removenonEEZcatch) were removed. These were landings that occurred in an unknown \gls{inpfc} area noted as XX or \gls{pacfin} area 02 or 4A.

The split at \CapeM required finding a method for splitting data within the Eureka (ERA) port-group complex. Data with the port of landing of Shelter Cove, a port within ERA, were assigned to the southern model and data from all other ports within ERA, i.e., r paste0(knitr::combine_words(catch_comm_CA_ERA), ",") were assigned to the northern model. If landings were not assigned to a port-group complex, then the physical location of the port of landing was used to assign an area. Lastly, if both port-group complex and port of landing were unknown, then the landings were assigned to an area based on a vessel's historical behavior while landing r spp. This ad hoc method that was used to assign area of landing to a trivial amount of landings without spatial information (r sprintf("(%.2f mt)", sum(catch_comm_CA_noport[["sumMT"]]))) has almost zero influence on the resulting catch streams and was done largely as an exercise to document how to do this should the stock boundary change again and the numbers would need to be recalculated. For example, if more than half of a vessel's landings of r utils_name("common") were in ERA or CCA, then all of their landings without an assigned area were assigned to the northern area.

Recreational landings

Recreational fleet structure

Recreational data were first compiled to the state level as a single fleet with all gear types. Then, data from northern California was added to the northern model as its own fleet rather than being combined with Oregon recreational data because each fishery is subject to unique regulations that affect selectivity (Figure \@ref(fig:catch-rec-ts)). The remaining data from California was used to model recreational fisheries in the southern area as a single fleet.

```{asis, catch-rec-WA-text, opts.label = 'north', echo = TRUE}

Washington recreational landings

The time series of recreational catches (numbers of fish) were provided by \gls{wdfw} and included information on fishing within Washington's coastal areas sampled by the Ocean Sampling Program in Marine Areas 1-4. Catches that were landed within the Strait of Juan de Fuca, i.e., Bonilla-Tatoosh line to the mouth of the Sekiu River, and sampled by the Puget Sound Sampling Program were included in the landings because they were potentially caught in ocean waters even though they were landed in coastal waters. All other landings from the Puget Sound Sampling Program were excluded because they are from the Puget Sound rather than the Pacific Ocean.

First, we calculated the mean length r sprintf("(%0.2f mm)", bio_rec_recfin_meanlength %>% dplyr::filter(STATE_NAME == "WASHINGTON") %>% dplyr::pull(mean_length)) of fish landed within Washington recreational fishery (Figure \@ref(fig:catch-rec-tsrecfinmeanlength)) across all years r format_range(bio_rec_recfin %>% dplyr::filter(STATE_NAME == "WASHINGTON") %>% dplyr::pull(Year)) and sexes. Next, we used this mean length and the weight-length relationship as calculated from the most recent survey data to determine the mean weight. Finally, weight (mt) was determined from mean weight and numbers. Final landings include retained and 7\% of discarded fish to account for those that died. The previous assessment modeled Washington landings in terms of numbers of fish rather than weight (Figure \@ref(fig:catch-rec-ca-oldts)). A comparison among the two assessments showed similar patterns in the catch when compared in the same units by using the internally estimated values calculated within the final models which are based on the estimated values for growth, selectivity, and age composition in each year showed similar patterns (Figure \ref{fig:wa-rec-catch-comparison}).

```{asis, catch-rec-OR-text, opts.label = 'north', echo = TRUE}
#### Oregon recreational landings

The recreational fishery in Oregon likely began in the early 1950s or 1960s
but data on catches prior to
`r min(data_catch %>% dplyr::filter(fleet == "OR") %>% dplyr::pull("Year"))`
are not available.
ODFW provided the time series of catches that includes information on
shoreside activities, estuary boats, private ocean boats, and charter ocean boats.

<!-- Historical ocean boat landings (1974 – 2000) -->
Recently, \gls{odfw} undertook an effort to comprehensively reconstruct
all marine fish recreational ocean boat landings prior to 2001
(pers. comm., A. Whitman, \gls{odfw}).
Reconstructed catch estimates from \gls{orbs} improve upon estimates
from the federal \gls{mrfss},
which have known biases related to effort estimation and sampling
[@vanvoorhees2000]
that resulted in catch estimates considered implausible by \gls{odfw}.
However, the \gls{orbs} sample estimates are known to lack the
comprehensive spatial and temporal coverage of \gls{mrfss}.
Addressing this coverage issue is a major part of this reconstruction.
In general, the base data and methodology for these reconstructed
estimates are consistent with recent assessments for other nearshore species
[@dick2016china; @dick2018bluedeacon; @haltuch2019lingcod; @cope2019cabezon].

Prior to 2001, \gls{orbs} monitored marine species in
both multi-species categories, such as rockfish, flatfish, and
other miscellaneous fishes, and
as individual species, such as lingcod or halibut.
For this comprehensive reconstruction,
four species categories were selected to reconstruct,
including rockfish, lingcod, flatfish, and miscellaneous,
which constitute the bulk of the managed marine fish species.
Lingcod have been a single species category throughout this time period.

Category-level estimates were expanded to account for gaps in sampling coverage
in two separate pathways.
First, estimates from five major ports were expanded to include unsampled
winter months in years lacking complete coverage.
Expansions were based on available year-round sampling data and
excluded years where regulations may have impacted the temporal distribution of catch.
Second, all other minor port estimates were expanded to include
seasonal estimates in years lacking any sampling based on
the amount of minor port catch as compared to all major port estimates.
A subset of landings were sampled by \gls{orbs}
for species compositions within these categories.
Once category-level landings were comprehensive in space and time,
species compositions were applied for the three multi-species categories,
including rockfish, flatfish, and miscellaneous fish.
Borrowing rules for species compositions were specific to the category and
determined based on a series of regression tree analyses
that detailed the importance of each domain
(year, month, port, and fishing mode) to variability in compositions.

Ocean boat estimates from 1979 – 2000 in numbers of fish of
`r utils_name("common")` from the above described methods were
converted to biomass using biological samples from \gls{mrfss}
(pers. comm., A. Whitman, \gls{odfw}).
\Gls{mrfss} biological data are available from
1980 – 1989 and 1993 – 2000.
An annual average weight was applied to the total annual number of fish
to obtain an annual landings estimate.
Several years missing biological data (1979, 1990 – 1992)
were filled in using neighboring years or interpolation.
Landings from 1974 – 1978 were estimated using similar
methods to the above reconstruction,
but are not currently part of this reconstructed time period,
and were not updated from the 2017 assessment [@haltuch2019lingcod].
Updated landings include only those from 1979 – 2000.
These landings in biomass were provided by \gls{odfw}
and do not include an estimate of discards.
Landings during this time period fluctuate,
with a peak of 237.2 mt in 1993,
and fluctuate between approximately 50 and 200 mt
following that peak (1979 – 2000).

<!-- Modern ocean boat landings (2001 – 2020) -->
Recreational landings for ocean boat modes from 2001 – 2020
are available from \gls{recfin} (extracted on 03/19/2021).
Both retained and released estimates of mortality are included,
though retained mortality contributes the vast majority to total mortality.
Release mortality is estimated from angler-reported release rates and
the application of discard mortality rates from the PFMC.
From 2001 – 2020, landings averaged 139.4 mt,
ranging from 60.1 mt in 2001 to 229.1 mt in 2015.
In 2020, Oregon ocean boat landings were 172.5 mt.

<!-- Shore and estuary landings (1980 – 2020) -->
\Gls{odfw} provided reconstructed estimates of shore and estuary
landings for `r spp` from 1980 – 2020,
using methodology similar to recent assessments
[@berger2015kelpgreenling; @dick2018bluedeacon; @cope2019cabezon].
Data sources include \gls{mrfss} and \gls{sebs}.
Numbers of fish were provided by \gls{mrfss}
from 1980 – 1989 and 1993 – June 2003, and
by \gls{sebs} from July 2003 – June 2005.
An annual mode-specific average weight was applied to numbers of
`r utils_name("common")` from 1980 – 1989 and 1993 – 2005.
Separate weights were calculated for shore and estuary boat modes and
excluded extreme outliers and imputed values.
This reconstruction also applied two scaling factors to
remove bias towards freshwater sampling and
underestimation of estuary boats [@dick2018bluedeacon].
To estimate `r utils_name("common")` landings from July – December 2005,
an expansion was developed using the three year average of the ratio between
the first six months of the year and the total annual landings from
\gls{mrfss} and \gls{sebs} landings from 2002 - 2004.
Separate expansions were developed for shore mode and estuary boat modes.

The \gls{odfw} does not currently sample shore and estuary boat fishing trips,
so a 10 year average landing (1996 – 2005; 6.7 mt/year)
was used to estimate shore and estuary boat landings during 2006 – 2020.
Shore and estuary boat landings combined fluctuate but gradually increase
until peaking in 1993 at 16.7 mt.
Shore and estuary boat landings average 6.5 mt annually from 1980 – 2003.
Shore and estuary landings were combined with the ocean boat landings
for the total Oregon recreational landings (1974 – 2020).

California recreational landings

California recreational r utils_name("common") catches since r min(catch_rec_1980[["Year"]]) - 1 are available within the \gls{mrfss} database (Figure \@ref(fig:catch-rec-CA-ts)) and stored in the \gls{recfin} database. Data were provided in terms of weight, and thus, converting from numbers to weight using mean length ((Figure \@ref(fig:catch-rec-tsrecfinmeanlength))) was not needed. The first year of data are typically not used because of the lack of standardization within the sampling protocols which led to vastly different estimates of catches compared to later years. Thus, r min(catch_rec_1980[["Year"]]) is used as the first year of \gls{mrfss} data. Data were provided by John Field for years prior to r min(catch_rec_1980[["Year"]]) and these data have been unchanged since the 2009 assessment of r utils_name("common") [@hamel2009lingcod].

For this assessment, we had to split the historical data provided by John Field and \gls{mrfss} to area. This was accomplished using data from Albin et al. (1993) that includes county-specific estimates of landings. Area-specific landings were informative about the proportion of landings in Del Norte and Humboldt county relative to the rest of the California coast. A catch-weighted mean proportion for the years r knitr::combine_words(unique(data_albin$Year)) was used to split coast-wide recreational landings to area.

Between 1990 and 1992, the \gls{mrfss} sampling program ceased because of budgetary limitations. Additionally, there was a transition year between \gls{mrfss} and \gls{crfs} without data. Thus, linear interpolation was used to provide proxy estimates for r knitr::combine_words(catch_rec_CA %>% dplyr::filter(source == "interpolate") %>% dplyr::pull(Year) %>% unique) (Figure \@ref(fig:catch-rec-CA-ts)).

Sampling under \gls{crfs} started in January of r catch_rec_2000$Year %>% min and data are still currently being collected. Information includes data on port group that was used to partition landings to area. Redwood was assigned to the northern area for all years since 2005, even though Redwood in 2005 through 2007 also contained landings from Shelter Cove. This time series also includes landings from Mexico and Canada that were excluded from this analysis.

Comparison of current catch stream to previous model

The current time series of catches were aggregated using the previous model structure, i.e., the northern area included Washington and Oregon and the southern area included California, for comparison purposes (Figure \@ref(fig:catch-comm-state)). This comparison revealed differences in the current catch stream for the commercial FG and TW fleets compared to what was previously used for the southern model. It was determined that the 2017 catch stream for the California TW fleet included catches from both TW and FG fleets. This has since been corrected, and the current catch stream represents what is known to be the best available data.

{asis, opts.label = 'north', echo = TRUE} Differences in the reconstructions for commercial catches are the result of linearly interpolating historical data rather than assuming it was equal to the previous year (Figure \@ref(fig:catch-comm-state)). Differences in the early reconstruction for California commercial catches are the result of using landings from @sette1928 instead of ramping up catches from zero (Figure \@ref(fig:catch-comm-state)).

There were small differences in the recreational catches from California and Oregon but the main differences were in Washington recreational catches. As mentioned previously, Washington recreational catches were converted from numbers to weight prior to fitting them as data, whereas the previous assessment fit them as numbers (Figure \ref{fig:catch-rec-ca-oldts}). Differences in Oregon recreational catches were the result of \gls{odfw}'s efforts to create comprehensive sport fishing reconstructions. These reconstructions are noted to be more robust than information from \gls{mrfss} and explain the differences in the trajectories prior to 2000.

iantaylor-NOAA/Lingcod_2021 documentation built on Oct. 30, 2024, 6:42 p.m.