StockStatus2021.Rmd

---
title: "2021 SOE stock status indicator"
author: "Sarah Gaichas"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    code_fold: hide
    toc: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

remotes::install_github("https://github.com/NOAA-EDAB/stocksmart")

library(here)
library(tidyverse)
library(DT)
library(stocksmart)
library(ggrepel)
library(ecodata)
library(ggiraph)

```


## Pull from `stocksmart` for 2022 SOE reports

Andy renamed the assessmentdata package [stocksmart](https://noaa-edab.github.io/stocksmart/) based on [Stock SMART](https://www.st.nmfs.noaa.gov/stocksmart?app=homepage). 

Two data frames are in the package, `stockAssessmentData` and `stockAssessmentSummary`.

In `stockAssessmentData` we have time series. Columns are `r names(stockAssessmentData)` and the reported metrics are `r unique(stockAssessmentData$Metric)`. 


```{r}

datatable(head(stockAssessmentData), rownames = FALSE)

```

In `stockAssessmentSummary` we have assessment metadata. Columns are `r (names(stockAssessmentSummary))`.

```{r}

datatable(head(stockAssessmentSummary), rownames = FALSE, options = list(scrollX = TRUE))

```

Build ecodata input spreadsheet from `stockAssessmentSummary` and use the ecodata code to make the dataset for plotting:

```{r make-2021assess}

assess2021 <- stockAssessmentSummary %>%
  filter(`Science Center` == "NEFSC") %>%
  select(c(`Stock Name`, Jurisdiction, FMP, `Science Center`, 
           `Stock Area`, `Assessment Year`, `Last Data Year`,
           `F Year`, `Estimated F`, Flimit, Fmsy, `F/Flimit`, 
           `F/Fmsy`, Ftarget, `F/Ftarget`, `B Year`, `Estimated B`,
           `B Unit`, Blimit, Bmsy, `B/Blimit`, `B/Bmsy`)) %>%
  arrange(Jurisdiction, `Stock Name`, FMP, `Assessment Year`) %>%
  rename(Entity.Name = `Stock Name`,
         Assessment.Year = `Assessment Year`,
         F.Fmsy = `F/Fmsy`,
         B.Bmsy = `B/Bmsy`)

write.csv(assess2021, here("assess.csv"))

# from get_stocks.R, ecodata 2020
  #assess <- read.csv(file.path(data.dir, "2019assess.csv"))
  assess <- assess2021
  #decode <- read.csv(file.path(data.dir, "2019decoder.csv"))
  decode <- read.csv(here("2020decoder.csv"))
  
write.csv(decode, here("decoder.csv"))

stock_status_stockSMART <-
    assess %>%
    dplyr::group_by(Entity.Name) %>%
    dplyr::filter(Assessment.Year == max(Assessment.Year)) %>%
  #Find last year assessment occurred for each stock
    dplyr::ungroup() %>%
    dplyr::left_join(.,decode, by = "Entity.Name") %>% #Join in list of managed species
    dplyr::select(Entity.Name, Assessment.Year, F.Fmsy, B.Bmsy, Council, Code) %>%
  #select column variables to keep
    dplyr::mutate(id = 1:length(Entity.Name)) %>%
    tidyr::gather(.,Var, Value,-id,-Entity.Name,-Assessment.Year,-Council,-Code) %>%
  #wide to long
    dplyr::select(-id) %>%
    dplyr::rename(`Last assessment` = Assessment.Year,
                  Stock = Entity.Name) %>% #rename variables for clarity
    dplyr::mutate(Units = "unitless") #%>%
    #dplyr::mutate(Value = replace(Value, which(Code == "N Windowpane" & Var == "F.Fmsy"), NA))



```

Then test to see if we see the updates relative to 2020 ecodata. I'm leaving out all the plot annotations for unknown status. 

### Comparisons {.tabset}

#### StockSMART 2021 source, Mid-Atlantic

```{r, MAFMCplot}

