get_prep_fishing.qmd

---
title: "Downloading and preparing fishing for ISIMIP 3a Baltic mizer"
author: "Max Lindmark, following Phoebe.Woodworth-Jefcoats@noaa.gov"
date: today
date-format: iso
toc: true
format: 
  html:
    page-layout: full
    embed-resources: true
knitr: 
  opts_chunk:
    fig.align: center
    out-width: 100%
editor: source
---

This document describes the workflow for preparing the ISIMIP 3a fishing

## Data Access
This script is used to prepare the fishing forcing using the data provided from the FishMIP coordinators (as opposed to the
data provided on the ISIMIP server).

## Data downloading and wrangling

```{r}
#| message: false
library(ggsidekick)
library(tidyverse); theme_set(theme_sleek())
library(tidylog)
library(patchwork)

home <- here::here()
```

```{r}
#| message: false
all_effort <- read_csv("effort_histsoc_1841_2010_regional_models.csv") %>% 
  filter(region == "Baltic.Sea")

ggplot(all_effort, aes(Year, NomActive, fill = Gear)) +
  geom_bar(stat = "identity") +
  theme(legend.position = "bottom") +
  coord_cartesian(expand = 0)
```

```{r}
# Test if catches also have spikes
# read_csv(paste0(home, "/baltic/FishMIP/FishingForcing/DataFiles/calibration_catch_histsoc_1850_2004_regional_models.csv")) %>%
#   filter(Sector == "industrial",
#          region == "Baltic.Sea",
#          FGroup == "demersal30-90cm") %>% 
#   summarise(rep_tot = sum(Reported), .by = "Year") %>% 
#   ggplot(aes(Year, rep_tot)) + 
#   geom_line()
```

This plot shows the effort by gear in the Baltic Sea across all Fgroups. One thing to notice is the reconstructed effort has a very steep increase, essentially going from zero around 1940's. Early exploration also shows that effort trends very much follow each other across these gear groups. Also, the composition species groups are very similar in bottom trawl and all other gears. For instance, even in the bottom trawl (cod panel above), pelagic is the majority group in most years. Therefore, I filter by both F group and the gears that are realistically catching these species. To know that, I check Swedish Loggbook data from fisherman and read the ICES working group WGBFAS reports.

![](FishMIP_gear_sp.png)

Filter the effort and Fgroup data to match the mizer model (main fisheries for cod, sprat and herring)

```{r}
#| message: false
#| fig-height: 9
# sort(unique(all_effort$Sector))
# sort(unique(all_effort$SAUP))
# sort(unique(all_effort$Gear))
# sort(unique(all_effort$FGroup))

names(all_effort)
unique(all_effort$Gear)

all_effort_sum <- all_effort %>% 
  filter(Sector == "Industrial" &
           Gear %in% c("Trawl_Bottom", "Trawl_Midwater_or_Unsp") &
           FGroup %in% c("demersal30-90cm", "pelagic<30cm")
         ) %>%
  mutate(sp = ifelse(FGroup %in% c("demersal30-90cm") & Gear %in% c("Trawl_Bottom"),
                     "Cod", NA),
         sp = ifelse(FGroup %in% c("pelagic<30cm") & Gear %in% c("Trawl_Midwater_or_Unsp"),
                     "Sprat/Herring", sp)) %>% 
  group_by(sp, Year) %>% 
  summarise(sum_NomActive = sum(NomActive)) %>% 
  ungroup() %>% 
  drop_na()

ggplot(all_effort_sum, aes(Year, sum_NomActive)) + 
  #geom_bar(stat = "identity") +
  geom_line() +
  coord_cartesian(xlim = c(1945, 2014), expand = 0) + 
  facet_wrap(~sp, scales = "free", ncol = 1) +
  scale_fill_brewer(palette = "Set3") +
  theme(legend.position = "bottom")
```

