formatting catches for ss3, starting pacfin comps

Rscripts/pacfin_ashop_comps.R

@@ -0,0 +1,37 @@
+# PacFIN comps ------------------------------------------------------------
+
+load(here('Data/Confidential/Commercial/pacfin_12-11-2024/PacFIN.YTRK.bds.11.Dec.2024.RData'))
+
+comm_catch <- readRDS('Data/processed/ss3_landings_2023.rds') |>
+  filter(fleet == 1, year >= 1981) |>
+  select(year, catch)
+
+bds_clean <- PacFIN.Utilities::cleanPacFIN(bds.pacfin, keep_age_method = c('B', 'S')) |>
+  filter(age_method != 'L', # n = 1
+         age_method != 'T', # n = 2,
+         age_method != '', # n = 33
+         age_method != 'S' | state != 'WA' # n = 99, keep OR surface reads for now
+  )
+
+w_l_pars <- read.csv('Data/processed/W_L_pars.csv')
+
+# length expansions
+pacfin_exp1 <- PacFIN.Utilities::getExpansion_1(Pdata = bds_clean, 
+                                                fa = w_l_pars$A[w_l_pars$Sex == 'F'],
+                                                fb = w_l_pars$B[w_l_pars$Sex == 'F'], 
+                                                ma = w_l_pars$A[w_l_pars$Sex == 'M'],
+                                                mb = w_l_pars$B[w_l_pars$Sex == 'M'],
+                                                ua = w_l_pars$A[w_l_pars$Sex == 'B'],
+                                                ub = w_l_pars$B[w_l_pars$Sex == 'B'])
+
+pacfin_exp2 <- pacfin_exp1 |>
+  mutate(fleet = 'catch') |>
+  PacFIN.Utilities::getExpansion_2(Catch = comm_catch, 
+                                   Units = 'MT', 
+                                   stratification.cols = 'fleet') |>
+  mutate(Final_Sample_Size = PacFIN.Utilities::capValues(Expansion_Factor_1_L * Expansion_Factor_2, 
+                                                         maxVal = 0.80)) # where does the 0.8 come from?
+
+comps <- PacFIN.Utilities::getComps(pacfin_exp2, 
+                                    strat = 'fleet', 
+                                    Comps = 'LEN')

docs/data_summary_doc.qmd

@@ -111,17 +111,6 @@ comm_draft_landings <- yt_n_catch_pacfin |>
   group_by(YEAR, AGENCY_CODE) |>
   summarise(catch_mt = sum(catch_mt))
 
-# 2017 assessment landings
-comm2017 <- read.csv(here("Data/Raw_not_confidential/commercial_catch_2017.csv")) |> 
-  as_tibble() |>
-  tidyr::replace_na(list(Mirick.reconstructed.POP1.URCK = 0)) |>
-  mutate(O = OR.Recent..MT. + Mirick.reconstructed.POP1.URCK + OR.Hist..MT.,
-         W = as.numeric(stringr::str_remove(WA, ','))) |> 
-  select(Year, W, C = CA.North..MT., O) |>
-  filter(Year >= 1981) |>
-  tidyr::pivot_longer(-Year, names_to = 'AGENCY_CODE', values_to = 'catch_mt') |>
-  mutate(model = '2017')
-
 # landings by state
 comm_draft_landings |>
   ggplot() +
@@ -150,6 +139,18 @@ yt_n_catch_pacfin |>
   labs(fill = "PacFIN gear group\n(MDT = midwater trawl)", x = 'YEAR') #+
   # viridis::scale_fill_viridis(discrete = TRUE) # colors were hard to see in viridis
 
+# 2017 assessment landings
+comm2017 <- read.csv(here("Data/Raw_not_confidential/commercial_catch_2017.csv")) |> 
+  as_tibble() |>
+  tidyr::replace_na(list(Mirick.reconstructed.POP1.URCK = 0)) |>
+  mutate(O = OR.Recent..MT. + Mirick.reconstructed.POP1.URCK + OR.Hist..MT.,
+         W = as.numeric(stringr::str_remove(WA, ',')),
+         C = CA.North..MT. + CA.North.Hist..MT.) |> 
+  select(Year, W, C, O) |>
+  # filter(Year >= 1981) |>
+  tidyr::pivot_longer(-Year, names_to = 'AGENCY_CODE', values_to = 'catch_mt') |>
+  mutate(model = '2017') 
+
 # compare 2017 landings to current landings, by state
 comm_draft_landings |>
   rename(Year = YEAR) |>
@@ -177,7 +178,13 @@ WDFW has alerted the STAT these do not include all tribal catches in recent year
 
 Oregon landings were provided directly to the STAT from ODFW and come from PacFIN and a UPOP/URCK reconstruction. 
 