stock_status <- 
  stock_status_stockSMART %>%
  mutate(Code = recode(Code, "Dogfish" = "Sp. Dogfish" )) %>% 
  spread(.,Var,Value) %>% 
  filter(Council %in% c("MAFMC","Both")) %>% 
  group_by(Stock) %>% 
  mutate(score = case_when(
    (B.Bmsy <0.5) ~"a",
    (F.Fmsy >1) ~ "a", 
    (F.Fmsy < 1 & B.Bmsy > 0.5 & B.Bmsy < 1) ~ "b",
    (F.Fmsy < 1 & B.Bmsy > 1) ~ "c"))
#Plot constants
y.max <- 2.1 #2.0 mackerel cut off F/Fmsy is 2.08
x.max <- 2.6
#A dataframe that defines custom legend for stocks with unknown status
# unknown <- data.frame(text = c("Unknown Status", "Longfin Squid",
#                               "Shortfin Squid", "N. Goosefish", "S. Goosefish"),
#                     x = rep(0.9*x.max,5), y = seq(0.93*y.max,1.5,-.1))

# Custom Color
custom_color<- c("#56B4E9", "#009E73", "#0072B2")
#Plotting code
ggplot(data = stock_status) +
  geom_vline(xintercept = 1, linetype = "dotted")+
  geom_vline(xintercept = 0.5, linetype = "dashed")+
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_point(aes(x = B.Bmsy,
                 y = F.Fmsy,
                 shape = Council,
                 color = score)) +
  geom_text_repel(aes(x = B.Bmsy, #geom_text_repel auto-jitters text around points
                      y = F.Fmsy,
                      label = Code, 
                      color = score), 
                  show.legend = FALSE, nudge_y = -0.01, nudge_x = 0.05) +
  scale_color_brewer(palette = "Dark2",
                     breaks = stock_status$score) +
  ylim(0,y.max) +
  xlim(0,x.max) +
  # geom_text(data = unknown, aes(x = x, y = y, label = text), #Custom legend for unknown stock status
  #           size = c(4.75,rep(4,4))) +
  # annotate("rect", xmin = 0.8*x.max,
  #          xmax = x.max,
  #          ymin = 0.65*y.max,
  #          ymax = 0.90*y.max,
  #          alpha = 0.1) +
  xlab(expression(~B/B[msy])) +
  ylab(expression(~F/F[msy])) +
  guides(color = FALSE) +
  theme_ts()
```

#### SOE 2020 ecodata source, Mid-Atlantic

```{r, MAFMC2020plot}
stock_status <- 
  ecodata::stock_status %>%
  mutate(Code = recode(Code, "Dogfish" = "Sp. Dogfish" )) %>% 
  spread(.,Var,Value) %>% 
  filter(Council %in% c("MAFMC","Both")) %>% 
  group_by(Stock) %>% 
  mutate(score = case_when(
    (B.Bmsy <0.5) ~"a",
    (F.Fmsy >1) ~ "a", 
    (F.Fmsy < 1 & B.Bmsy > 0.5 & B.Bmsy < 1) ~ "b",
    (F.Fmsy < 1 & B.Bmsy > 1) ~ "c"))
#Plot constants
y.max <- 2.0 #1.75 mackerel cut off F/Fmsy is 1.8
x.max <- 2.6
#A dataframe that defines custom legend for stocks with unknown status
# unknown <- data.frame(text = c("Unknown Status", "Longfin Squid",
#                               "Shortfin Squid", "N. Goosefish", "S. Goosefish"),
#                     x = rep(0.9*x.max,5), y = seq(0.93*y.max,1.4,-.1))

