From a133974f54550e48d9740cfc28adb6b484a84a37 Mon Sep 17 00:00:00 2001 From: okenk Date: Thu, 12 Dec 2024 16:00:06 -0800 Subject: [PATCH] start observer index, minor data updates --- Rscripts/logbook_analys.R | 83 +++++++++++++++++++++++++++++++++++++++ docs/data_summary_doc.qmd | 30 ++++---------- docs/index.html | 15 +++---- 3 files changed, 96 insertions(+), 32 deletions(-) create mode 100644 Rscripts/logbook_analys.R diff --git a/Rscripts/logbook_analys.R b/Rscripts/logbook_analys.R new file mode 100644 index 0000000..242794a --- /dev/null +++ b/Rscripts/logbook_analys.R @@ -0,0 +1,83 @@ +library(dplyr) + +observer <- read.csv('data/confidential/wcgop_logbook/ward_oken_observer_data_2024-09-10.csv') |> + as_tibble() |> + filter(AVG_LAT > 40+1/6, + sector == 'Midwater Rockfish', + DATATYPE == 'Analysis Data') + +observer_yt <- observer|> + group_by(HAUL_ID) |> + summarise(has_yt = 'Yellowtail Rockfish' %in% species, + across(AVG_DEPTH:YEAR, first)) |> + filter(!has_yt) |> + mutate(species = 'Yellowtail Rockfish', MT = 0) |> + select(-has_yt) |> + bind_rows( + filter(observer, species == 'Yellowtail Rockfish') + ) |> + mutate(DRVID = as.integer(DRVID)) + +em <- read.csv('data/confidential/wcgop_logbook/ward_oken_em_data_2024-12-03.csv') |> + as_tibble() |> + filter(AVG_LAT > 40+1/6, sector == 'Midwater Rockfish EM') + +em_yt <- em |> + group_by(HAUL_ID) |> + summarise(has_yt = 'Yellowtail Rockfish' %in% species, + across(AVG_DEPTH:YEAR, first)) |> + filter(!has_yt) |> + mutate(species = 'Yellowtail Rockfish', MT = 0) |> + select(-has_yt) |> + bind_rows( + filter(em, species == 'Yellowtail Rockfish') + ) + +all_hauls_yt <- bind_rows(list(obs = observer_yt, em = em_yt), .id = 'program') |> + filter(HAUL_DURATION > 0, HAUL_DURATION < 10) |> + mutate(depth_std = (AVG_DEPTH - mean(AVG_DEPTH))/sd(AVG_DEPTH), + fyear = factor(YEAR), log_duration = log(HAUL_DURATION), + cpue = MT/HAUL_DURATION, month = lubridate::month(as.Date(SET_DATE)), + fDRVID = factor(DRVID), fport = factor(R_PORT)) # note D_PORT is NA in some cases, R_PORT is always defined (fish ticket) + +yt_mesh <- sdmTMB::make_mesh(all_hauls_yt, xy_cols = c('AVG_LONG', 'AVG_LAT'), n_knots = 50) +plot(yt_mesh) +# this does not look like a st model is a good idea. very clear "hotspots" of fishing effort. +# but we need it for sdmTMB to work (not actually used) + +mod_simple <- sdmTMB::sdmTMB(formula = cpue ~ fyear + s(depth_std) + s(month, bs = 'cc') + program + (1|fDRVID) + (1|fport) - 1, family = sdmTMB::delta_lognormal(), + data = all_hauls_yt, spatial = 'off', spatiotemporal = 'off', time = 'fyear', mesh = yt_mesh) + +summary(mod_simple) + + +resids <- residuals(mod_simple, model = 2) # the presence/absence model looks consistently good +qqnorm(resids) +qqline(resids) + +sdmTMB::visreg_delta(mod_simple, xvar = "depth_std", model = 1, gg = TRUE) +sdmTMB::visreg_delta(mod_simple, xvar = "depth_std", model = 2, gg = TRUE) +sdmTMB::visreg_delta(mod_simple, xvar = "month", model = 1, gg = TRUE) +sdmTMB::visreg_delta(mod_simple, xvar = "month", model = 2, gg = TRUE) + +preds <- predict(mod_simple, + newdata = tibble(fyear = factor(2012:2023), depth_std = 0, month = 6, program = 'obs', fDRVID = '546053', fport = 'NEWPORT'), + return_tmb_object = TRUE, re_form = NA, re_form_iid = NA) + +ind <- sdmTMB::get_index(preds, bias_correct = TRUE) + +ind |> + ggplot(aes(x = as.numeric(fyear), y = est, ymin = lwr, ymax = upr)) + + geom_line() + + geom_errorbar() + + +all_hauls_yt |> + ggplot() + + geom_point(aes(x = HAUL_DURATION, y = MT)) + geom_density(aes(x = HAUL_DURATION, col = YEAR, group = YEAR)) + +### general observation: as I have added terms to the model, the signal in the index went from well-defined to noisy and meaningless. +### the port (8 levels) and vessel (46 levels) effects really did it in. should think more carefully about what terms to include. +### note 3962 hauls, 2762 positive hauls over 12 years. ~330 hauls/yr, 230 positive hauls/yr. effort varies a lot by year. +### (increases 2012-2018) \ No newline at end of file diff --git a/docs/data_summary_doc.qmd b/docs/data_summary_doc.qmd index 4db0fbd..f97fa4a 100644 --- a/docs/data_summary_doc.qmd +++ b/docs/data_summary_doc.qmd @@ -32,7 +32,7 @@ yt_n_catch_pacfin <- catch.pacfin |> filter( AGENCY_CODE == 'W' | AGENCY_CODE == 'O' | COUNTY_NAME == 'DEL NORTE' | COUNTY_NAME == 'HUMBOLDT', # OR, WA, N. CA LANDING_YEAR < 2024, - !(CATCH_AREA_CODE %in% c('23A', '23C', '26A')) # exclude 48 Puget Sound fish tickets + !(CATCH_AREA_CODE %in% c('23A', '23C', '26A')) # exclude Puget Sound fish tickets ) # OR landings time series @@ -45,7 +45,7 @@ or_comm_catch <- read.csv(here("Data/Confidential/commercial/Oregon Commercial l # WA URCK/UPOP reconstruction, use prior to 2000 wa_urck <- read.csv(here('Data/Raw_not_confidential/WA_URCK_UPOP_recon_1981-2016.csv')) |> filter(Species_AfterComps %in% c('YTRK', 'YTR1'), - Year < 2000) |> + Year < 2001) |> group_by(Year) |> summarise(catch_mt = sum(RoundPounds_AfterComps) * 0.453592 / 1000) |> mutate(AGENCY_CODE = 'W') |> @@ -54,7 +54,7 @@ wa_urck <- read.csv(here('Data/Raw_not_confidential/WA_URCK_UPOP_recon_1981-2016 # Put it all together comm_draft_landings <- yt_n_catch_pacfin |> filter(AGENCY_CODE != 'O', - AGENCY_CODE != 'W' | PACFIN_YEAR >= 2000 + AGENCY_CODE != 'W' | PACFIN_YEAR >= 2001 ) |> group_by(PACFIN_YEAR, AGENCY_CODE)|> summarise(catch_mt = sum(ROUND_WEIGHT_MTONS)) |> @@ -115,25 +115,7 @@ comm_draft_landings |> geom_line(aes(x = Year, y = catch_mt, col = AGENCY_CODE, linetype = model), linewidth = 1) + viridis::scale_color_viridis(discrete = TRUE) + labs(title = 'Comparing commercial landings in 2017 model vs current', x = 'YEAR') - -# compare 2017 landings to current landings, WA only -comm_draft_landings |> - rename(Year = YEAR) |> - mutate(model = 'PacFIN') |> - filter(AGENCY_CODE == 'W') |> - bind_rows( - filter(comm2017, AGENCY_CODE == 'W') - )|> - bind_rows( - mutate(wa_urck, model = 'URCK/UPOP') - ) |> - ggplot() + - geom_line(aes(x = Year, y = catch_mt, col = model), linewidth = 1) + - viridis::scale_color_viridis(discrete = TRUE) + - labs(title = 'Possible sources for WA landings (vert. line = 2000)', x = 'YEAR') + - geom_vline(xintercept = 2000) ``` -URCK/UPOP = sum of RoundPounds_AfterComps for Species_AfterComps = YTR1 or YTRK, converted to mt. This is used *instead of* PacFIN catches prior to 2000. Table of commercial landings by state, 1981 onward. @@ -147,9 +129,11 @@ comm_draft_landings |> knitr::kable(align = 'l', digits = 1) ``` -WDFW has alerted the STAT these do not include all tribal catches in recent years. +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 + +Oregon landings were provided directly to the STAT from ODFW and come from PacFIN and a UPOP/URCK reconstruction. -Note this is for catch landed into Del Norte and Humboldt counties only. +California landings are for catch landed into Del Norte and Humboldt counties only. ### At-Sea landings diff --git a/docs/index.html b/docs/index.html index 25c41bc..635ac01 100644 --- a/docs/index.html +++ b/docs/index.html @@ -6,7 +6,7 @@ - + Data summary