The above plot shows the total effort for cod and sprat + herring. I've assumed that if the gear is bottom trawl, gillnets or lines, and the Fgroup is demersal < 90 cm,  it's cod effort. If gear is midwater trawl and Fgroup is pelagic <30, it's pelagic effort 

Following the Hawaii model, I now sum *NomActive* to get two time series between 1961--2010. Then I calculate the average effort between 1992--2004, because the original model calibrated to average conditions between 1992--2002, but here we are asked to use up to year 2004 for calibrating. I next find the maximum fishing effort in the calibration 

```{r}
#| message: false
all_effort_agg <- all_effort_sum %>% 
  filter(Year >= 1961 & Year <= 2010) %>% 
  group_by(sp, Year) %>% 
  summarise(agg_NomActive = sum(sum_NomActive)) %>%
  ungroup()

# Now I need to split sprat and herring into two groups
all_effort_agg_cod <- all_effort_agg %>% filter(sp == "Cod")

all_effort_agg <- all_effort_agg_cod %>% 
  bind_rows(all_effort_agg %>% filter(sp == "Sprat/Herring") %>% mutate(sp = "Sprat")) %>% 
  bind_rows(all_effort_agg %>% filter(sp == "Sprat/Herring") %>% mutate(sp = "Herring"))

cal_tp_avg <- all_effort_agg %>% 
  filter(Year >= 1992 & Year <= 2004) %>% 
  group_by(sp) %>% 
  summarise(max_NomActive = max(agg_NomActive)) %>% 
  ungroup()

# Plot time series of effort by species group
all_effort_agg <- all_effort_agg %>% 
  left_join(cal_tp_avg, by = "sp") %>% 
  mutate(scaled_NomActive_max = agg_NomActive/max_NomActive)

# Plot the max-scaled efforts
ggplot(all_effort_agg, aes(Year, scaled_NomActive_max, color = sp, linetype = sp)) + 
  geom_rect(xmin = 1992, xmax = 2004, ymin = -Inf, ymax = Inf, color = NA,
            fill = "grey90", alpha = 0.05,
            inherit.aes = FALSE) + 
  geom_line(linewidth = 0.9) + 
  facet_wrap(~sp, scales = "free", ncol = 2) +
  scale_linetype_manual(values = c(3, 1, 2)) +
  scale_color_brewer(palette = "Dark2", name = "") +
  guides(linetype = "none") +
  theme(legend.position = "bottom") +
  labs(caption = "FishMIP summed efforts before scaling, calibration time period in rectangle, horisontal lines are averages")
```

Read in the fishing efforts I used in the calibration and compare with the FishMIP effort:

```{r}
#| fig-height: 7
#| message: false

# Load effort
projectEffort <- read_csv("projectEffort.csv")[, 2:4]

# Add in Year
projectTemp_old <- read_csv("projectTemp.csv")[, 2:3]

projectEffort <- projectEffort %>% 
  mutate(year = projectTemp_old$Year) %>% 
  filter(year >= 1961 & year <= 2010) %>% 
  pivot_longer(cols = c("Cod", "Herring", "Sprat"), names_to = "sp", values_to = "effort") %>% 
  mutate(source = "assessment")

# Add in FishMip effort
all_effort_agg_assess_max <- all_effort_agg %>%
  rename(effort = scaled_NomActive_max,
         year = Year) %>% 
  dplyr::select(sp, year, effort) %>% 
  mutate(source = "FishMIP_assessment F max")

projectEffort <- bind_rows(projectEffort, all_effort_agg_assess_max)

# Now scale the maximum to match the magnitude of the assessment F
projectEffort %>% 
  ggplot(aes(year, effort, color = source)) + 
  geom_line() + 
  facet_wrap(~sp, ncol = 2, scales = "free") +
  theme_sleek(base_size = 9) + 
  scale_color_brewer(palette = "Set1", name = "") +
  labs(x = "Year", y = "Effort") + 
  NULL

# What's the best scaling factor? Fit a regression
d <- projectEffort %>% 
  dplyr::select(year, sp, effort, source) %>% 
  pivot_wider(names_from = "source", values_from = "effort") %>% 
  janitor::clean_names()

fit <- lm(assessment - fish_mip_assessment_f_max ~ 0 + sp, data = d) %>% 
  broom::tidy()

fit

# Re-scale
projectEffort %>% 
  filter(source == "FishMIP_assessment F max") %>% 
  arrange(effort)

projectEffort <- projectEffort %>% 
  mutate(effort = ifelse(source == "FishMIP_assessment F max" & sp == "Cod",
                         fit$estimate[fit$term == "spCod"] + effort,
                         effort),
         effort = ifelse(source == "FishMIP_assessment F max" & sp == "Sprat",
                         fit$estimate[fit$term == "spSprat"] + effort,
                         effort),
         effort = ifelse(source == "FishMIP_assessment F max" & sp == "Herring",
                         fit$estimate[fit$term == "spHerring"] + effort,
                         effort)
         ) %>% 
  # Because we estimate the correction factor without forcing it to be 0 we can end up with negative efforts in some years for herring. Fix this!
  mutate(effort = ifelse(effort < 0, 0, effort))

# This is for trimming quantiles
# ggplot(projectEffort, aes(effort)) + 
#   geom_histogram() + 
#   facet_wrap(~source, ncol = 1)
# 
# upr <- projectEffort %>% 
#   filter(source == "FishMIP_assessment F" & sp == "Cod") %>% 
#   summarise(upr = quantile(effort, probs = 0.9), .by = sp)
# 
# projectEffort <- projectEffort %>% 
#   mutate(effort = effort, 
#          # effort = ifelse(sp == "Cod" &
#          #                   source == "FishMIP_assessment F" &
#          #                   effort > as.numeric(upr$upr),
#          #                 as.numeric(upr$upr), effort)
#          )

# Check correlation
projectEffort <- projectEffort %>% 
  filter(!source == "FishMIP_assessment F") %>% 
  mutate(source2 = ifelse(source == "assessment", "Regional assessment F", "Global"))

#devtools::install_github("rensa/stickylabeller")
library(stickylabeller)

cors <- projectEffort %>% 
  dplyr::select(-source) %>% 
  pivot_wider(names_from = source2, values_from = effort) %>% 
  summarise(cor = as.numeric(cor.test(`Regional assessment F`, `Global`)$estimate),
            cor_p = as.numeric(cor.test(`Regional assessment F`, `Global`)$p.value),
            .by = "sp") %>% 
  mutate(cor = round(cor, digits = 3))

cors 

projectEffort <- projectEffort %>% left_join(cors, by = "sp")

ggplot(projectEffort, aes(year, effort, color = source2)) + 
  geom_line() + 
  facet_wrap(
    ncol = 3,
    ~ sp + cor,
    labeller = label_glue(
      '{sp}\nCor. = {cor}')) +
  scale_color_brewer(palette = "Set1", name = "Effort source") +
  labs(x = "Year", y = "Scaled effort") + 
  theme_sleek(base_size = 10) +
  theme(aspect.ratio = 1,
        legend.position.inside = c(0.8, 0.84),
        axis.text.x = element_text(size = 7)) +
  guides(color = guide_legend(position = "inside", title.hjust = 0)) + 
  NULL

ggsave("scaled_effort.png", width = 17, height = 7, unit = "cm")
```

In summary, it seems that the FishMIP effort is quite different from the species-specific F derived from assessments for cod, but not that much for sprat and herring. After filtering relevant Fgroups and fishing gears, I get good match in trends and absolute values for the pelagics, but less so for cod. However, I might be able to use this since the magnitude is probably within reasonable bounds (slightly high, though the high values are spiky and perhaps it isn't for long enough to crash the stock)