# Custom Color
custom_color<- c("#56B4E9", "#009E73", "#0072B2")
#Plotting code
ggplot(data = stock_status) +
  geom_vline(xintercept = 1, linetype = "dotted")+
  geom_vline(xintercept = 0.5, linetype = "dashed")+
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_point(aes(x = B.Bmsy,
                 y = F.Fmsy,
                 shape = Council,
                 color = score)) +
  geom_text_repel(aes(x = B.Bmsy, #geom_text_repel auto-jitters text around points
                      y = F.Fmsy,
                      label = Code, 
                      color = score), 
                  show.legend = FALSE, nudge_y = -0.01, nudge_x = 0.05) +
  scale_color_brewer(palette = "Dark2",
                     breaks = stock_status$score) +
  ylim(0,y.max) +
  xlim(0,x.max) +
  # geom_text(data = unknown, aes(x = x, y = y, label = text), #Custom legend for unknown stock status
  #           size = c(4.75,rep(4,4))) +
  # annotate("rect", xmin = 0.8*x.max,
  #          xmax = x.max,
  #          ymin = 0.65*y.max,
  #          ymax = 0.90*y.max,
  #          alpha = 0.1) +
  xlab(expression(~B/B[msy])) +
  ylab(expression(~F/F[msy])) +
  guides(color = FALSE) +
  theme_ts()
```


#### StockSMART 2021 source, New England

```{r, NEFMCtest}

stock_status <- stock_status_stockSMART %>%
  mutate(Code = recode(Code, "Dogfish" = "Sp. Dogfish" )) %>% 
  spread(.,Var,Value) %>% 
  filter(Council %in% c("NEFMC","Both")) %>% 
  group_by(Stock) %>% 
  mutate(score = case_when(
    (B.Bmsy <0.5) ~"a",
    (F.Fmsy >1) ~ "a", 
    (F.Fmsy < 1 & B.Bmsy > 0.5 & B.Bmsy < 1) ~ "b",
    (F.Fmsy < 1 & B.Bmsy > 1) ~ "c"))
#Plot constants
y.max <- 1.5
x.max <- 10
all_missing <- stock_status %>%
  filter(is.na(B.Bmsy),is.na(F.Fmsy)) %>% 
  dplyr::select(Code, Council)
b_missing <- stock_status %>%
  filter(is.na(B.Bmsy), !is.na(F.Fmsy)) %>% 
  dplyr::select(Code, Council)
f_missing <- stock_status %>%
  filter(is.na(F.Fmsy), !is.na(B.Bmsy)) %>% 
  dplyr::select(Code, Council)
#A dataframe that defines custom legend for stocks with unknown status
# all.df <- data.frame(text = all_missing$Code,
#                     x = rep(x.max*0.9,length(all_missing$Code)),
#                     #y = seq(1.45,1.05, length.out = 7))
#                     y = seq(1.45,1.05, length.out = length(all_missing$Code)))
# b.df <- data.frame(text = b_missing$Code,
#                     x = rep(x.max*0.7,length(b_missing$Code)),
#                     y = c(1.45,2.15, length.out = length(b_missing$Code)))
# f.df <- data.frame(text = f_missing$Code,
#                     x = rep(x.max*0.5,length(f_missing$Code)),
#                     y = seq(1.45,1.0, length.out = length(f_missing$Code)))

# Custom Color
custom_color<- c("#56B4E9", "#009E73", "#0072B2")

#Plotting code
ggplot(data = stock_status) +
  geom_vline(xintercept = 1, linetype = "dotted", color = "grey60")+
  geom_vline(xintercept = 0.5, linetype = "dashed", color = "grey60")+
  geom_hline(yintercept = 1, linetype = "dashed", color = "grey60") +
  geom_point(aes(x = B.Bmsy,
                 y = F.Fmsy,
                 color = stock_status$score)) +
  geom_text_repel(aes(x = B.Bmsy, #geom_text_repel auto-jitters text around points
                      y = F.Fmsy,
                      label = Code,
                      color = stock_status$score), show.legend = FALSE,nudge_y = -0.01, nudge_x = 0.05) +
  ylim(0,y.max) +
  xlim(0,x.max*1.1) +
  # geom_text(data = all.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  # geom_text(data = b.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  # geom_text(data = f.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  # scale_color_brewer(palette = "Dark2", #Change legend labels for clarity
  #                    breaks = stock_status$score) +
  # annotate("rect", xmin = 0.924*x.max,
  #          xmax = 1.08*x.max,
  #          ymin = 0.645*y.max,
  #          ymax = 0.98*y.max,
  #          alpha = 0.01) +
  # annotate("text", x = 9, y = 1.5, label = "F and B missing", fontface =2, size = 3)+
  # annotate("rect",  
  #            xmin = 0.70*x.max,
  #            xmax = 0.85*x.max,
  #            ymin = 0.90*y.max,
  #            ymax = 1.8,
  #            alpha = 0.01) +
  # annotate("text", x = 7, y = 1.5, label = "B missing", fontface =2, size = 3)+
  # annotate("rect", xmin = 0.509*x.max,
  #          xmax = 0.681*x.max,
  #          ymin = 0.65*y.max,
  #          ymax = 0.98*y.max,
  #          alpha = 0.01) +
  # annotate("text", x = 5, y = 1.5, label = "F missing", fontface =2, size = 3)+
  xlab(expression(~B/B[msy])) +
  ylab(expression(~F/F[msy])) +
  guides(color = FALSE) +
  theme_ts()