-California landings are for catch landed into Del Norte and Humboldt counties only.
+California landings are for catch landed into Del Norte and Humboldt counties only. Data sources are:
+
+- Ralston et al. (2010) + supplementary catch from OR waters: 1916 - 1968
+- CALCOM: 1969 - 1980
+- PacFIN: 1981 - present
+
+In the Ralston et al. (2010) reconstruction, the northernmost region (region 2) includes Crescent City, Eureka AND Fort Bragg port area groups. Catch in the Crescent City and Eureka port area groups (i.e., Del Norte and Humboldt counties) is estimated based on the fraction of catch in Crescent City and Eureka divided by catch in all of region 2 during the first five years of CALCOM. The vast majority of the catch is in the two northern port areas (`r round(calcom_ratio$pct_n_40_10 * 100, 0)`\%). All supplementary catches in the Ralston reconstruction caught off the coast of Oregon, but landed into California are included. 
 
 ### At-Sea landings
 
@@ -188,7 +195,7 @@ ashop_catch <- readxl::read_excel(
   mutate(catch_mt = EXPANDED_SumOfEXTRAPOLATED_2SECTOR_WEIGHT_KG / 1000)
 
 # compare to 2017
-dat2017 <- r4ss::SS_readdat(here('Model_Runs/base_2017_3.30.23/YTRK.NORTH.data.ss'))
+dat2017 <- r4ss::SS_readdat(here('Model_Runs/1.02_base_2017_3.30.23/YTRK.NORTH.data.ss'))
 
 ashop2017 <- dat2017$catch |>
   filter(fleet == 2) |> 
@@ -362,7 +369,7 @@ ralston_rec <- readxl::read_excel(
   rename(RECFIN_YEAR = Year)
 
 # Put it all together
-bind_rows(
+rec_landings <- bind_rows(
   or_rec_catch, 
   wa_modern, 
   tidyr::pivot_longer(wa_historical, cols = -RECFIN_YEAR, names_to = 'State', values_to = 'Dead_Catch'),
@@ -371,7 +378,26 @@ bind_rows(
   ca_mrfss_catch, 
   ca_mrfss_interpolate,
   ralston_rec
-  ) |>
+  ) 
+
+# rec_landings |>
+#   filter(State != 'WA (num)') |>
+#   mutate(fleet = ifelse(State == 'WA (mt)', 'WA', 'OR/CA')) |>
+#   group_by(RECFIN_YEAR, fleet) |>
+#   summarise(Dead_Catch = sum(Dead_Catch)) |>
+#   mutate(model = 'current') |>
+#   bind_rows({
+#     filter(dat2017$catch, fleet %in% 3:4) |>
+#       mutate(fleet = ifelse(fleet == 3, 'OR/CA', 'WA'), 
+#              model = '2017') |>
+#       select(RECFIN_YEAR = year, Dead_Catch = catch, fleet, model)
+#   }) |>
+#   ggplot() +
+#   geom_line(aes(x = RECFIN_YEAR, y = Dead_Catch, col = fleet, linetype = model),
+#             linewidth = 1)
+# This is bizarre, it is as if WA was actually entered in mt, even though it is specified in numbers.
+
+rec_landings |>
   tidyr::pivot_wider(names_from = 'State', values_from = 'Dead_Catch', values_fill = 0) |>
   rename(YEAR = 'RECFIN_YEAR') |>
   arrange(YEAR) |>
@@ -413,6 +439,40 @@ For CA rec catch estimates for the assessment:
 
 In addition, Albin ratio (1) is used to calculate what fraction of N. CA catches from the Ralston recreational reconstruction are in the assessment area.
 
+```{r ss3-landings}
+rec_landings_ss3 <- rec_landings |>
+  filter(State != 'WA (num)') |>
+  group_by(RECFIN_YEAR) |>
+  summarise(catch = sum(Dead_Catch)) |>
+  rename(YEAR = RECFIN_YEAR) 
+
+comm_landings_ss3 <- comm2017 |>
+  filter(AGENCY_CODE == 'W', Year < 1981) |>
+  rename(YEAR = Year) |>
+  bind_rows(comm_draft_landings) |> 
+  select(-model) |>
+  group_by(YEAR) |>
+  summarise(catch = sum(catch_mt))
+
+ss3_landings <- bind_rows(list(
+  comm_landings_ss3,
+  select(ashop_catch, YEAR, catch = catch_mt),
+  rec_landings_ss3
+  ), .id = 'fleet') |>
+  mutate(seas = 1, catch_se = 0.05) |>
+  rename(year = YEAR) |>
+  select(names(dat2017$catch)) |>
+  as.data.frame()
+
+ss3_landings |>
+  saveRDS(file = here('Data/Processed/ss3_landings_2023.rds'))
+
+ss3_landings |>
+  filter(year < 2017) |>
+  saveRDS(file = here('Data/Processed/ss3_landings_2016.rds'))
+```
+
+
 ## Comp Data
 
 ### Commercial
@@ -421,7 +481,10 @@ Initial length sample sizes after running `PacFIN.Utilities::cleanPacFIN()`:
 
 ```{r comm-bio, cache=TRUE}
 load(here('Data/Confidential/Commercial/pacfin_12-11-2024/PacFIN.YTRK.bds.11.Dec.2024.RData'))
-bds_clean <- PacFIN.Utilities::cleanPacFIN(bds.pacfin, verbose = FALSE) 
+bds_clean <- PacFIN.Utilities::cleanPacFIN(bds.pacfin)#, verbose = FALSE) 
+
+# to figure out, why is clean pacfin filtering out OR ages before 1981.
+
 bds_clean |>
   filter(state == 'WA' | state == 'OR' | PACFIN_GROUP_PORT_CODE == 'CCA' |
            PACFIN_GROUP_PORT_CODE == 'ERA', geargroup != 'XXX') |> 
@@ -452,10 +515,6 @@ bds_clean |>
   knitr::kable(align = 'l', digits = 1)
 ```
 
-Note that California lengths are not yet available for 2023.
-
-While running `PacFIN.Utilities::cleanPacFIN()`, it noted that 20 Washington ages in 1991 were double read, the second reader was not recorded, the first and second ages differed, and no final age was determined.
-
 #### Comparing midwater vs bottom trawl lengths and ages
 
 The length and age distribution associated with the midwater trawl (MDT) code in PacFIN is is very similar to bottom trawl catch (everything else in the TWL gear group). Note that the sample sizes for ages in 2022-2024 are lower than other years, especially midwater in 2022 which only has 14 ages.

docs/index.html

@@ -6,7 +6,7 @@
 
 <meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
 
-<meta name="dcterms.date" content="2024-12-19">
+<meta name="dcterms.date" content="2024-12-26">
 
 <title>Data summary</title>
 <style>
@@ -77,7 +77,7 @@ <h1 class="title">Data summary</h1>
     <div>
     <div class="quarto-title-meta-heading">Published</div>
     <div class="quarto-title-meta-contents">
-      <p class="date">December 19, 2024</p>
+      <p class="date">December 26, 2024</p>
     </div>
   </div>
 
@@ -93,7 +93,7 @@ <h2 class="anchored" data-anchor-id="landings">Landings</h2>
 <section id="commercial-landings" class="level3">
 <h3 class="anchored" data-anchor-id="commercial-landings">Commercial landings</h3>
 <p>Some summary figures:</p>
-<div class="cell" data-hash="data_summary_doc_cache/html/pacfin-catch_5385a239b4438f689e1ca1c596ff6665">
+<div class="cell" data-hash="data_summary_doc_cache/html/pacfin-catch_9f5ec1d8276585739cdbad259ee8c05e">
 <div class="cell-output-display">
 <p><img src="data_summary_doc_files/figure-html/pacfin-catch-1.png" class="img-fluid" width="672"></p>
 </div>
@@ -918,7 +918,13 @@ <h3 class="anchored" data-anchor-id="commercial-landings">Commercial landings</h
 </div>
 <p>WDFW has alerted the STAT these do not include all tribal catches in recent years. 1981-2000 are from the UPOP/URCK reconstruction, summing catches for YTRK and YTR1. 2001 onward comes from PacFIN</p>
 <p>Oregon landings were provided directly to the STAT from ODFW and come from PacFIN and a UPOP/URCK reconstruction.</p>
-<p>California landings are for catch landed into Del Norte and Humboldt counties only.</p>
+<p>California landings are for catch landed into Del Norte and Humboldt counties only. Data sources are:</p>
+<ul>
+<li>Ralston et al.&nbsp;(2010) + supplementary catch from OR waters: 1916 - 1968</li>
+<li>CALCOM: 1969 - 1980</li>
+<li>PacFIN: 1981 - present</li>
+</ul>
+<p>In the Ralston et al.&nbsp;(2010) reconstruction, the northernmost region (region 2) includes Crescent City, Eureka AND Fort Bragg port area groups. Catch in the Crescent City and Eureka port area groups (i.e., Del Norte and Humboldt counties) is estimated based on the fraction of catch in Crescent City and Eureka divided by catch in all of region 2 during the first five years of CALCOM. The vast majority of the catch is in the two northern port areas (96%). All supplementary catches in the Ralston reconstruction caught off the coast of Oregon, but landed into California are included.</p>
 </section>
 <section id="at-sea-landings" class="level3">
 <h3 class="anchored" data-anchor-id="at-sea-landings">At-Sea landings</h3>
@@ -1827,14 +1833,14 @@ <h3 class="anchored" data-anchor-id="recreational-landings">Recreational landing
 <h4 class="anchored" data-anchor-id="method-for-converting-historical-wa-catches-from-numbers-to-weight">Method for converting historical WA catches from numbers to weight</h4>
 <p>We are planning to use WA (mt) for all years, but show WA (num) to compare with raw data.</p>
 <p>From WA sport biodata, calculate average length of samples collected prior to 1990. This included 48 samples from 1979, 1981, and 1982. The next available samples are in 1995, so we did not use length data from samples in years shortly after the historical reconstruction applies.</p>
-<p>Calculate average weight using the unsexed version of the weight-length relationship calculated from WCGBTS data. Could try triennial instead?</p>
+<p>Calculate average weight using the unsexed version of the weight-length relationship calculated from WCGBTS data. (Could use triennial instead?)</p>
 <p>Would WDFW prefer we</p>
 <ol type="1">
 <li>linearly interpolate between 1967 and 1975,</li>
 <li>linearly interpolate between 1967 and average of 1975-1977 (similar to the approach for missing CA MRFSS years), or</li>
 <li>leave those years as having zero catch?</li>
 </ol>
-<p>Current table of catches uses approach (2).</p>
+<p>Current table of catches uses approach (2). If WDFW has no strong opinion, this is what we plan to do.</p>
 </section>
 <section id="method-for-ca-mrfss-catches" class="level4">
 <h4 class="anchored" data-anchor-id="method-for-ca-mrfss-catches">Method for CA MRFSS catches:</h4>
@@ -1852,6 +1858,9 @@ <h4 class="anchored" data-anchor-id="method-for-ca-mrfss-catches">Method for CA
 <li>1990-1992: Interpolate between average of 1987-1989 and 1993-1995</li>
 </ul>
 <p>In addition, Albin ratio (1) is used to calculate what fraction of N. CA catches from the Ralston recreational reconstruction are in the assessment area.</p>
+<div class="cell">
+
+</div>
 </section>
 </section>
 </section>
@@ -1860,7 +1869,17 @@ <h2 class="anchored" data-anchor-id="comp-data">Comp Data</h2>
 <section id="commercial" class="level3">
 <h3 class="anchored" data-anchor-id="commercial">Commercial</h3>
 <p>Initial length sample sizes after running <code>PacFIN.Utilities::cleanPacFIN()</code>:</p>
-<div class="cell" data-hash="data_summary_doc_cache/html/comm-bio_2668abac0d099bc72a3c1078e7275134">
+<div class="cell" data-hash="data_summary_doc_cache/html/comm-bio_eebe1cf2f1e9f284a668a3671d0a636d">
+<div class="cell-output cell-output-stdout">
+<pre><code>
+   HKL    MSC    NET    POT    TLS    TWL    TWS    XXX 
+ 19247      7   2325      3    805 217645   8951     45 </code></pre>
+</div>
+<div class="cell-output cell-output-stdout">
+<pre><code>
+    CA     OR     WA 
+ 28026  94946 126206 </code></pre>
+</div>
 <div class="cell-output cell-output-stdout">
 <pre><code>  AGE_METHOD1 AGE_METHOD2 AGE_METHOD3 Age method for best age      n
 1           B           B           B                       B   2346
@@ -1873,6 +1892,13 @@ <h3 class="anchored" data-anchor-id="commercial">Commercial</h3>
 8           T        &lt;NA&gt;        &lt;NA&gt;                       T      2
 9        &lt;NA&gt;        &lt;NA&gt;        &lt;NA&gt;                    &lt;NA&gt;  82197</code></pre>
 </div>
+<div class="cell-output cell-output-stdout">
+<pre><code>  strings            region    n
+1      3D       CAN (VNCVR)  441
+2      4A Sound and Straits 1015
+3      5A               CAN 7618
+4      5B               CAN 6194</code></pre>
+</div>
 <div class="cell-output-display">
 <table class="table table-sm table-striped">
 <thead>
@@ -2691,8 +2717,6 @@ <h3 class="anchored" data-anchor-id="commercial">Commercial</h3>
 </table>
 </div>
 </div>
-<p>Note that California lengths are not yet available for 2023.</p>
-<p>While running <code>PacFIN.Utilities::cleanPacFIN()</code>, it noted that 20 Washington ages in 1991 were double read, the second reader was not recorded, the first and second ages differed, and no final age was determined.</p>
 <section id="comparing-midwater-vs-bottom-trawl-lengths-and-ages" class="level4">
 <h4 class="anchored" data-anchor-id="comparing-midwater-vs-bottom-trawl-lengths-and-ages">Comparing midwater vs bottom trawl lengths and ages</h4>
 <p>The length and age distribution associated with the midwater trawl (MDT) code in PacFIN is is very similar to bottom trawl catch (everything else in the TWL gear group). Note that the sample sizes for ages in 2022-2024 are lower than other years, especially midwater in 2022 which only has 14 ages.</p>