Pull from stocksmart for 2022 SOE reports

Andy renamed the assessmentdata package stocksmart based on Stock SMART.

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

In stockAssessmentData we have time series. Columns are StockName, Year, Value, Metric, Description, Units, AssessmentYear, Jurisdiction, FMP, CommonName, ScientificName, ITIS, UpdateType, StockArea, RegionalEcosystem and the reported metrics are Catch, Fmort, Recruitment, Abundance, Index.

datatable(head(stockAssessmentData), rownames = FALSE)

In stockAssessmentSummary we have assessment metadata. Columns are Stock Name, Jurisdiction, FMP, Science Center, Regional Ecosystem, FSSI Stock?, ITIS Taxon Serial Number, Scientific Name, Common Name, Stock Area, Assessment Year, Assessment Month, Last Data Year, Update Type, Review Result, Assessment Model, Model Version, Lead Lab, Citation, Final Assessment Report 1, Final Assessment Report 2, Point of Contact, Life History Data, Abundance Data, Catch Data, Assessment Level, Assessment Frequency, Assessment Type, Model Category, Catch Input Data, Abundance Input Data, Biological Input Data, Ecosystem Linkage, Composition Input Data, F Year, Estimated F, F Unit, F Basis, Flimit, Flimit Basis, Fmsy, Fmsy Basis, F/Flimit, F/Fmsy, Ftarget, Ftarget Basis, F/Ftarget, B Year, Estimated B, B Unit, B Basis, Blimit, Blimit Basis, Bmsy, Bmsy Basis, B/Blimit, B/Bmsy, MSY, MSY Unit.

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:

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

StockSMART 2021 source, Mid-Atlantic

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

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

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

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:

Northeast Shelf, Alaska Ecosystem Complex, Atlantic Highly Migratory, Pacific Highly Migratory, Pacific Islands Ecosystem Complex, California Current, Gulf of Mexico, Southeast Shelf, Southeast Shelf / Gulf of Mexico, California Current / Alaska Ecosystem Complex, Caribbean Sea, NA

But the issue is converting to common units before summing.

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")
## List of 1
##  $ legend.position: chr "none"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE
ggiraph(code=print(p))
NEcatchtt <- NEcatch %>%
  filter(Units=="Thousand Metric Tons")

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

Metric Tons, Thousand Metric Tons

Which species have catch in thousands of metric tons:

Ocean pout - Northwestern Atlantic Coast, Windowpane - Southern New England / Mid-Atlantic

These were converted to tons in the plot above.