```


#### SOE 2020 ecodata source, New England

```{r, NEFMC2020plot}

stock_status <- ecodata::stock_status %>%
  mutate(Code = recode(Code, "Dogfish" = "Sp. Dogfish" )) %>% 
  spread(.,Var,Value) %>% 
  filter(Council %in% c("NEFMC","Both")) %>% 
  group_by(Stock) %>% 
  mutate(score = case_when(
    (B.Bmsy <0.5) ~"a",
    (F.Fmsy >1) ~ "a", 
    (F.Fmsy < 1 & B.Bmsy > 0.5 & B.Bmsy < 1) ~ "b",
    (F.Fmsy < 1 & B.Bmsy > 1) ~ "c"))
#Plot constants
y.max <- 1.5
x.max <- 10
all_missing <- stock_status %>%
  filter(is.na(B.Bmsy),is.na(F.Fmsy)) %>% 
  dplyr::select(Code, Council)
b_missing <- stock_status %>%
  filter(is.na(B.Bmsy), !is.na(F.Fmsy)) %>% 
  dplyr::select(Code, Council)
f_missing <- stock_status %>%
  filter(is.na(F.Fmsy), !is.na(B.Bmsy)) %>% 
  dplyr::select(Code, Council)
#A dataframe that defines custom legend for stocks with unknown status
# all.df <- data.frame(text = all_missing$Code,
#                     x = rep(x.max*0.9,length(all_missing$Code)),
#                     y = seq(1.45,1.05, length.out = length(all_missing$Code)))
# b.df <- data.frame(text = b_missing$Code,
#                     x = rep(x.max*0.7,length(b_missing$Code)),
#                     y = c(1.45,2.15, length.out = length(b_missing$Code)))
# f.df <- data.frame(text = f_missing$Code,
#                     x = rep(x.max*0.5,length(f_missing$Code)),
#                     y = seq(1.45,1.0, length.out = length(f_missing$Code)))

# Custom Color
custom_color<- c("#56B4E9", "#009E73", "#0072B2")

#Plotting code
ggplot(data = stock_status) +
  geom_vline(xintercept = 1, linetype = "dotted", color = "grey60")+
  geom_vline(xintercept = 0.5, linetype = "dashed", color = "grey60")+
  geom_hline(yintercept = 1, linetype = "dashed", color = "grey60") +
  geom_point(aes(x = B.Bmsy,
                 y = F.Fmsy,
                 color = stock_status$score)) +
  geom_text_repel(aes(x = B.Bmsy, #geom_text_repel auto-jitters text around points
                      y = F.Fmsy,
                      label = Code,
                      color = stock_status$score), show.legend = FALSE,nudge_y = -0.01, nudge_x = 0.05) +
  ylim(0,y.max) +
  xlim(0,x.max*1.1) +
  # geom_text(data = all.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  # geom_text(data = b.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  # geom_text(data = f.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  # scale_color_brewer(palette = "Dark2", #Change legend labels for clarity
  #                    breaks = stock_status$score) +
  # annotate("rect", xmin = 0.924*x.max,
  #          xmax = 1.08*x.max,
  #          ymin = 0.645*y.max,
  #          ymax = 0.98*y.max,
  #          alpha = 0.01) +
  # annotate("text", x = 9, y = 1.5, label = "F and B missing", fontface =2, size = 3)+
  # annotate("rect",  
  #            xmin = 0.70*x.max,
  #            xmax = 0.85*x.max,
  #            ymin = 0.90*y.max,
  #            ymax = 1.8,
  #            alpha = 0.01) +
  # annotate("text", x = 7, y = 1.5, label = "B missing", fontface =2, size = 3)+
  # annotate("rect", xmin = 0.509*x.max,
  #          xmax = 0.681*x.max,
  #          ymin = 0.65*y.max,
  #          ymax = 0.98*y.max,
  #          alpha = 0.01) +
  # annotate("text", x = 5, y = 1.5, label = "F missing", fontface =2, size = 3)+
  xlab(expression(~B/B[msy])) +
  ylab(expression(~F/F[msy])) +
  guides(color = FALSE) +
  theme_ts()


```


### {-}

### Issues

MAFMC looks reasonable based on 2021 updates to assessments.

NEFMC may a problem; silver hake stocks look better in 2021 than 2020. No update on the stock assessment report site but maybe hake assessments don't appear there?

Both hakes were assessed in 2020. The 2020 S Silver Hake assessment shows them as not overfished. So the updated stocksmart info looks correct. 

There was an error in the entity name field in the preliminary stocksmart dataset that caused S Silver Hake 2020 assessment results to not merge properly... we were showing 2017 in the previous version of ecodata. The new stocksmart data is correct for silver hake. 

N Windowpane appears because I didn't erase he Fmsy info in stocksmart as the previous ecodata get_stocks() function does. This may be appropriate; it did not appear in the 2020 plot but I don't remember why. Need to recheck.

#### Decisions

Leave all current stocksmart data in the 2021assess.csv and submit it. Ask for review by PDB prior to publishing.

#### 2021 local stock assessment source; ensure that stockSMART is up to date

Used the form here: https://apps-nefsc.fisheries.noaa.gov/saw/sasi/sasi_report_options.php

Files should be up-to-date as of the day of my query, November 16, 2021, according to the website.

Checked the boxes: Year--2021

Species with 2021 updates: Bluefish, Atlantic cod, Atlantic mackerel, Black sea bass, Golden tilefish, Scup, Summer flounder. 

The only ones that have stock status are Mid-Atlantic species and they look like reasonable changes. 

## More fun with stocksmart

### What is total (federally managed) catch from the ecosystem?

I think we should be able to summarize catches from `stocksmart::stockAssessmentData` by RegionalEcosystem:

`r unique(stocksmart::stockAssessmentData$RegionalEcosystem)`

But the issue is converting to common units before summing.

```{r NEUScatch}

NEcatch <- stocksmart::stockAssessmentData %>%
  filter(RegionalEcosystem == "Northeast Shelf",
         Metric == "Catch") %>%
  group_by(StockName) %>%
  filter(AssessmentYear == max(AssessmentYear)) %>%
  ungroup() %>%
  mutate(Catchmt = case_when(Units == "Thousand Metric Tons" ~ Value*1000,
                             TRUE ~ Value))
  
p <- ggplot(NEcatch, aes(x=Year, y=Catchmt, fill = StockName)) +
  #geom_bar(stat="identity") +
  geom_bar_interactive(width = 0.95, stat = "identity", show.legend = FALSE,
                       aes(tooltip = StockName, data_id = StockName))
  theme(legend.position = "none")
  
ggiraph(code=print(p))
  

NEcatchtt <- NEcatch %>%
  filter(Units=="Thousand Metric Tons") 


```

Units of `Catch` in the Northeast Shelf regional ecosystem, most recent assessments only:

`r unique(NEcatch$Units)`

Which species have catch in thousands of metric tons:

`r unique(NEcatchtt$StockName)`

These were converted to tons in the plot above.