SOE-MAFMC-2025.Rmd

---
title: 
csl: plos.csl
fontsize: 10pt
output:
  bookdown::word_document2:
    toc: false
  bookdown::html_document2: 
    df_print: paged
  bookdown::pdf_document2:
    includes:
      in_header: latex/header1.tex
    keep_tex: yes
    toc: false
    number_sections: false
link-citations: yes
geometry: left=2cm, right=2cm, top=2cm, bottom=3cm, footskip = .5cm
subparagraph: yes
bibliography: SOE2024.bib
urlcolor: blue
always_allow_html: true
---

```{r setup, include=FALSE}

# library(tint)
# # invalidate cache when the package version changes
# knitr::opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tint'))
# options(htmltools.dir.version = FALSE)

#Default Rmd options
knitr::opts_chunk$set(echo = FALSE,
                      message = FALSE,
                      dev = "cairo_pdf",
                      warning = FALSE,
                      fig.width = 4,
                      fig.asp = 0.45,
                      fig.align = 'center'
                      ) #allows for inserting R code into captions

#Plotting and data libraries
#remotes::install_github("noaa-edab/ecodata@0.1.0") #change to 2020 ecodata version for release
library(tidyverse)
library(ecodata)
library(here)
library(flextable)
library(patchwork)

set_flextable_defaults(
  font.size = 9,
  padding.bottom = 0,  
  padding.top = 1, 
  line_spacing = 1,
  font.family = "Palatino")

```

\setcounter{page}{4}
\thispagestyle{fancy}

# Introduction

## About This Report

This report is for the Mid-Atlantic Fishery Management Council (MAFMC). The purpose of this report is to synthesize ecosystem information to allow the MAFMC to better meet fishery management objectives, and to update the MAFMC’s Ecosystem Approach to Fishery Management (EAFM) risk assessment. The major messages of the report are synthesized on pages 1 and 2, with highlights of 2024 ecosystem events on page 3. The information in this report is organized into two main sections; [performance measured against ecosystem-level management objectives](#performance-relative-to-fishery-management-objectives) (Table \ref{tab:management-objectives}), and potential [risks to meeting fishery management objectives](#risks-to-meeting-fishery-management-objectives) ((Table \ref{tab:management-risks}): [climate change](#climate-and-ecosystem-change) and [other ocean uses](#other-ocean-uses-offshore-wind)). A final section highlights [notable 2024 ecosystem observations](#highlights). 

## Report structure

A glossary of terms^[https://noaa-edab.github.io/tech-doc/glossary.html], detailed technical methods documentation^[https://NOAA-EDAB.github.io/tech-doc], indicator data^[https://noaa-edab.github.io/ecodata/], and detailed indicator descriptions^[https://noaa-edab.github.io/catalog/index.html]  are available online. We recommend new readers first review the details of standard figure formatting (Fig. \ref{fig:docformat}a), categorization of fish and invertebrate species into feeding guilds (Table \ref{tab:species-groupings}), and definitions of ecological production units (EPUs, including the Mid-Atlantic Bight, MAB; Fig. \ref{fig:docformat}b) provided at the end of the document. 

The two main sections contain subsections for each management objective or potential risk. Within each subsection, we first review observed trends for indicators representing each objective or risk, including the status of the most recent data year relative to a threshold (if available) or relative to the long-term average. Second, we identify potential drivers of observed trends, and synthesize results of indicators related to those drivers to outline potential implications for management. For example, if there are multiple drivers related to an indicator trend, do indicators associated with the drivers have similar trends, and can any drivers be affected by management action(s)? We emphasize that these implications are intended to represent testable hypotheses at present, rather than “answers,” because the science behind these indicators and syntheses continues to develop.


```{r management-objectives, ft.arraystretch = 1}

mng_obj <- data.frame("Objective Categories" = c("Seafood Production",
                                                 "Commercial Profits","Recreational Opportunities",
                                                 "Stability","Social & Cultural",
                                                 "Protected Species",
                                                 "Management","Biomass"
                                                 #,"Productivity",
                                                 #"Trophic structure","Habitat"
                                                 ),
"Indicators reported" = c("Landings; commercial total and by feeding guild; recreational harvest",
                               "Revenue decomposed to price and volume",
                               "Angler trips; recreational fleet diversity",
                               "Diversity indices (fishery and ecosystem)",
                               "Community fishing engagement and social vulnerability status",
                               "Bycatch; population (adult and juvenile) numbers; mortalities",
                               "Stock status; catch compared with catch limits",
                               "Biomass or abundance by feeding guild from surveys"#,
                               #"Condition and recruitment of managed species, primary productivity",
                               #"Relative biomass of feeding guilds, zooplankton",
                               #"Estuarine and offshore habitat conditions"
                          ))

# knitr::kable(mng_obj, linesep = "",
#       col.names = c("Objective Categories","Indicators reported"),
#       caption = "Ecosystem-scale fishery management objectives in the Mid-Atlantic Bight",
#       #align = 'c',
#       booktabs = T) %>%
#   kable_styling(latex_options = "hold_position", "scale_down") %>%
#  # column_spec(c(2), width = c("25em")) %>%
#   row_spec(0, bold = TRUE) %>%
#  # group_rows("Provisioning/Cultural", 1,4) %>%
#  # group_rows("Supporting/Regulating", 5,9)
#   pack_rows("Provisioning and Cultural Services", 1,6) %>%
#   pack_rows("Supporting and Regulating Services", 7,10)

mng_obj$service <- c(rep("Objectives: Provisioning and Cultural Services", 6), 
                     rep("Potential Drivers: Supporting and Regulating Services", 2))

flextable::as_grouped_data(mng_obj, groups = "service") %>% 
  flextable::as_flextable(hide_grouplabel=TRUE) %>% 
  flextable::align(i = ~ !is.na(service), align = "left") %>% 
  flextable::bold(i = ~ !is.na(service), bold = TRUE) %>%
  flextable::set_header_labels(Objective.Categories = "Objective categories",
                               Indicators.reported = "Indicators reported") %>%
  flextable::set_caption("Ecosystem-scale fishery management objectives in the Mid-Atlantic Bight") %>%
  flextable::autofit() 
  


```
 
```{r management-risks, ft.arraystretch = 1}

FitFlextableToPage <- function(ft, pgwidth = 6){

  ft_out <- ft %>% autofit()

  ft_out <- width(ft_out, width = dim(ft_out)$widths*pgwidth /(flextable_dim(ft_out)$widths))
  return(ft_out)
}

mng_obj <- data.frame("Risks" = c("Risks to Managing Spatially",
                                  "Risks to Managing Seasonally",
                                  "Risks to Setting Catch Limits",
                                  "Offshore Wind Risks"),
                      "Observation Indicators" = c("Managed species (fish and cetacean) distribution shifts",
                                                   "Managed species spawning and migration timing changes", 
                                                   "Managed species body condition and recruitment changes",
                                                   "Fishery revenue and landings from wind lease areas by species and port"),
                      "Potential Driver Indicators" = c("Benthic and pelagic forage distribution; ocean temperature, changes in currents and cold pool",
                                                        "Habitat timing: Length of ocean summer, cold pool seasonal persistence",
                                                        "Benthic and pelagic forage quality & abundance: ocean temperature & acidification ",
                                                        "Wind development speed; Protected species presence and  hotspots"
                      ))

# knitr::kable(mng_obj, linesep = "",
#       col.names = c("Objective Categories","Indicators reported"),
#       caption = "Ecosystem-scale fishery management objectives in the Mid-Atlantic Bight",
#       #align = 'c',
#       booktabs = T) %>%
#   kable_styling(latex_options = "hold_position", "scale_down") %>%
#  # column_spec(c(2), width = c("25em")) %>%
#   row_spec(0, bold = TRUE) %>%
#  # group_rows("Provisioning/Cultural", 1,4) %>%
#  # group_rows("Supporting/Regulating", 5,9)
#   pack_rows("Provisioning and Cultural Services", 1,6) %>%
#   pack_rows("Supporting and Regulating Services", 7,10)

mng_obj$service <- c(rep("Climate and Ecosystem Risks", 3), 
                     rep("Other Ocean Uses Risks", 1))

tab <- flextable::as_grouped_data(mng_obj, groups = "service") %>% 
  flextable::as_flextable(hide_grouplabel=TRUE) %>% 
  flextable::align(i = ~ !is.na(service), align = "left") %>% 
  flextable::bold(i = ~ !is.na(service), bold = TRUE) %>%
  flextable::set_header_labels(Risks = "Risk categories",
                               Observation.Indicators = "Observation indicators reported",
                               Potential.Driver.Indicators = "Potential driver indicators reported") %>%
  flextable::set_caption("Risks to meeting fishery management objectives in the Mid-Atlantic Bight") #%>%
  #flextable::autofit() 
  
FitFlextableToPage(tab)


```

 
# Performance Relative to Fishery Management Objectives

In this section, we examine indicators related to broad, ecosystem-level fishery management objectives. We also provide hypotheses on the implications of these trends—why we are seeing them, what’s driving them, and potential or observed regime shifts or changes in ecosystem structure. Identifying multiple drivers, regime shifts, and potential changes to ecosystem structure, as well as identifying the most vulnerable resources, can help managers determine whether anything needs to be done differently to meet objectives and how to prioritize upcoming issues/risks.

## Seafood Production 

### Indicators: Landings; commercial and recreational 
This year, we present updated indicators for total [commercial landings](https://noaa-edab.github.io/catalog/comdat.html), U.S. seafood landings, and Council-managed U.S. seafood landings. Total commercial landings within the Mid-Atlantic have declined over the long term, and total U.S. Mid-Atlantic seafood landings are near their all time low (Fig. \ref{fig:total-landings}).

```{r, total-landings, fig.asp=.8, fig.cap = "Total commercial landings (black), total U.S. seafood landings (blue), and Mid-Atlantic managed U.S. seafood landings (red), with significant decline (purple) in total landings."}
#code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-comdat-total-landings.R")
ecodata::plot_comdat(report = "MidAtlantic", varName = "landings", n=10) +
  ggplot2::ylab("")
```

Commercial landings by guild include all species and all uses, and are reported as total for the guild and the MAFMC managed species within the [guild](https://noaa-edab.github.io/catalog/species_groupings.html). Landings of benthos have been below the long term average since 2010, primarily driven by surf clam and ocean quahog, with scallops now contributing to the decline as well. Total landings of planktivores is presenting a significant downward trend, primarily due to decreases in species not managed by the MAFMC (Atlantic herring and Atlantic menhaden; Fig. \ref{fig:comm-landings}).

```{r, comm-landings, fig.asp=1.6, fig.cap ="Total commercial landings in the Mid-Atlantic Bight (black) and MAFMC-managed U.S seafood landings (red) by feeding guild, with significant declines (purple) in total planktivore landings."}
#code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-comdat-commercial-landings.R"),

ecodata::plot_comdat(report = "MidAtlantic", plottype = "guild", n=10)

```

Although total [recreational harvest](https://noaa-edab.github.io/catalog/recdat.html) (retained fish presumed to be eaten) has increased from a historic low in 2018, there is a long-term decline in the MAB (Fig. \ref{fig:rec-landings}). 

```{r, rec-landings, fig.cap ="Total recreational seafood harvest (millions of pounds, black, significant decrease, purple) in the Mid-Atlantic region."}
#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-recdat-landings.R")

ecodata::plot_recdat(report = "MidAtlantic", varName ="landings", n=10)
```

[Recreational shark landings](https://noaa-edab.github.io/catalog/rec_hms.html) show an increase in pelagic sharks over the past decade, with a sharp decrease in 2018 - 2019 persisting through 2022 (Fig \ref{fig:rec-hms}). This is likely influenced by regulatory changes implemented in 2018 intended to rebuild shortfin mako stocks. 

```{r rec-hms, fig.cap="Recreational shark landings from Large Pelagics Survey.",  fig.width=5, fig.asp=.35}
#code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-rec_hms.R"),

ecodata::plot_rec_hms(shadedRegion = c(2014,2023), report="MidAtlantic")  +
  ggplot2::scale_color_discrete(        
    limits = c("LargeCoastal", "Pelagic", "Prohibited", "SmallCoastal"),         
    labels = c("Large Coastal", "Pelagic", "Prohibited", "Small Coastal")
    ) 
```

Aquaculture production is not yet included in total seafood landings. Available [aquaculture production](https://noaa-edab.github.io/catalog/aquaculture.html) of oysters for a subset of Mid-Atlantic states indicates a decline in recent years.

### Implications

Declining commercial (total and seafood) and recreational landings can be driven by many interacting factors, including combinations of ecosystem and stock production, management actions, market conditions, and environmental change. While we cannot evaluate all possible drivers at present, here we evaluate the extent to which stock status and system biomass trends may play a role.

#### Stock Status and Catch Limits
Single species [management objectives](https://noaa-edab.github.io/catalog/stock_status.html) (1. maintaining biomass above minimum thresholds and 2. maintaining fishing mortality below overfishing limits) are being met for all but three MAFMC-managed species (Fig. \ref{fig:stock-status}), though the status of six stocks is unknown (Table \ref{tab:unkstocks}). 

```{r stock-status, fig.width = 7.5,  fig.asp = 0.5, fig.cap = "Summary of single species status for MAFMC and jointly federally managed stocks (Spiny dogfish and both Goosefish). The dotted vertical line is the target biomass reference point of $B_{MSY}$. The dashed lines are the management thresholds of one half $B_{MSY}$ (vertical) or $F_{MSY}$. (horizontal). Stocks in orange are below the biomass threshold (overfished) or have fishing mortality above the limit (subject to overfishing), so are not meeting objectives. Stocks in purple  are above the biomass threshold but below the biomass target with fishing mortality within the limit. Stocks in green are above the biomass target, with fishing mortality within the limit."}

# code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-stock-status.R"), 

a <- ecodata::plot_stock_status(report = "MidAtlantic")

a$p + ggplot2::coord_cartesian(xlim=c(0,2), ylim=c(0,2)) 

#+
#  ggplot2::annotation_custom(gridExtra::tableGrob(a$unknown,
#                             theme = gridExtra::ttheme_default(base_size = 7), 
#                             rows=NULL), 
#                             xmin=0.8, xmax=1.8, ymin=1.5, ymax=2)

```
**ADD BLUELINE TILEFISH, CHUB MACKEREL**
```{r unkstocks}

flextable::flextable(a$unknown) |>
  flextable::set_header_labels(F.Fmsy = "F/Fmsy",
                               B.Bmsy = "B/Bmsy") |>
  flextable::colformat_num(na_str = "-") |>
  flextable::set_caption("Unknown or partially known stock status for MAFMC and jointly managed species.") |>
  flextable::autofit() 


```


Stock status affects catch limits established by the Council, which in turn may affect landings trends. Summed across all MAFMC managed species, total Acceptable Biological Catch or Annual Catch Limits [(ABC or ACL)](https://noaa-edab.github.io/catalog/abc_acl.html) have been relatively stable 2012-2022 (Fig. \ref{fig:abcacl-stacked}). The recent total ABC or ACL is lower relative to 2012-2013, with much of that decrease due to declining Atlantic mackerel ABC. This is true even with the addition of blueline tilefish management in 2017 contributing an additional ABC and ACL to the total 2017-2022, due to that fishery’s small relative size.

```{r abcacl-stacked, fig.cap="Sum of catch limits across all MAFMC managed commercial (C) and recreational (R) fisheries.", fig.width = 6, fig.asp = 0.4}
#, code =  readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-abcacl-stacked.R")

ecodata::plot_abc_acl(report = "MidAtlantic", plottype = "Stacked")

```

Nevertheless, the percentage caught (landings and discards) for each stock’s ABC/ACL suggests that these catch limits are not generally constraining as most species are well below the 1/1 ratio (Fig. \ref{fig:abcacl-catch}). Therefore, stock status and associated management constraints are unlikely to be driving decreased landings for the majority of species.

```{r abcacl-catch, fig.cap = "Catch divided by ABC/ACL for MAFMC managed fisheries. High points are recreational black sea bass (up to 2021) and scup (2022). Red line indicates the median ratio across all fisheries."}
#, code =  readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-abcacl-catch.R")

ecodata::plot_abc_acl(report = "MidAtlantic", plottype = "Catch")

```

#### System Biomass
Although [aggregate biomass](https://noaa-edab.github.io/catalog/aggregate_biomass.html) trends derived from scientific resource surveys are mostly stable in the MAB, spring piscivores, fall benthivores, and fall benthos show long-term increases (Fig. \ref{fig:nefsc-biomass-mab}). While managed species make up varying proportions of aggregate biomass, trends in landings are not mirroring shifts in the overall trophic structure of survey-sampled fish and invertebrates. Therefore, major shifts in feeding guilds or ecosystem trophic structure are unlikely to be driving the decline in landings.
```{r nefsc-biomass-mab, fig.cap = "Spring (left) and fall (right) surveyed biomass in the Mid-Atlantic Bight. Data from the NEFSC Bottom Trawl Survey are shown in black, with the nearshore NEAMAP survey shown in red. Significant increases (orange lines) are present for spring piscivore and fall benthivore and benthos biomass. The shaded area around each annual mean represents 2 standard deviations from the mean.", fig.width=8, fig.asp = 0.75, results='hide'}

#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/macrofauna_MAB.Rmd-aggregate-biomass.R")

ecodata::plot_aggregate_biomass(report = "MidAtlantic", n=10) 

```

#### Effect on Seafood Production
Stock status is above the minimum threshold for all but three stocks, and aggregate biomass trends appear stable, so the decline in managed commercial seafood landings is most likely driven by market dynamics affecting the landings of surfclams and ocean quahogs, as landings have been below quotas for these species. In addition, regional availability of scallops has contributed to the decline of benthos landings not managed by the MAFMC, with some of the most productive grounds currently closed through rotational management. The long term decline in total planktivore landings is largely driven by Atlantic menhaden fishery dynamics, including a consolidation of processors leading to reduced fishing capacity between the 1990s and mid-2000s. 

Climate change also seems to be shifting the distribution of surfclams and ocean quahogs, resulting in areas with overlapping distributions and increased mixed landings. Given the regulations governing mixed landings, this could become problematic in the future and is currently being evaluated by the Council.

The decline in recreational seafood harvest stems from other drivers. Some of the decline, such as that for recreational shark landings, is driven by management intended to reduce fishing mortality on mako sharks. However, NOAA Fisheries’ Marine Recreational Information Program survey methodology was updated in 2018, so it is unclear whether the lower than average landings for species other than sharks since 2018 are driven by changes in fishing behavior or the change in the survey methodology. Nevertheless, the recreational harvest seems to be stabilizing at a lower level than historical estimates.

[PLACEHOLDER AGGREGATE COMMUNITY LANDING CLIMATE VULNERABIITY]

Other environmental changes require monitoring as they may become important drivers of landings in the future:
 
- Climate is trending into uncharted territory. Globally, 2024 was the warmest year on record (see [2024 Highlights section](#highlights)). 
- Stocks are shifting their distributions, moving towards the northeast and into deeper waters throughout the Northeast US Large Marine Ecosystem (see [Climate Risks section](#climate-and-ecosystem-change)).
- Some ecosystem composition and production changes have been observed (see [Stability section](#stability)). 
- Some fishing communities are affected by socioeconomic vulnerabilities (see [Community Social and Climate Vulnerability section](#community-social-and climate-vulnerability)).  

## Commercial Profits 

### Indicators: revenue (a proxy for profits)  
Total [commercial revenue](https://noaa-edab.github.io/catalog/comdat.html) and MAFMC managed species revenue within the Mid-Atlantic Bight have declined over the past 20-30 years. In 2023, total revenue was at an all-time low, and revenue from MAFMC managed species was near an all-time low (Fig. \ref{fig:comm-revenue}). 
```{r comm-revenue, fig.width = 6, fig.asp = 0.45, fig.cap = "Revenue for the for the Mid-Atlantic region: total (black) and from MAFMC managed species (red)."}
#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-comdat-comm-revenue.R")

ecodata::plot_comdat(report="MidAtlantic", varName="revenue", n=10) + 
  ggplot2::theme(legend.position = "right",
                      legend.title = ggplot2::element_blank())

```

Revenue earned by harvesting resources is a function of both the quantity landed of each species and the prices paid for landings. Beyond monitoring yearly changes in revenue, it is even more valuable to determine what drives these changes: harvest levels, the mix of species landed, price changes, or a combination of these. The [Bennet Indicator](https://noaa-edab.github.io/catalog/bennet.html) decomposes revenue change into two parts, one driven by changing quantities (volumes), and a second driven by changing prices. All changes are in relation to a base year (1982).

In the Mid-Atlantic region revenues were above the 1982 baseline for all years in the series until 2022 (Fig. \ref{fig:bennet}).  Both increasing prices and volumes contributed to the positive revenue change in most years.  In terms of prices, since 2000 Benthos contributed the most to increasing prices (Fig.  \ref{fig:bennet-all}). Beginning in the 1990s, in most years benthivores contributed the most to increasing volumes (Fig.  \ref{fig:bennet-all}). 
```{r bennet, fig.cap = "Revenue change from 1982 values in 2022 dollars (black); Price (PI), and Volume Indicators (VI) for total commercial landings in the Mid-Atlantic Bight.", fig.width=6, fig.asp=.35}
#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-bennet.R")

ecodata::plot_bennet(report = "MidAtlantic", varName = "total")

```

```{r bennet-all, fig.cap = "Total price and volume indicators in 2022 dollars (black) for commercial landings, and individual guild contributions to each indicator, in the Mid-Atlantic Bight.",  fig.width = 7.5, fig.asp = 0.35}
#code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-bennet-all.R"),

ecodata::plot_bennet(report="MidAtlantic") #+ ggplot2::facet_wrap(~Var, ncol=1)
```

### Implications 
Although the Mid-Atlantic region shows declining revenue since 2016, inflation-adjusted revenue from harvested species was still greater than 1982 levels until 2022. In a similar manner to seafood landings, the results here are driven in large part by market dynamics affecting the landings of surfclams and ocean quahogs, as landings have been below quotas for these species, as well as lower quotas for Atlantic scallops. The declining Benthos category since 2012 may be partially caused by decreases in surfclam and ocean quahogs in the southern part of their range as harvest have shifted northward. Changes in other indicators, particularly those driving landings and those related to climate change, require monitoring as they may become important drivers of revenue in the future; for example:

-  Surfclams, ocean quahogs, and scallops are sensitive to warming ocean temperatures and ocean acidification.   
-  [Multiple stressors](https://noaa-edab.github.io/catalog/ocean_acidification) are interacting in Mid-Atlantic shellfish habitats.

[PLACEHOLDER AGGREGATE COMMUNITY REVENUE CLIMATE VULNERABIITY]
  

## Recreational Opportunities 

### Indicators: Angler trips, fleet diversity  
[Recreational effort](https://noaa-edab.github.io/catalog/recdat.html) (angler trips) in 2022 has increased and is above the long-term average (Fig. \ref{fig:rec-op}). in the MAB. However, recreational fleet diversity (i.e., effort by shoreside, private boat, and for-hire anglers) has declined over the long term (Fig. \ref{fig:rec-div}). 
```{r rec-op, fig.cap = "Recreational effort (number of trips, black) in the Mid-Atlantic, with significant increase (orange line)."}
# code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-recdat-effort.R")

ecodata::plot_recdat(report = "MidAtlantic", varName = "effort", n=10)
```

```{r rec-div, fig.cap = "Recreational fleet effort diversity (black) in the Mid-Atlantic, with significant decrease (purple line)."}
#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-recdat-diversity.R")

ecodata::plot_recdat(report = "MidAtlantic", varName = "effortdiversity", n=10)
```

### Implications
While the overall number of recreational opportunities in the MAB is above the long-term average, the continuing decline in recreational fleet effort diversity suggests a potentially reduced range of recreational fishing options, despite the slight increase in this indicator’s value between 2020 and 2022.

The downward effort diversity trend is driven by party/charter contraction (down from 2% in 2021 to 1.4% in 2023), and a shift toward shorebased angling, which currently makes up 59% of angler trips. Effort in private boats has increased slightly to 40% of trips from 37% in 2021.

Changes in recreational fleet diversity can be considered when managers seek options to maintain recreational opportunities. Shore anglers will have access to different species than vessel-based anglers, and when the same species is accessible both from shore and from a vessel, shore anglers typically have access to smaller individuals. Many states have developed shore-based regulations where the minimum size is lower than in other areas and sectors to maintain opportunities in the shore angling sector.


## Stability 

### Indicators: fishery fleet and catch diversity, ecological component diversity 
While there are many potential metrics of stability, we use diversity indices to evaluate overall stability in fisheries and ecosystems. In general, diversity that remains constant over time suggests a similar capacity to respond to change over time. A significant change in diversity over time does not necessarily indicate a problem or an improvement, but does indicate a need for further investigation. We examine diversity in commercial fleet and species catch, recreational species catch (with fleet effort diversity discussed above), zooplankton, adult fishes, and fish traits.

#### Fishery Stability 
[Diversity](https://noaa-edab.github.io/catalog/commercial_div.html)  estimates have been developed for fleets landing managed species, and species landed by commercial vessels with Mid-Atlantic permits. Commercial fishery fleet count has declined while fleet diversity has been stable over time in the MAB, with current values near the long-term average. 

```{r, fig.width = 7.5, fig.asp = 0.3}
a <- ecodata::plot_commercial_div(report = "MidAtlantic", varName = "Fleet count", n=21)

b <- ecodata::plot_commercial_div(report = "MidAtlantic", varName = "Fleet diversity in revenue", n=21)

a+b

```


This indicates different commercial fleet composition but similar species targeting opportunities over time. Commercial fisheries are relying on fewer species relative to the mid-90s, although current species revenue diversity has recovered in the last year (Fig. \ref{fig:commercial-div-species-div}). 

CAVEAT ABOUT CLAM DATA?

```{r commercial-div-species-div,  fig.cap = paste0("Species revenue diversity in the Mid Atlantic.")}
#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-commercial-div-species-div.R")

ecodata::plot_commercial_div(report = "MidAtlantic", varName = "Permit revenue species diversity", n=21)
```

As noted [above](#recreational-opportunities), [recreational fleet effort diversity](https://noaa-edab.github.io/catalog/recdat.html) is declining (Fig. \ref{fig:rec-div}), suggesting a shift in recreational fishing opportunities. However, recreational species catch diversity has no long term trend so is considered stable, and has been at or above the long term average in 8 of the last 10 years (Fig. \ref{fig:recdat-div-catch}). 

```{r recdat-div-catch, fig.cap = paste0("Diversity of recreational catch in the Mid Atlantic.")}
#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-recdat-div-catch.R")

ecodata::plot_recdat(report = "MidAtlantic", varName = "catchdiversity", n=10)
```

#### Ecological Stability

Total [primary production](https://noaa-edab.github.io/catalog/chl_pp.html) (PP) is a measure of the total energy input into a system per year. Total primary production in the Mid Atlantic has no clear trend (Fig. \ref{fig:totpp}), suggesting stability in energy at the base of the food web. 

```{r totpp, fig.cap="Total areal annual primary production for the MAB. The dashed line represents the long-term (1998-2023) annual mean."}

ecodata::plot_annual_chl_pp(report = "MidAtlantic", varName = "pp", plottype = "total") +
  ggplot2::coord_cartesian(ylim = c(1e+07, 4e+07), xlim = c(1998, 2023)) + 
  ggplot2::ggtitle("MAB Primary Production") +
  ggplot2::ylab("Carbon (mt)")
```

Ecological diversity indices show mixed trends. [Zooplankton diversity](https://noaa-edab.github.io/catalog/zoo_diversity.html) is increasing in the MAB (Fig. \ref{fig:zoo-diversity}). [Adult fish diversity](https://noaa-edab.github.io/catalog/exp_n.html) is measured as the expected number of species in a standard number of individuals sampled from the NEFSC bottom trawl survey. Adult fish diversity indices appear stable over time, with current values within one standard deviation from most historic estimates (Fig. \ref{fig:exp-n}). 

```{r zoo-diversity, fig.cap = "Zooplankton diversity in the Mid-Atlantic Bight, Shannon diversity index (black) with significant increase (orange line)."}
#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/LTL_MAB.Rmd-zoo-diversity.R")

ecodata::plot_zoo_diversity(report = "MidAtlantic", n=10)
```


```{r exp-n, fig.cap = "Adult fish diversity in the Mid-Atlantic Bight, based on expected number of species. Results from survey vessels Albatross and Bigelow are reported separately due to catchability differences.", fig.width=5, fig.asp=.4}
#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/macrofauna_MAB.Rmd-exp-n.R")

ecodata::plot_exp_n(report = "MidAtlantic", varName = "fall", n=10)
```

Functional traits, such as length at maturity, asymptotic body size, or fecundity, can synthesize change across complex, diverse communities. Monitoring changes in functional trait distributions can provide a means of assessing ecosystem-scale resilience. There was evidence for long term change in trait distributions in the mid-Atlantic Bight. The spring finfish community in the mid-Atlantic Bight is showing long-term shifts towards slower life history strategies with lower fecundity.

[PLACEHOLDER FIG]


### Implications
Fleet diversity indices are used by the MAFMC to evaluate stability objectives as well as risks to fishery resilience and maintaining equity in access to fishery resources. Instability in one commercial fleet diversity metric suggests potentially lower capacity to respond to the current range of fishing opportunities. Commercial species diversity is relatively stable but comparisons are limited by missing clam fishery data prior to 2003.

Declining recreational fleet effort diversity, as noted above, indicates that the party/charter boat sector continues to contract, with shoreside angling becoming more important as a percentage of recreational angler trips. Stability in recreational species catch diversity has been maintained by a different set of species over time. A recent increase in Atlantic States Marine Fisheries Commission (ASMFC) and South Atlantic Fishery Management Council (SAFMC) managed species in recreational catch is helping to maintain diversity in the same range that MAFMC and New England Fishery Management Council (NEFMC) managed species supported in the 1990s.

Production at the base of the food web is variable, but stable over time. Ecological diversity indices can provide insight into ecosystem structure. Changes in ecological diversity over time may indicate altered ecosystem structure with implications for fishery productivity and management. Stable adult fish diversity indicates the same overall number and evenness over time, but doesn’t rule out species substitutions (e.g., warm-water replacing cold-water). In addition, the change in survey vessels complicates interpretation of long-term fish diversity trends.

In the MAB, existing diversity indicators suggest some instability in the fisheries and ecosystem components examined. In addition, declining recreational fleet diversity suggests a potential loss in the range of recreational fishing opportunities. Increasing zooplankton diversity (due to increases in abundance of several taxa and stable or declining dominance of an important copepod species) suggests a shift in the zooplankton community that warrants continued monitoring to determine if managed species are affected. In addition, the species diversity in commercial landings warrants continued attention given its relatively low index value indicating average reliance on a small number of species for revenue.

## Community Social and Climate Vulnerability 

Providing for sustained participation of fishing communities, and avoiding adverse economic impacts to fishing communities are objectives of fishery management. We report the top ten communities most engaged in, and/or reliant upon, commercial and recreational fisheries and the degree to which these communities may be vulnerable to environmental justice issues (i.e., Poverty, Population Composition, and Personal Disruption) using data for the most recent available year (2022). We also compare these results with those presented in previous SOE reports to highlight changes in community status. 

### Indicators: Fishing Engagement and Community Social Vulnerability

The [engagement](https://noaa-edab.github.io/catalog/engagement.html) indices demonstrate the importance of commercial and recreational fishing to a given community relative to other coastal communities in a region. Social vulnerability indicators measure social factors that shape a community’s ability to adapt to change. A subset of these factors can be used to assess potential environmental justice issues. Similarly to the engagement and reliance indicators, the environmental justice indices characterize different facets and levels of social vulnerability in a given community relative to other coastal communities in a region. 

Changes in fishing activity between years changed community engagement and reliance rankings. The largest change from last year’s report is that Hatteras and Hobucken, NC are no longer listed as top ten commercial fishing communities, replaced by Hampton, VA; Swan Quarter, NC; Bowers and Little Creek, DE (Fig.\ref{fig:commercial-engagement}).  Manteo, Vandemere, and Hobuken, NC are no longer listed as top ten recreational communities, replaced by Cape May and Barnegat Light, NJ; Orient, NY; Topsail Beach, Avon and Rodanthe, NC (Fig.\ref{fig:recreational-engagement}). 

 
```{r commercial-engagement, fig.cap= "Commercial engagement, reliance, and environmental justice vulnerability for the top commercially engaged and reliant fishing communities in the Mid-Atlantic.", fig.width = 7, fig.asp = 0.85, results='hide'}
#Communities ranked medium-high or above for one or more of the environmental justice indicators are highlighted in bright orange. *Community scored high (1.00 and above) for both commercial engagement and reliance indicators.
#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-commercial-engagement.R")
#fig.width = 6.5, fig.asp = 0.75
#(* Scored high (1.00 and above)) for both commercial engagement and reliance indicators)

#ecodata::plot_engagement(report = "MidAtlantic", varName = "Commercial") 

dat <- read.csv(here::here("data-raw/fishdata_2022_CSVI_1-23-2025.csv"))

dat$GEO_NAME <- replace(dat$GEO_NAME,
        dat$GEO_NAME=="Port Clyde-Tenants Harbor/Saint George/Spruce Head, ME",
        "Port Clyde-Tenants Harbor, ME")


dat$GEO_NAME <- replace(dat$GEO_NAME,
                        dat$GEO_NAME=="Sandwich/East Sandwich/Forestdale, MA",
                        "Sandwich, MA")


dat$GEO_NAME <- replace(dat$GEO_NAME,
                        dat$GEO_NAME=="Reedville/District 5 (Northumberland County), VA",
                        "Reedville, VA*")

##reorder ranks
dat$ComEng_ct <- 
  factor(dat$ComEng_ct, levels=c("high", "med high", "med", "low"))
dat$ComRel_ct <- 
  factor(dat$ComRel_ct, levels=c("high", "med high", "med", "low"))
dat$RecEng_ct <- 
  factor(dat$RecEng_ct, levels=c("high", "med high", "med", "low"))
dat$RecRel_ct <- 
  factor(dat$RecRel_ct, levels=c("high", "med high", "med", "low"))
dat$poverty_rank <- 
  factor(dat$poverty_rank, levels=c("high", "med high", "med", "low"))
dat$pop_composition_rank <- 
  factor(dat$pop_composition_rank, levels=c("high", "med high", "med", "low"))
dat$personal_disruption_rank <- 
  factor(dat$personal_disruption_rank, levels=c("high", "med high", "med", "low"))

#make data numeric
dat$ComEng <- as.numeric(dat$ComEng)
dat$ComRel <- as.numeric(dat$ComRel)
dat$RecEng <- as.numeric(dat$RecEng)
dat$RecRel <- as.numeric(dat$RecRel)

#remove NA
dat_com_e<-dat[!is.na(dat$ComEng),]
dat_com_r<-dat[!is.na(dat$ComRel),]
dat_rec_e<-dat[!is.na(dat$RecEng),]
dat_rec_r<-dat[!is.na(dat$RecRel),]

#remove low scores
dat_com_e_nolow<-dat_com_e[!(dat_com_e$ComEng_ct=='low'),]
dat_com_r_nolow<-dat_com_r[!(dat_com_r$ComRel_ct=='low'),]
dat_rec_e_nolow<-dat_rec_e[!(dat_rec_e$RecEng_ct=='low'),]
dat_rec_r_nolow<-dat_rec_r[!(dat_rec_r$RecRel_ct=='low'),]

#subset Mid Atlantic

dat_com_e_MA <- dat_com_e[(dat_com_e$STATEABBR=='NY' | 
                             dat_com_e$STATEABBR=='PA'| 
                             dat_com_e$STATEABBR=='NJ'| 
                             dat_com_e$STATEABBR=='DE'| 
                             dat_com_e$STATEABBR=='MD'| 
                             dat_com_e$STATEABBR=='VA'| 
                             dat_com_e$STATEABBR=='NC'), ]

dat_com_e_nolow_MA <- dat_com_e_nolow[(dat_com_e_nolow$STATEABBR=='NY' | 
                                         dat_com_e_nolow$STATEABBR=='PA'| 
                                         dat_com_e_nolow$STATEABBR=='NJ'| 
                                         dat_com_e_nolow$STATEABBR=='DE'| 
                                         dat_com_e_nolow$STATEABBR=='MD'| 
                                         dat_com_e_nolow$STATEABBR=='VA'| 
                                         dat_com_e_nolow$STATEABBR=='NC'),]

dat_com_r_MA <- dat_com_r[(dat_com_r$STATEABBR=='NY' | 
                             dat_com_r$STATEABBR=='PA'| 
                             dat_com_r$STATEABBR=='NJ'| 
                             dat_com_r$STATEABBR=='DE'| 
                             dat_com_r$STATEABBR=='MD'| 
                             dat_com_r$STATEABBR=='VA'| 
                             dat_com_r$STATEABBR=='NC'),]

dat_com_r_nolow_MA <- dat_com_r_nolow[(dat_com_r_nolow$STATEABBR=='NY' | 
                                         dat_com_r_nolow$STATEABBR=='PA'| 
                                         dat_com_r_nolow$STATEABBR=='NJ'| 
                                         dat_com_r_nolow$STATEABBR=='DE'| 
                                         dat_com_r_nolow$STATEABBR=='MD'| 
                                         dat_com_r_nolow$STATEABBR=='VA'| 
                                         dat_com_r_nolow$STATEABBR=='NC'),]


dat_rec_e_MA <- dat_rec_e[(dat_rec_e$STATEABBR=='NY' | 
                             dat_rec_e$STATEABBR=='PA'| 
                             dat_rec_e$STATEABBR=='NJ'| 
                             dat_rec_e$STATEABBR=='DE'| 
                             dat_rec_e$STATEABBR=='MD'| 
                             dat_rec_e$STATEABBR=='VA'| 
                             dat_rec_e$STATEABBR=='NC'),]

dat_rec_e_nolow_MA <- dat_rec_e_nolow[(dat_rec_e_nolow$STATEABBR=='NY' | 
                                         dat_rec_e_nolow$STATEABBR=='PA'| 
                                         dat_rec_e_nolow$STATEABBR=='NJ'| 
                                         dat_rec_e_nolow$STATEABBR=='DE'| 
                                         dat_rec_e_nolow$STATEABBR=='MD'| 
                                         dat_rec_e_nolow$STATEABBR=='VA'| 
                                         dat_rec_e_nolow$STATEABBR=='NC'),]

dat_rec_r_MA <- dat_rec_r[(dat_rec_r$STATEABBR=='NY' | 
                             dat_rec_r$STATEABBR=='PA'| 
                             dat_rec_r$STATEABBR=='NJ'| 
                             dat_rec_r$STATEABBR=='DE'| 
                             dat_rec_r$STATEABBR=='MD'| 
                             dat_rec_r$STATEABBR=='VA'| 
                             dat_rec_r$STATEABBR=='NC'),]

dat_rec_r_nolow_MA <- dat_rec_r_nolow[(dat_rec_r_nolow$STATEABBR=='NY' | 
                                         dat_rec_r_nolow$STATEABBR=='PA'| 
                                         dat_rec_r_nolow$STATEABBR=='NJ'| 
                                         dat_rec_r_nolow$STATEABBR=='DE'| 
                                         dat_rec_r_nolow$STATEABBR=='MD'| 
                                         dat_rec_r_nolow$STATEABBR=='VA'| 
                                         dat_rec_r_nolow$STATEABBR=='NC'),]



# logic: category gets the max indicator rank, I think this is how EJ works
# ignore this not using it right now
commsocvulrank <- dat |>
  #dplyr::filter(REGION=="Northeast") |>
  dplyr::mutate(across(25:32, ~ case_when(   #stringr::str_detect("_rank")
                  . == "low" ~ "1",
                  . == "med" ~ "2",
                  . == "med high" ~ "3",
                  . == "high" ~ "4",
                  TRUE ~ NA_character_
                ))) |>
  dplyr::rowwise() |>
  dplyr::mutate(EJrank = max(as.numeric(personal_disruption_rank), 
                             as.numeric(pop_composition_rank), 
                             as.numeric(poverty_rank)),
                Econrank = max(as.numeric(labor_force_str_rank),
                               as.numeric(housing_characteristics_rank)),
                Gentrank = max(as.numeric(housing_disrupt_rank),
                               as.numeric(retiree_migration_rank),
                               as.numeric(urban_sprawl_index_rank)))


MA_com_plot <- dat |>
  dplyr::filter(STATEABBR %in% c("DE", "MD", "NJ", "NY", "NC", "PA", "VA")) |>
  ggplot2::ggplot(ggplot2::aes(x=ComEng, y=ComRel)) + 
  ggplot2::geom_point(shape=1, 
             size=1.4, stroke=0.9, alpha=0.6, 
             ggplot2::aes(color=ComEng_ct))+
  ggplot2::scale_color_manual(values = c("high" = "navy",
                                "med high"="slateblue3",
                                "med"="skyblue",
                                "low"="gray"))+
             #aes(color=as.factor(EJrank)))+
  #ggplot2::scale_colour_viridis_d(direction = -1) +
  # ggplot2::scale_color_manual(values = c(4 ="firebrick4",
  #                               3="darkorange1",
  #                               2="goldenrod1",
  #                               1="gray")) +
  ggplot2::theme_bw()+
  ggplot2::ylim(-3,20)+
  ggplot2::xlim(-1,10)+
  ggplot2::theme(#text = element_text(size=10),
        #plot.title = element_text(size=9, vjust=-7.5, hjust=0.01),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        legend.position = "bottom")+
  ggplot2::ylab(expression(paste("Relative Commercial Engagement Index", italic(" (factor score)"))))+
  ggplot2::xlab(expression(paste("Commercial Engagement Index", italic(" (factor score)"))))+
  ggplot2::guides(color=guide_legend(title="Engagement Score Rank"))+
  ggrepel::geom_label_repel(aes(label = ifelse(ComEng>2,as.character(GEO_NAME),ifelse(ComRel>1.28,as.character(GEO_NAME),''))),
                   size=2,
                   box.padding   = 0.4, 
                   point.padding = 0,
                   segment.color = 'grey50',
                   min.segment.length = 0.01,
                   max.overlaps=Inf,
                   label.size = NA,
                   fill = alpha(c("white"),0.1))+
  ggplot2::ggtitle("engagement index")

MA_com_poverty <- ggplot(dat_com_e_MA %>% filter(poverty_rank == "high"|
                                                   poverty_rank == "med high"|
                                                   poverty_rank =="med"), 
                         aes(x=ComEng, y=ComRel)) + 
  geom_point(shape=1, 
             size=1.4, stroke=0.9, alpha=0.6, 
             aes(color=poverty_rank, ))+
  scale_color_manual(values = c("high" = "firebrick4",
                                "med high"="darkorange1",
                                "med"="goldenrod1",
                                "low"="skyblue2")) +
  theme_bw()+
  ylim(-3,20)+
  xlim(-1,10)+
  theme(text = element_text(size=9),
        plot.title = element_text(size=8, vjust=-5, hjust=0.01),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        legend.position="none")+
  ggrepel::geom_label_repel(aes(label = ifelse(ComEng>2,as.character(GEO_NAME),ifelse(ComRel>1.28,as.character(GEO_NAME),''))),
                   size=1.6,
                   force=2.5,
                   box.padding   = 0.4, 
                   point.padding = 0,
                   segment.color = 'grey50',
                   min.segment.length = 0.01,
                   max.overlaps=Inf,
                   label.size = NA,
                   fill = alpha(c("white"),0.1))+
  ggtitle("poverty index")

MA_com_pers <- ggplot(dat_com_e_MA %>% filter(personal_disruption_rank == "high"|personal_disruption_rank == "med high"|personal_disruption_rank =="med"), aes(x=ComEng, y=ComRel)) + 
  geom_point(shape=1, 
             size=1.4, stroke=0.9, alpha=0.6, 
             aes(color=personal_disruption_rank, ))+
  scale_color_manual(values = c("high" = "firebrick4",
                                "med high"="darkorange1",
                                "med"="goldenrod1",
                                "low"="skyblue2")) +
  theme_bw()+
  ylim(-3,20)+
  xlim(-1,10)+
  theme(text = element_text(size=9),
        plot.title = element_text(size=8, vjust=-7.5, hjust=0.01),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        legend.position="none")+
  ggrepel::geom_label_repel(aes(label = ifelse(ComEng>2,as.character(GEO_NAME),ifelse(ComRel>1.28,as.character(GEO_NAME),''))),
                   size=1.6,
                   force=2.5,
                   box.padding   = 0.4, 
                   point.padding = 0,
                   segment.color = 'grey50',
                   min.segment.length = 0.01,
                   max.overlaps=Inf,
                   label.size = NA,
                   fill = alpha(c("white"),0.1))+
  ggtitle("personal disruption index")


MA_com_pop <- ggplot(dat_com_e_MA %>% filter(pop_composition_rank == "high"| 
                                               pop_composition_rank == "med high"|
                                               pop_composition_rank =="med"), 
                     aes(x=ComEng, y=ComRel)) + 
  geom_point(shape=1, 
             size=1.4, stroke=0.9, alpha=0.6, 
             aes(color=pop_composition_rank, ))+
  scale_color_manual(values = c("high" = "firebrick4",
                                "med high"="darkorange1",
                                "med"="goldenrod1",
                                "low"="skyblue2")) +
  theme_bw()+
  ylim(-3,20)+
  xlim(-1,10)+
  theme(text = element_text(size=9),
        plot.title = element_text(size=8, vjust=-7.5, hjust=0.01),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        legend.position="right")+
  guides(color=guide_legend(title="Vulnerability\nIndicator\nScore Rank"))+
  ggrepel::geom_label_repel(aes(label = ifelse(ComEng>2,as.character(GEO_NAME),ifelse(ComRel>1.28,as.character(GEO_NAME),''))),
                   size=1.6,
                   force=2.5,
                   box.padding   = 0.4, 
                   point.padding = 0,
                   segment.color = 'grey50',
                   min.segment.length = 0.01,
                   max.overlaps=Inf,
                   label.size = NA,
                   fill = alpha(c("white"),0.1))+
  ggtitle("population composition index")


MA_com_plot + (MA_com_poverty / MA_com_pers / MA_com_pop) +
  plot_layout(ncol=2, widths = c(5,2))
```

Scores for environmental justice concerns remain similar for communities based on 2020 and 2021 data, with top commercial ports (Fig. \ref{fig:commercial-EJ}) showing more concerns than top recreational ports overall (Fig. \ref{fig:recreational-EJ}).  Atlantic City, NJ ranks highest for all three environmental justice concerns. There is variability in the specific issues facing communities with environmental justice concerns.  Higher scores in population composition indicate community vulnerability related to the presence of non-white, non-English speaking, and younger populations. 

```{r recreational-engagement, fig.cap= "Recreational engagement and reliance, and environmental justice vulnerability, for the top recreationally engaged and reliant fishing communities in the Mid-Atlantic.", fig.width = 7, fig.asp = 0.85, results='hide'}
#Communities ranked medium-high or above for one or more of the environmental justice indicators are highlighted in bright orange. *Community scored high (1.00 and above) for both recreational engagement and reliance indicators.
#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-recreational-engagement.R")
#fig.width = 6.5, fig.asp = 0.75

#(* Scored high (1.00 and above) for both recreational engagement and reliance indicators;

#ecodata::plot_engagement(report = "MidAtlantic", varName = "Recreational")

####Mid-Atlantic Recreational Engagement####

MA_rec_plot <- ggplot(dat_rec_e_MA, aes(x=RecEng, y=RecRel)) + 
  geom_point(shape=1, 
             size=1.4, stroke=0.9, alpha=0.6, 
             aes(color=RecEng_ct))+
  scale_color_manual(values = c("high" = "navy",
                                "med high"="slateblue3",
                                "med"="skyblue",
                                "low"="gray")) +
  theme_bw()+
  ylim(-2,25)+
  xlim(-1,16)+
  theme(text = element_text(size=10),
        plot.title = element_text(size=9, vjust=-7.5, hjust=0.01),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        legend.position = "bottom")+
  ylab(expression(paste("Population Relative Recreational Engagement Index", italic(" (factor score)"))))+
  xlab(expression(paste("Recreational Engagement Index", italic(" (factor score)"))))+
  guides(color=guide_legend(title="Engagement Score Rank"))+
  ggrepel::geom_label_repel(aes(label = ifelse(RecEng>7,as.character(GEO_NAME),ifelse(RecRel>7,as.character(GEO_NAME),''))),
                   size=2,
                   box.padding   = 0.4, 
                   point.padding = 0,
                   segment.color = 'grey50',
                   min.segment.length = 0.01,
                   max.overlaps=Inf,
                   label.size = NA,
                   fill = alpha(c("white"),0.1))+
  ggtitle("engagement index")


MA_rec_poverty <- ggplot(dat_rec_e_MA %>% filter(poverty_rank == "high"|
                                                   poverty_rank == "med high"|
                                                   poverty_rank =="med"), aes(x=RecEng, y=RecRel)) + 
  geom_point(shape=1, 
             size=1.4, stroke=0.9, alpha=0.6, 
             aes(color=poverty_rank, ))+
  scale_color_manual(values = c("high" = "firebrick4",
                                "med high"="darkorange1",
                                "med"="goldenrod1",
                                "low"="skyblue2")) +
  theme_bw()+
  ylim(-2,25)+
  xlim(-1,16)+
  theme(text = element_text(size=9),
        plot.title = element_text(size=8, vjust=-5, hjust=0.01),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        legend.position="none")+
  ggrepel::geom_label_repel(aes(label = ifelse(RecEng>7,as.character(GEO_NAME),ifelse(RecRel>7,as.character(GEO_NAME),''))),
                   size=1.6,
                   force=2.5,
                   box.padding   = 0.4, 
                   point.padding = 0,
                   segment.color = 'grey50',
                   min.segment.length = 0.01,
                   max.overlaps=Inf,
                   label.size = NA,
                   fill = alpha(c("white"),0.1))+
  ggtitle("poverty index")



MA_rec_pers <- ggplot(dat_rec_e_MA %>% filter(personal_disruption_rank == "high"|
                                                personal_disruption_rank == "med high"|
                                                personal_disruption_rank =="med"), 
                      aes(x=RecEng, y=RecRel)) + 
  geom_point(shape=1, 
             size=1.4, stroke=0.9, alpha=0.6, 
             aes(color=personal_disruption_rank, ))+
  scale_color_manual(values = c("high" = "firebrick4",
                                "med high"="darkorange1",
                                "med"="goldenrod1",
                                "low"="skyblue2")) +
  theme_bw()+
  ylim(-2,25)+
  xlim(-1,16)+
  theme(text = element_text(size=9),
        plot.title = element_text(size=8, vjust=-7.5, hjust=0.01),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        legend.position="none")+
  ggrepel::geom_label_repel(aes(label = ifelse(RecEng>7,as.character(GEO_NAME),ifelse(RecRel>7,as.character(GEO_NAME),''))),
                   size=1.6,
                   force=8,
                   box.padding   = 0.4, 
                   point.padding = 0,
                   segment.color = 'grey50',
                   min.segment.length = 0.01,
                   max.overlaps=Inf,
                   label.size = NA,
                   fill = alpha(c("white"),0.1))+
  ggtitle("personal disruption index")



MA_rec_pop <- ggplot(dat_rec_e_MA %>% filter(pop_composition_rank == "high"| 
                                               pop_composition_rank == "med high"|
                                               pop_composition_rank =="med"), 
                     aes(x=RecEng, y=RecRel)) + 
  geom_point(shape=1, 
             size=1.4, stroke=0.9, alpha=0.6, 
             aes(color=pop_composition_rank, ))+
  scale_color_manual(values = c("high" = "firebrick4",
                                "med high"="darkorange1",
                                "med"="goldenrod1",
                                "low"="skyblue2")) +
  theme_bw()+
  ylim(-2,25)+
  xlim(-1,16)+
  theme(text = element_text(size=9),
        plot.title = element_text(size=8, vjust=-7.5, hjust=0.01),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        legend.position="right")+
  guides(color=guide_legend(title="Vulnerability\nIndicator\nScore Rank"))+
  ggrepel::geom_label_repel(aes(label = ifelse(RecEng>7,as.character(GEO_NAME),ifelse(RecRel>7,as.character(GEO_NAME),''))),
                   size=1.6,
                   force=2.5,
                   box.padding   = 0.4, 
                   point.padding = 0,
                   segment.color = 'grey50',
                   min.segment.length = 0.01,
                   max.overlaps=Inf,
                   label.size = NA,
                   fill = alpha(c("white"),0.1))+
  ggtitle("population composition index")





#put them together
#tiff("MA_rec_1-23-25.tiff", units="in", width=8, height=5, res=300)

MA_rec_plot + (MA_rec_poverty / MA_rec_pers / MA_rec_pop) +
  plot_layout(ncol=2, widths = c(5,2))

#dev.off()

```

Both commercial and recreational fishing are important activities in Montauk, NY, Cape May, Barnegat Light and Point Pleasant Beach, NJ; and Rodanthe, NC, meaning these communities may be impacted simultaneously by commercial and recreational regulatory changes. However, in all but Point Pleasant Beach NJ,  environmental justice may not be a major concern in these communities given the index scores (Figs \ref{fig:commercial-EJ} and \ref{fig:recreational-EJ})). Point Pleasant Beach, NJ scored medium for the personal disruption index, indicating that environmental justice may  be a moderate concern in Point Pleasant Beach.


### Indicators: Community Climate Vulnerability in the Mid-Atlantic

Coastal fishing communities worldwide have or are likely to experience social, economic, and cultural impacts from climate change, both negative (e.g., loss of infrastructure, fish stock decline) and positive (e.g., increased abundance of valuable species). Changes in marine fisheries as a consequence of climate change will require adaptation by coastal fishing communities and fisheries managers alike. The Community Climate Change Vulnerability Indicators were developed to help examine trends in climate change vulnerability in U.S. coastal fishing communities in the Northeast Region using indicators developed to understand fishing community level risk to climate change as based on species dependency.

```{r commvulex, fig.width = 7.5, fig.asp=.5, fig.cap="Proportion of Mid-Atlantic communities at each revenue climate vulnerability level over time."}

#CCCVI <- readxl::read_excel("CCCVI-NE/NE_CCCVI_00_22.xlsx", sheet = "NE_CCCVI_Final00_22")



# CCCVImid <- ecodata::community_climate_vulnerability |> 
#   dplyr::filter(State %in% c("MD", "NJ", "VA","NY", "NC", "DE")) 
# 
# CCCVImid[CCCVImid == 0] <- NA

vulquants <- data.frame(ymin = c(1,1.6,2.6,3.6),
                        ymax = c(1.5,2.5,3.5,4),
                        category = c("Low",
                                     "Moderate",
                                     "High",
                                     "Very High"))

risk = c("darkgreen", "yellow3", "orange", "red")
names(risk) <- c("Low","Moderate","High","Very High")


props <- ecodata::community_climate_vulnerability |>
  dplyr::filter(!is.na(Value), #remove communities with no ranking in a given year
                Var %in% "vulnSum") |> # hardcoded for total vulnerability
  #dplyr::select(Time, Port, State, Var) |>
  dplyr::group_by(Time, EPU) |>
  dplyr::mutate(commn = n(),
                vulCat = dplyr::case_when(Value <= 1.5 ~ "Low",
                                          Value > 3.5 ~ "Very High",
                                          Value > 1.5 & Value <= 2.5 ~ "Moderate",
                                          Value > 2.5 & Value <= 3.5 ~ "High")) |>
  dplyr::group_by(Time, EPU, vulCat) |>
  dplyr::mutate(catn = n(),
                propcat = catn/commn)

propsonly <- props |>
  dplyr::select(Time, vulCat, propcat) |>
  dplyr::distinct() |>
  dplyr::filter(EPU == "MAB")

ggplot2::ggplot() + #, group=PORT_NAME
  #ggplot2::geom_line(aes(group = PORT_NAME), alpha = 0.05) +
  # ggplot2::geom_rect(data = vulquants,
  #                    aes(ymin = ymin, ymax = ymax, fill = category, xmin = -Inf, xmax = Inf),
  #                    alpha = .3) +
  ggplot2::scale_fill_manual(values=risk) +
  ggplot2::geom_point(data = propsonly, ggplot2::aes(x=Time, y=propcat, color=vulCat)) +
  ggplot2::geom_line(data = propsonly, ggplot2::aes(x=Time, y=propcat, color=vulCat)) +
  ggplot2::scale_color_manual(values=risk, 
                              breaks=c('Low', 'Moderate', 'High', 'Very High'),
                              name="Climate Vulnerability") +
  ggplot2::facet_wrap(~EPU) +
  ecodata::geom_lm(n=22) + #this doesnt work yet
  #ggplot2::ylim(1,4) +
  ggplot2::theme_bw() + 
  #ggplot2::ggtitle("Mid-Atlantic") +
  #ggplot2::theme(legend.position = "blank")
  ggplot2::ylab("Proportion communities")

```


### Implications
A range of socioeconomic and climate vulnerability concerns are found throughout Mid-Atlantic fishing communities. However, index scores for socioeconomic concerns are higher overall in the top commercial communities relative to the top recreational communities. 

These indicators provide a snapshot of the presence of socioeconomi issues in the most highly engaged and most highly reliant commercial and recreational fishing communities in the Mid-Atlantic. These communities may be especially vulnerable to changes in fishing patterns due to regulations and/or climate change. Some changes occurred among the top fishing communities due to shifts in fishing activities, both commercial and recreational. Many of these communities, especially top commercial fishing communities, demonstrated medium to high environmental justice vulnerability, indicating that they may be at a disadvantage when responding to change.

There is evidence that a majority of Mid-Atlantic communities have high to very high revenue climate vulnerability. Coastal fishing communities are greatly affected by climate change, both because of their physical location and because of their frequent social, cultural, and economic dependence on fishing. These impacts are expected to become more pressing as climatic changes become more extensive. Changes in ocean temperature and acidification affecting marine life have the potential to directly impact fisheries and fishery dependent livelihoods. 


## Protected Species
Fishery management objectives for protected species generally focus on reducing threats and on habitat conservation/restoration. Protected species include marine mammals protected under the Marine Mammal Protection Act, endangered and threatened species protected under the Endangered Species Act, and migratory birds protected under the Migratory Bird Treaty Act. In the Northeast U.S., endangered/threatened species include Atlantic salmon, Atlantic and shortnose sturgeon, all sea turtle species, and five baleen whales. Protected species objectives include managing bycatch to remain below potential biological removal (PBR) thresholds, recovering endangered populations, and monitoring unusual mortality events (UMEs). Here we report on performance relative to these objectives with available indicator data, as well as indicating the potential for future interactions driven by observed and predicted ecosystem changes in the Northeast U.S. 

### Indicators: bycatch, population (adult and juvenile) numbers, mortalities
Average indices for both [harbor porpoise](https://noaa-edab.github.io/catalog/harborporpoise.html) (Fig. \ref{fig:harborporpoise}) and [gray seal](https://noaa-edab.github.io/catalog/grayseal.html) bycatch (Fig. \ref{fig:grayseal}) are below current PBR thresholds, meeting management objectives. 
```{r harborporpoise, fig.cap="Harbor porpoise average bycatch estimate for Mid-Atlantic and New England gillnet fisheries (blue) and the potential biological removal (red).", fig.width=5}
# code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/macrofauna_MAB.Rmd-harborporpoise.R")
#fig.width=6}

ecodata::plot_harborporpoise()
```

The annual estimate for gray seal bycatch has declined since 2019, in part driven by declining gillnet landings. In addition, estimates since 2019 have greater uncertainty stemming from low observer coverage since 2019. The rolling mean confidence interval remains just below the PBR threshold.
```{r grayseal, fig.cap="Gray Seal average bycatch estimate for gillnet fisheries (blue) and and the potential biological removal (red).", fig.width=5}
# code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/macrofauna_MAB.Rmd-grayseal.R")
#fig.width=6}

ecodata::plot_grayseal()
```

The [North Atlantic right whale population](https://noaa-edab.github.io/catalog/narw.html) was on a recovery trajectory until 2010, but has since declined (Fig. \ref{fig:narw-abundance}). The sharp decline observed from 2015-2020 appears to have slowed, although the right whale population continues to experience annual mortalities above recovery thresholds. Reduced survival rates of adult females lead to diverging abundance trends between sexes. It is estimated that there are fewer than 70 adult females remaining in the population.  
```{r narw-abundance, fig.cap = "Estimated North Atlanic right whale abundance on the Northeast Shelf."}
#code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/macrofauna_MAB.Rmd-narw-abundance.R")

ecodata::plot_narw(varName = "adult", n=10) + ggplot2::ggtitle("North Atlantic right whale abundance")
```
North Atlantic right whale [calf counts](https://noaa-edab.github.io/catalog/narw.html) have generally declined after 2009 to the point of having zero new calves observed in 2018 (Fig. \ref{fig:NARW-calf-abundance}). However, since 2019, we have seen more calf births each year with 15 births in 2022.
```{r NARW-calf-abundance, fig.cap = "Number of North Atlantic right whale calf births, 1990 - 2022."}
#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/macrofauna_MAB.Rmd-NARW-calf-abundance.R")

ecodata::plot_narw(varName = "calf", n=10) + ggplot2::ggtitle("North Atlantic right whale calf abundance")
```

This year, the Unusual Mortality Event (UME) for North Atlantic right whales continued. From 2017 through 16 February 2024, the total UME right whale mortalities includes 38 dead stranded whales, 17 in the US and 21 in Canada. When alive but seriously injured whales (34) and sublethal injuries or ill whales (51) are taken into account, 123 individual whales are included in the UME. Recent research suggests that many mortalities go unobserved and the true number of mortalities are about three times the count of the observed mortalities. The primary cause of death is “human interaction” from entanglements or vessel strikes.

A UME continued from previous years for humpback whales (2016-present) and Atlantic minke whales (2018-present); suspected causes include human interactions. A UME for Northeast pinnipeds that began in 2018 for infectious disease is pending closure as of February 2024.

### Implications
Bycatch management measures have been implemented to maintain bycatch below PBR thresholds. The downward trend in harbor porpoise bycatch could also be due to a decrease in harbor porpoise abundance in U.S. waters, reducing their overlap with fisheries, and a decrease in gillnet effort. The increasing trend in gray seal bycatch may be related to an increase in the gray seal population ([U.S. pup counts](https://noaa-edab.github.io/catalog/seal_pups.html)), supported by the dramatic rise over the last three decades in observed numbers of gray seal pups born at U.S. breeding sites plus an increase in adult seals at the breeding sites, some of which are supplemented by Canadian adults. 

Strong evidence exists to suggest that interactions between right whales and both the fixed gear fisheries in the U.S. and Canada and vessel strikes in the U.S. are contributing substantially to the decline of the species. Further, right whale distribution has changed since 2010. [New research](https://noaa-edab.github.io/catalog/narw.html) suggests that recent climate driven changes in ocean circulation have resulted in right whale distribution changes driven by increased warm water influx through the Northeast Channel, which has reduced the primary right whale prey (the copepod Calanus finmarchicus) in the central and eastern portions of the Gulf of Maine. Additional potential stressors include offshore wind development, which overlaps with important habitat areas used year-round by right whales, including mother and calf migration corridors and foraging habitat. This area is also the only known right whale winter foraging habitat. Additional information can be found in the  [offshore wind risks section](#other-ocean-uses-offshore-wind).

The UMEs are under investigation and are likely the result of multiple drivers. For all large whale UMEs, human interaction appears to have contributed to increased mortalities, although investigations are not complete.

A climate vulnerability assessment is published for Atlantic and Gulf of Mexico marine mammal populations.

# Risks to Meeting Fishery Management Objectives

## Climate and Ecosystem Change

Regulations and measures designed to meet fishery management objectives are often based on historical information about stocks, their distribution in space and time, and their overall productivity. Large scale climate related changes in the ecosystem can lead to changes in important habitats and ecological interactions, altering distributions and productivity. With large enough ecosystem changes, management measures may be less effective, and management objectives may not be met.

This year, we have restructured this section to focus on three categories of management decisions and the risk posed to them by climate and ecosystem change: spatial management, seasonal management, and quota setting or rebuilding depleted stocks. In each section, we describe potential risks to a management category, highlight indicators of observed changes that contribute to those risks, review possible biological and environmental drivers and the ways they may explain the observed indicators, and raise potential future implications if these trends persist or change. 

### Risks to Spatial Management

Shifting species distributions (changes in spatial extent or center of gravity) alter both species interactions and fishery interactions. In particular, shifting species distributions can affect expected management outcomes from spatial allocations and bycatch measures based on historical fish and protected species distributions. Species availability to surveys can also change as distributions shift within survey footprints.

#### Indicators: Fish and protected species distribution shifts
As noted in the [Landings Implications section above](#implications), the center of [distribution](https://noaa-edab.github.io/catalog/species_dist.html) for a suite of 48 commercially or ecologically important fish species along the entire Northeast Shelf continues to show movement towards the northeast and generally into deeper water (Fig. \ref{fig:species-dist}). 
  
```{r species-dist, fig.cap = "Aggregate species distribution metrics for species in the Northeast Large Marine Ecosystem: along shelf distance with increasing trend (orange), and depth with decreasing trend indicating deeper water (purple).", fig.width = 8, fig.asp=0.3}

#, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/macrofauna_MAB.Rmd-species-dist.R")

a <- ecodata::plot_species_dist(varName = "along", n=10) + ggplot2::coord_cartesian(xlim = c(1969,2021))
b <- ecodata::plot_species_dist(varName = "depth", n=10) + ggplot2::coord_cartesian(xlim = c(1969,2021)) 

a+b
```  

[Habitat model-based species richness](https://noaa-edab.github.io/catalog/habitat_diversity.html) suggests shifts of both cooler and warmer water species to the northeast. Similar patterns have been found for [marine mammals](https://noaa-edab.github.io/catalog/HMS_species_distribution.html), with multiple species shifting northeast between 2010 and 2017 in most seasons (Fig. \ref{fig:protectedspp-dist-shifts}).

```{r protectedspp-dist-shifts, fig.cap="Direction and magnitude of core habitat shifts, represented by the length of the line of the seasonal weighted centroid for species with more than 70 km difference between 2010 and 2017 (tip of arrow).", fig.width=7, fig.asp=0.8}
#, code=readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/macrofauna_MAB.Rmd-protectedspp-dist-shifts.R")

ecodata::plot_cetacean_dist() + ggplot2::ggtitle("Whale and Dolphin Distribution Shifts")
```


#### Drivers: 

Mobile populations shift distributions to maintain suitable temperature and prey fields, possibly expanding ranges if new suitable habitat exists.   Changes in managed species distribution is related, in part, to the [distribution of forage biomass](https://noaa-edab.github.io/catalog/forage_index.html). Since 1982, the fall center of gravity of forage fish (20 species combined) has moved to the north and east (Fig. \ref{fig:forageshifts}). Spring forage fish center of gravity has been more variable over time. 

```{r forageshifts, fig.cap="Eastward (left) and northward (right) shifts in the center of gravity for 20 forage fish species on the Northeast U.S. Shelf, with increasing trend (orange) for fall eastward and northward center of gravity.", fig.width = 8, fig.asp = .3}

ecodata::plot_forage_index(varName = "cog", n=10)
```

ADD COG FOR BENTHOS AND ZOOPLANKTON?

In contrast, macrobenthos center of gravity has shifted westward. Macrobenthos are small bottom-dwelling invertebrates including polychaete worms, small crustaceans, bivalves (non-commercial), gastropods, nemerteans, tunicates, cnidarians, brittle stars, sea cucumbers, and sand dollars. Large copepods have a similar pattern to macrobenthos (link to catalog).
```{r macrobenthosshifts, fig.alt="Eastward (left) and northward (right) shifts in the center of gravity for macrobenthos species on the Northeast U.S. Shelf", fig.width = 8,fig.asp = .3}
#, fig.asp = .3
ecodata::plot_benthos_index(plottype = "cog", varName = "Macrobenthos", n=10)+ ggplot2::coord_cartesian(xlim = c(1982, 2022))
```

While small Copepods center of gravity from ECOMON shows generally similar shifts to forage fish, [link to catalog dont need to show this].
```{r Smcopeshifts, fig.alt="Eastward (left) and northward (right) shifts in the center of gravity for small copepod species on the Northeast U.S. Shelf", fig.width = 8,fig.asp = .3}
#, fig.asp = .3
ecodata::plot_zooplankton_index(plottype = "cog", varName = "Smallcopeall", n=10)+ ggplot2::coord_cartesian(xlim = c(1982, 2022))
```

Ocean temperatures influence the distribution, seasonal timing, and productivity of managed species (see sections below). The Northeast US shelf, including the Mid-Atlantic, has experienced a continued warming trend for both the [long term](https://noaa-edab.github.io/catalog/long_term_sst.html) (Fig. \ref{fig:long-term-sst}) and [recent surface](https://noaa-edab.github.io/catalog/seasonal_oisst_anom.html) and [bottom](https://noaa-edab.github.io/catalog/bottom_temp_comp.html) in all seasons.

```{r long-term-sst, fig.cap="Northeast US annual sea surface temperature (SST, black), with increasing trend (orange)."}

ecodata::plot_long_term_sst(n=10)
```


Species suitable habitat can expand or contract when changes in temperature and major oceanographic conditions alter distinct water mass habitats.The variability of the Gulf Stream is a major driver of the predominant oceanographic conditions of the Northeast U.S. continental shelf. As the [Gulf Stream](https://noaa-edab.github.io/catalog/gsi.html) has become less stable and shifted northward in the last decade (Fig. \ref{fig:GSI}), warmer ocean temperatures have been observed on the northeast shelf and a higher proportion of [Warm Slope Water](https://noaa-edab.github.io/catalog/slopewater.html) has been present in the Northeast Channel.

```{r GSI, fig.cap = "Index representing changes in the location of the western Gulf Stream north wall (black). Positive values represent a more northerly Gulf Stream position, with increasing trend (orange)."}
#, code=readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/LTL_MAB.Rmd-gsi.R")

ecodata::plot_gsi(varName = "westgsi", n=10)
```


Changes in ocean temperature and circulation alter habitat features such as the seasonal [cold pool](https://noaa-edab.github.io/catalog/cold_pool.html), a band of relatively cold near-bottom water from spring to fall over the  northern MAB. The cold pool represents essential fish spawning and nursery habitat, and affects fish distribution and behavior. The cold pool has been getting warmer and smaller over time (Fig. \ref{fig:cold-pool-size}). The spatial extent (or area) of the seasonal cold pool is decreasing over time, yet the interannual variability of cold pool area has increased.

```{r cold-pool-size, fig.cap="Seasonal cold pool mean temperature (left) and spatial extent index (right), based on bias-corrected ROMS-NWA (open circles) and GLORYS (closed circles), with declining trends (purple).", fig.width = 7.5, fig.asp = 0.3, results='hide'}
#,  code=readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/LTL_MAB.Rmd-cold_pool.R")

a <- ecodata::plot_cold_pool(varName = "cold_pool", n=10)
b <- ecodata::plot_cold_pool(varName = "extent", n=10)

a+b
```

#### Future Considerations 

Distribution shifts caused by changes in thermal habitat are likely to continue as long as long-term temperature trends persist. Near-term oceanographic forecasts are currently in development and may inform how future warming impacts species distributions.

Distribution patterns associated with climate-driven changes in ocean circulation are also unlikely to be reversed to historical ranges in the short term. Increased oceanographic variability needs to be captured by regional ocean models and linked to species distribution processes to better understand potential future distributions. Species with high mobility or short lifespans react differently from immobile or long lived species.

Adapting management to changing stock distributions and dynamic ocean processes will require continued monitoring of populations in space and evaluating management measures against a range of possible future spatial distributions. Processes like the [East Coast Climate Scenario Planning](https://www.mafmc.org/climate-change-scenario-planning) can help coordinate management.


### Risks to Seasonal Management

The effectiveness of seasonal management actions (fishing seasons or area opening/closing) depends on a proper alignment with the seasonal life cycle events (phenology) of fish stocks (e.g. migration timing and spawning). Changes in the timing of these biological cycles can reduce the effectiveness of management measures if not accounted for. The timing of seasonal patterns can also change the interactions between fisheries and non-target species thus influencing the amount of bycatch and the availability of species to surveys. 

#### Indicators: Timing shifts

[Spawning timing](https://noaa-edab.github.io/catalog/spawn_timing.html) is shifting earlier for multiple stocks, including haddock and yellowtail flounder. Spawning of both haddock stocks occurred earlier in the year, as indicated by more resting (post-spawning) stage fish in the 2010s as compared to earlier in the time series (Fig. \ref{fig:spawntiming}). The northern (CC/GOM) stock shows earlier active spawning in recent years with a decline in pre-spawning resting females.The recent increase in resting females in the southern (SNE) stock also indicates a shift to earlier spawning (i.e. more post-spawn fish). Yellowtail flounder spawning is related to bottom temperature, week of year, and decade sampled for each of the three stocks.

```{r spawntiming, fig.cap="Percent resting stage (non-spawning) mature female fish (black) with significant increases (orange) and decreases (purple) from two haddock and three yellowtail flounder stocks: CC = Cape Cod Gulf of Maine, GOM = Gulf of Maine, GB = Georges Bank, SNE = Southern New England. ", fig.width = 7, fig.asp=.5}

ecodata::plot_spawn_timing(n=10) + ggplot2::ggtitle("Spring Resting Maturity Stage")

```

[Migration timing](https://noaa-edab.github.io/catalog/timing_shifts.html) of some tuna and large whale migrations has changed. For example, tuna were caught in recreational fisheries 50 days earlier in the year in 2019 compared to 2002. In Cape Cod Bay, peak spring habitat use by right and humpback whales has shifted 18-19 days later over time. 

Understanding whether seasonal patterns are changing for stocks requires regular observations throughout the year.  Despite the importance of understanding seasonal patterns,  we have few indicators that directly assess timing shifts of species. We plan on incorporating more indicators of timing shifts and phenology in future reports.


#### Drivers: 

The drivers of timing shifts in managed stocks are generally coupled to shifts in environmental or biological conditions, since these can result in changes in habitat quality or food availability within the year. Changes in the timing of fall phytoplankton blooms and seasonal shifts in zooplankton communities are thought to be critical indicators of changes in seasonal food availability to stocks. 

Along with the overall warming trends in the Mid Atlantic, ocean summer conditions have been lasting longer, as shown by the later [transition](https://noaa-edab.github.io/catalog/trans_dates.html) from warm stratified summer conditions to well mixed cool fall conditions (Fig. \ref{fig:transition}). Changes in the broad seasonal cycles of their environment can lead to changes in species biological processes (migrations, spawning, etc.) that are triggered by seasonal events.  

```{r transition, fig.cap="Ocean summer length in the MAB: the annual total number of days between the spring thermal transition date and the fall thermal transition date (black), with an increasing trend (orange)."}
#, code=readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/LTL_MAB.Rmd-sumlength.R")
# The transition dates are defined as the day of the year when surface temperatures changeover from cool to warm conditions in the spring and back to cool conditions in the fall.

ecodata::plot_trans_dates(varName = "length", n=10) + ggplot2::ggtitle("Time between spring and fall transition: MAB") + theme(
  strip.background = element_blank(),
  strip.text.x = element_blank()
)
```

The [cold pool](https://noaa-edab.github.io/catalog/cold_pool.html) is a seasonal feature within the MAB that creates seasonally suitable habitat for many species. In 8 of the past 10 years, cold pool persistence has been well below average, so this habitat was available for a shorter portion of the year  (Fig. \ref{fig:cold-pool-time}). A change in the timing of the cold pool may impact the recruitment of species that depend on it for juvenile habitat, such as yellowtail flounder.

```{r cold-pool-time, fig.cap="Cold pool persistence index based on bias-corrected ROMS-NWA (open circles) and GLORYS (closed circles).", fig.width = 3.5, fig.asp=.5, results='hide'}
#, fig.width = 5,  code=readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/LTL_MAB.Rmd-cold_pool.R")

ecodata::plot_cold_pool(varName = "persistence", n=10)

```

#### Future Considerations

For stocks reliant on environmental processes to dictate the timing of their behavior (e.g. phytoplankton bloom timing, thermal transition, or the duration of the cold pool), it is possible that some changes are episodic and have interannual variability, while other effects on timing can change on scales of years to decades. However, other species may rely on the general seasonal succession of their environment, which exhibits long-term trends unlikely to reverse in coming years. For those species, timing shifts in migration or spawning may continue. Management actions that rely on effective alignment of fisheries availability and biological processes should continue to evaluate whether prior assumptions on seasonal timings still hold, and new indicators should be developed to monitor timing shifts for stocks.

### Risks to Quota Setting/Rebuilding

The efficacy of short-term stock projections and rebuilding plans rely on an accurate understanding of processes affecting stock growth, reproduction, and natural mortality. These biological processes are often driven by underlying environmental change. When observed environmental change occurs, there is a risk that established stock-level biological reference points may no longer reflect the current population. 

#### Indicators: Fish productivity and condition shifts
Indicators of [fish productivity](https://noaa-edab.github.io/catalog/productivity_anomaly.html) are derived from observations (surveys) or models (stock assessments). Fish productivity has been declining in the Mid-Atlantic since the early 2000s, as described by the small-fish-per-large-fish anomaly indicator (derived from NEFSC bottom trawl survey) (Fig. \ref{fig:productivity-anomaly}). This decline in fish productivity is also shown by a similar analysis based on stock assessment model outputs (recruitment per spawning stock biomass anomaly).

```{r productivity-anomaly, fig.cap = "Fish productivity measures. Left: Small fish per large fish survey biomass anomaly in the Mid-Atlantic Bight. Right: assessment recruitment per spawning stock biomass anomaly for stocks mainly in the Mid-Atlantic. The summed anomaly across species is shown by the black line, drawn across all years with the same number of stocks analyzed.", fig.width=8, fig.asp=0.6}
#out.width='49%', fig.show='hold',
a <- ecodata::plot_productivity_anomaly(report = "MidAtlantic") + 
  ggplot2::guides(fill=guide_legend(ncol=2)) + 
  ggplot2::theme(legend.position = "bottom",
                      legend.title = ggplot2::element_blank(),
                 plot.title =element_text(size = 11))

b <- ecodata::plot_productivity_anomaly(report = "MidAtlantic", varName = "assessment")+ 
  ggplot2::guides(fill=guide_legend(ncol=2)) +
  ggplot2::theme(legend.position = "bottom",
                      legend.title = ggplot2::element_blank(),
                 plot.title =element_text(size = 11))

a + b

```

The health of individual fish (i.e. fish condition) can contribute to population productivity through improved growth, reproduction and survival. [Fish condition](https://noaa-edab.github.io/catalog/condition.html) in the MAB was generally good prior to 2000, poor from 2001-2010 (concurrent with declines in productivity, Fig. \ref{fig:productivity-anomaly}), and a mix of good and poor since 2011. In 2023, condition was mixed, with general improvement since a relatively low condition year in 2021 (Fig. \ref{fig:mab-cf}). Preliminary analyses show that changes in temperature, zooplankton, fishing pressure, and population size influence the condition of different fish species.

```{r mab-cf, fig.cap = "Condition factor for fish species in the MAB based on fall NEFSC bottom trawl survey data. MAB data are missing for 2017 due to survey delays, and no survey was conducted in 2020.", fig.width=9}

ecodata::plot_condition() +
  theme(#legend.position = 'bottom',
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 11),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 10),
        plot.title = element_text(size = 12))
```

#### Drivers: 

Fish productivity and condition are affected by increasing metabolic demands from increasing temperature, combined with changes in the availability and quality of prey.  Long-term environmental trends and episodic extreme temperatures, ocean acidification, and low oxygen events represent multiple stressors that can affect growth rates, reproductive success, recruitment, and cause mortality.

##### Biological Drivers: Forage quality and abundance
The amount of forage fish available in the ecosystem combined with the energy content of the forage species determines the amount of energy potentially available to predators in the ecosystem. Changes in the forage base can drive managed and protected species production. 

The [energy content](https://noaa-edab.github.io/catalog/energy_density.html) of juvenile and adult forage fish as prey is related to forage fish growth and reproductive cycles, as well as environmental conditions.  The energy content of Atlantic herring from the NEFSC trawl surveys has increased recently (Fig. \ref{fig:energy-density}) but is still well below that observed in the 1980s and 1990s. Silver hake, longfin squid (Loligo in figure) and shortfin squid (Illex in figure) remain lower than previous estimates.
  
```{r energy-density, fig.cap="Forage fish energy density mean and standard deviation by season and year, compared with 1980s (solid line) and 1990s (dashed line) values.",  fig.width = 7.5, fig.asp = 0.5}
#code=readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/macrofauna_MAB.Rmd-energy-density.R"),

ecodata::plot_energy_density()
```

Changes in the overall abundance of forage fish can influence managed species productivity as it relates to changes in food availability. A spatially-explicit [forage index](https://noaa-edab.github.io/catalog/forage_index.html) for the Mid-Atlantic shows a long term declining trend in fall, with higher forage biomass in fall than spring (Fig. \ref{fig:foragebio}). Forage biomass was highest during fall in the early-1980s.   
```{r foragebio, fig.cap = "Forage fish index in the MAB for spring (blue) and fall (red) surveys, with a decline (purple) in fall. Index values are relative to the maximum observation within a region across surveys.", fig.asp=.6}
#, code=readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/macrofauna_MAB.Rmd-forage-index.R")

ecodata::plot_forage_index(n=10)
```

ADD BENTHOS FORAGE HERE?

Benthic invertebrates are extremely important forage for some managed species (Black Sea Bass,...) Macrobenthos indices show long term declines in spring. In contrast, Megabenthos indices show long term increases in spring in the MAB.


```{r benthos, fig.alt="Changes in benthos abundance in the MAB for megabenthos (left) and macrobenthos (right).", fig.width=7.5, fig.asp=.45}

c <- ecodata::plot_benthos_index(report = "MidAtlantic", varName = "Megabenthos", n=10) +  ggplot2::theme(legend.position = "bottom")
d <- ecodata::plot_benthos_index(report = "MidAtlantic", varName = "Macrobenthos", n=10) +  ggplot2::theme(legend.position = "none")


c+d

```

##### Biological Drivers: Lower trophic levels

[Phytoplankton](https://noaa-edab.github.io/catalog/chl_pp.html) are the foundation of the food web and are the primary food source for zooplankton and filter feeders such as shellfish. Numerous environmental and oceanographic factors affect the abundance, [size composition](https://noaa-edab.github.io/catalog/phyto_size.html), spatial distribution, and productivity of phytoplankton. While changes in fish productivity (including forage) could result from changing primary productivity, total primary production in the Mid Atlantic has no clear trend (Fig. \ref{fig:totpp}).

[Zooplankton](https://noaa-edab.github.io/catalog/zoo_abundance_anom.html) communities in the Mid-Atlantic have increasing trends for krill (euphausiids; Fig. \ref{fig:zoopanom}) with recent decreasing trends for smaller bodied copepods. Small copepods are less energy-rich than Eupausiids (krill) or the larger-bodied copepod *Calanus finmarchicus*. A changing mix of zooplankton prey can impact forage fish energy content and abundance, as well as the prey field of filter feeding whales. 

```{r zoopanom, fig.cap="Changes in zooplankton abundance in the MAB for large (left) and small (canter) copepods, and Euphausiids (right), with significant decreases (purple) in small copepods and increases (orange) in Euphausiids.", fig.width=7.5, fig.asp=.5}
# a <- ecodata::plot_zoo_abundance_anom(report = "MidAtlantic", varName = "copepod") + 
#   ggplot2::facet_wrap(~EPU~Var, labeller = labeller(EPU = function(x) {rep("", length(x))}))
# b <- ecodata::plot_zoo_abundance_anom(report = "MidAtlantic", varName = "euphausid") + 
#   ggplot2::facet_wrap(~EPU~Var, labeller = labeller(EPU = function(x) {rep("", length(x))}))
# 
# a/b

a <- ecodata::plot_zooplankton_index(report = "MidAtlantic", varName = "Lgcopeall", n=10) +  ggplot2::theme(legend.position = "none")
b <- ecodata::plot_zooplankton_index(report = "MidAtlantic", varName = "Smallcopeall", n=10) +  ggplot2::theme(legend.position = "none")
c <- ecodata::plot_zooplankton_index(report = "MidAtlantic", varName = "Euph", n=10) +  ggplot2::theme(legend.position = "none")

a+b + c

```

##### Environmental Drivers

Fish production can also be directly related to the prevailing environmental conditions by altering metabolic processes (growth) and reproduction. Many species possess thermal tolerances and can experience stressful or lethal conditions if temperatures exceed certain levels. Extreme temperatures at both the [surface](https://noaa-edab.github.io/catalog/seasonal_oisst_anom.html) and [bottom](https://noaa-edab.github.io/catalog/bottom_temp_comp.html) can exceed [thermal tolerance](https://noaa-edab.github.io/catalog/thermal_habitat_persistence.html) limits for some fish. 

REPLACE WITH 2024 TEXT? 2024 didn't have a lot of hot water (Fig. \ref{fig:therm-hab-persist-2024}). OLD TEXT For example, 2015 had the [warmest summer and fall bottom temperatures](https://noaa-edab.github.io/catalog/bottom_temp_comp.html) in the Mid-Atlantic. A large proportion of the region had bottom temperatures above the 15℃ thermal tolerance for most groundfish, with some days exceeding the 24℃ potential mortality limit. Many Mid-Atlantic species have different thermal tolerance limits from groundfish, and we will work to include those next year.


```{r therm-hab-persist-2024, fig.cap = "The number of days in 2024 where bottom temperature exceeds 15℃ (left) and 24℃ (right) based on the GLORYS 1/12 degree grid.", fig.width=10, fig.asp=1}
# year = 2015
# p <- ecodata::plot_thermal_habitat_persistence(year=year) + ggplot2::ggtitle(paste(year))
# ggplot2::ggsave(here::here('images','thermal_habitat_2015.png'),p,width = 8,height = 8, units = 'in',dpi = 300)
# 
# magick::image_read(here::here('images','thermal_habitat_2015.png')) |>
#   magick::image_trim() |>
#   magick::image_write(here::here('images', 'thermal_habitat_2015_cropped.png'))

#knitr::include_graphics(here::here('images','thermal_habitat_2015_cropped.png'))

ecodata::plot_thermal_habitat_gridded()

```


REPLACE WITH NEW FIGS NEED WRITEUP POSTED
LOWEST ARAGONITE SATURATION OBSERVED SINCE MEASUREMENTS STARTED IN 2007
[Ocean acidification](https://noaa-edab.github.io/catalog/ocean_acidification.html) (OA) risks vary among species and include reduced survival, growth, reproduction, and productivity, where high OA risk indicates potential negative effects to species.  High OA risk conditions were observed for Atlantic sea scallop and longfin squid in Long Island Sound and the nearshore and mid shelf regions of the New Jersey shelf during summer of 2016, 2018, 2019, and 2023 (Fig. \ref{fig:mab-oa}).

```{r mab-oa, out.width = '100%', fig.cap = "Locations where bottom aragonite saturation state ($\\Omega_{Arag}$; summer only: June-August) were at or below the laboratory-derived sensitivity level for Atlantic sea scallop (left panel) and longfin squid (right panel) for the time periods 2007-2022 (dark cyan) and 2023 only (magenta). Gray circles indicate locations where bottom $\\Omega_{Arag}$ values were above the species specific sensitivity values.."}


#knitr::include_graphics("https://github.com/NOAA-EDAB/ecodata/raw/master/docs/images/Saba_Fig_SOE_MAFMC-GraceSaba.jpg")

knitr::include_graphics("https://github.com/NOAA-EDAB/ecodata/blob/dev/data-raw/workshop/images/SeasonalBottomAragoniteLongfin-GraceSaba_2025.png?raw=true")

#magick::image_read("https://github.com/NOAA-EDAB/ecodata/raw/master/docs/images/Saba_Fig_SOE_MAFMC - Grace Saba.jpg")

```

Biological and oceanographic processes can affect the amount of oxygen present in the water column. During low oxygen (hypoxic) events, species growth is negatively affected, and very low oxygen can result in mortality. The duration and extent of hypoxic events is being monitored, but long-term shelf-wide observations are not yet available. However, [hypoxic events](https://noaa-edab.github.io/catalog/ocean_acidification.html) were detected off the coast of New Jersey in 2023 and were potentially responsible for fish, lobster, and crab [mortalities](https://sebsnjaesnews.rutgers.edu/2023/12/rutgers-scientists-observe-unusual-ocean-conditions-possibly-linked-to-mortality-in-marine-life-off-new-jersey/).

#### Drivers: Predation
The abundance and distribution of predators can affect both the productivity and mortality rates on managed stocks. Predators can consume managed species or compete for the same resources, resulting in increased natural mortality or decreased productivity. The northeast shift in [whales and dolphins](https://noaa-edab.github.io/catalog/HMS_species_distribution.html)  (Fig. \ref{fig:protectedspp-dist-shifts}) indicates a change in the overlap between predators and prey. Since we also observe distribution shifts in managed species as well as forage species, the effect of changing predator distributions alone is difficult to quantify. 

Indicators for shark populations, combined with information on gray seals (see [Protected Species Implications section, above](#protected-species)), suggests predator populations range from stable ([sharks](https://noaa-edab.github.io/catalog/hms_cpue.html)) to increasing ([gray seals](https://noaa-edab.github.io/catalog/seal_pups.html))  in the MAB. [Stock status](https://noaa-edab.github.io/catalog/hms_stock_status.html) is mixed for Atlantic Highly Migratory Species (HMS) stocks (including sharks, swordfish, billfish, and tunas) occurring throughout the Northeast U.S. shelf. While there are several HMS species considered to be overfished or that have unknown stock status, the population status for some managed Atlantic sharks and tunas is at or above the biomass target, suggesting the potential for robust predator populations among these managed species. Stable predator populations suggest stable predation pressure on managed species, but increasing predator populations may reflect increasing predation pressure.

### Future Considerations
The processes that control fish productivity and mortality are dynamic, complex, and are the result of the interactions between multiple system drivers. There is a real risk that short-term predictions in assessments and rebuilding plans that assume unchanging underlying conditions will not be as effective, given the observed change documented in the prior sections in both ecological and environmental processes. Assumptions for species’ growth, reproduction, and natural mortality should continue to be evaluated for individual species. With observations of system-wide productivity shifts of multiple managed stocks, more research is needed to determine whether regime shifts or ecosystem reorganization are occurring, and how this should be incorporated into management.


## Other Ocean Uses: Offshore Wind

### Indicators: development timeline, revenue in lease areas, coastal community vulnerability
As of January 2024, 30 offshore [wind development](https://noaa-edab.github.io/catalog/wind_dev_speed.html) projects are proposed for construction over the next decade in the Northeast (timelines and project data for 2024 are based on the Ocean Wind 1 Offshore Wind Farm Final Environmental Impact Statement. Volume II: Appendix F). Offshore wind areas are anticipated to cover more than 2.3 million acres by 2030 in the Greater Atlantic region (Fig. \ref{fig:wind-proposed-dev}). All states will be able to reach their 2030 offshore wind goals with existing lease areas.  
```{r wind-proposed-dev, fig.cap='Total area proposed for wind develpment on the northeast shelf through 2030.', fig.width=5}
#, code=readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-wind-proposed-dev.R")

ecodata::plot_wind_dev_speed()
```

Just over 3,300 foundations and more than 12,000 miles of inter-array and offshore export cables are proposed to date (Fig. \ref{fig:wind-dev-cumul}). Based on current timelines, the areas affected would be spread out such that it is unlikely that any one particular area would experience full development at one time. Construction of two projects in Southern New England (South Fork Wind and Vineyard Wind 1) during 2023 affected fisheries managed by the Mid-Atlantic Fishery Management Council, while construction activities began for Revolution Wind in early 2024. It is likely that construction will begin on other projects in Southern New England and possibly the New York Bight during 2024 that will further affect regional fisheries. 

```{r wind-dev-cumul, fig.cap = "All Northeast Project areas by year construction ends (each project has 2 year construction period).", out.width='90%'}
#knitr::include_url("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/docs/images/All_2021128_needsgraph-01.jpg")
#knitr::include_graphics("https://github.com/NOAA-EDAB/ecodata/raw/master/docs/images/offshore_wind_timeline.png")

#magick::image_read("https://github.com/NOAA-EDAB/ecodata/raw/master/workshop/images/Cumulative_Timeline_FullRegion_SoE2024_v2_2024.png")

knitr::include_graphics("https://github.com/NOAA-EDAB/ecodata/raw/dev/data-raw/workshop/images/SoE2025_TimelineandFootprint_Sept2024-AngelaSilva_2025.pdf")

```

Based on federal vessel logbook data, [commercial fishery revenue](https://noaa-edab.github.io/catalog/wind_revenue.html) from trips in the current offshore wind lease areas, including the newly designated lease areas in the Central Atlantic, have varied annually from 2008-2022, with less than \$1 million in maximum annual revenue overlapping with these areas for most fisheries with the exception of the surfclam, monkfish, and longfin squid fisheries. Some fisheries see periodic spikes in revenue overlap with wind energy lease areas, including the surfclam (\$6.6 million), longfin squid (\$4.7 million), monkfish (\$4.3 million), and summer flounder (\$1.3 million) fisheries (Fig. \ref{fig:wea-spp-rev}). 

```{r wea-spp-rev, fig.cap="Fishery revenue in wind energy lease areas in the Mid-Atlantic.",  fig.width=5, fig.asp=.4}
#code=readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-wea-spp-rev.R"),

ecodata::plot_wind_revenue(report = "MidAtlantic", varName = "value", plottype = "nofacets")

```

Of MAFMC managed fisheries, the monkfish fishery would be the fishery most affected by offshore wind development, with a maximum of 20\% of annual regional fishery revenue occurring within existing and proposed wind lease areas and the Gulf of Maine Draft Wind Energy Area during 2008-2022 (see Table \ref{tab:wea-landings-rev}). Future fishery resource overlap with wind leases, especially surfclams and ocean quahogs, may change due to species distribution shifts attributable to climate change and recruitment and larval dispersion pattern changes caused by hydrodynamic flow disruptions from turbine foundations, which could also affect fishery landings/revenue.


```{r wea-landings-rev, ft.arraystretch = 1}
#fig.cap="Percent Landings and Revenue from wind energy areas. Data from GARFO VTR.",
#, out.width="90%", code=readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-wea-landings-rev.R")

# dt<- ecodata::wea_landings_rev[1:10, 1:5] %>% 
#   dplyr::select("GARFO and ASMFC Managed Species", "perc_landings_max"  ,"perc_landings_min",
#                 "perc_revenue_max","perc_revenue_min" ) %>% 
#   dplyr::rename("Maximum Percent Total Annual Regional Species Landings"="perc_landings_max",
#     "Maximum Percent Total Annual Regional Species Revenue"="perc_revenue_max",
#     "Minimum Percent Total Annual Regional Species Landings"="perc_landings_min",
#     "Minimum Percent Total Annual Regional Species Revenue"="perc_revenue_min",) 
# kable(dt,"latex", booktabs = TRUE,
#       caption = "Top ten species Landings and Revenue from Wind Energy Areas.") %>% 
#   kable_classic(full_width = F, html_font = "Cambria") %>%
#   column_spec(2:5, width = "10em") %>%
#   kable_styling(latex_options = "scale_down")

# dt<- ecodata::wea_landings_rev[1:20,] %>% 
#   dplyr::select("NEFMC, MAFMC, and ASMFC Managed Species", "perc_landings_max",
#                 "perc_revenue_max" ) %>%
#   dplyr::filter(`NEFMC, MAFMC, and ASMFC Managed Species` != "Ocean pout" ) %>%
#   dplyr::rename("Maximum Percent Total Annual Regional Species Landings"="perc_landings_max",
#     "Maximum Percent Total Annual Regional Species Revenue"="perc_revenue_max")
# 
# flextable::flextable(dt) %>%
#   flextable::set_caption("Top Species Landings and Revenue from Wind Energy Areas. * Landings and revenue for these species are likely underestimated due to limited coverage of these fisheries in historic reporting requirements for vessels issued federal permits by the NMFS Greater Atlantic Regional Fisheries Office. However, such limitations also suggest an inaccurately higher proportion of such landings and revenues in existing lease areas. ** Clearnose skates were reported separately from skates, which is presumed to include all skates managed under the Northeast skate complex. *** Based on comparison with other data sources, the high values for Illex squid are likely overestimates affected by the methods used to model logbook data to estimate spatial overlap of fishign operations with wind energy areas.") %>%
#   flextable::width(width = c(2,2,2)) #%>%
  # flextable::footnote(i=, j = 1,
  #              value = as_paragraph(
  #                c("* Landings and revenue for these species are likely underestimated due to limited coverage of these fisheries in historic reporting requirements for vessels issued federal permits by the NMFS Greater Atlantic Regional Fisheries Office. However, such limitations also suggest an inaccurately higher proportion of such landings and revenues in existing lease areas.",
  #                  "** Clearnose skates were reported separately from skates, which is presumed to include all skates managed under the Northeast skate complex.",
  #                  "*** Based on comparison with other data sources, the high values for *Illex* squid are likely overestimates affected by the methods used to model logbook data to estimate spatial overlap of fishign operations with wind energy areas.")))

# # table footnotes
#  * Landings and revenue for these species are likely underestimated due to limited coverage of these fisheries in historic reporting requirements for vessels issued federal permits by the NMFS Greater Atlantic Regional Fisheries Office. However, such limitations also suggest an inaccurately higher proportion of such landings and revenues in existing lease areas. 
# 
# ** Clearnose skates were reported separately from skates, which is presumed to include all skates managed under the Northeast skate complex. 
# 
# *** Based on comparison with other data sources, the high values for Illex squid are likely overestimates affected by the methods used to model logbook data to estimate spatial overlap of fishign operations with wind energy areas.

#ecodata::plot_wea_landings_rev(n=Inf)

dt<- ecodata::wea_landings_rev |> 
  dplyr::filter(Council %in% c("MAFMC", "MAFMC/NEFMC")) |>
  dplyr::select("NEFMC, MAFMC, and ASMFC Managed Species", "perc_landings_max",
                "perc_revenue_max" ) %>%
  dplyr::rename("Maximum Percent Total Annual Regional Species Landings"="perc_landings_max",
    "Maximum Percent Total Annual Regional Species Revenue"="perc_revenue_max")

flextable::flextable(dt) %>%
  flextable::set_caption("Mid-Atlantic managed species Landings and Revenue from Wind Energy Areas.") %>%
  flextable::width(width = c(2,2,2))

```

Proposed wind development areas interact with the region’s federal scientific surveys. Scientific surveys are impacted by offshore wind in four ways:  
1.  Exclusion of NOAA Fisheries’ sampling platforms from the wind development area due to operational and safety limitations.  
2.  Impacts on the random-stratified statistical design that is the basis for scientific assessments, advice, and analyses.  
3.  Alteration of benthic and pelagic habitats, and airspace in and around the wind energy development, requiring new designs and methods to sample new habitats.  
4.  Reduced sampling productivity through navigation impacts of wind energy infrastructure on aerial and vessel survey operations.  

Increased vessel transit between stations may decrease data collections that are already limited by annual days-at-sea day allocations. The total survey area overlap ranges from 1-70% for all Greater Atlantic federal surveys. The Gulf of Maine Cooperative Research Bottom Longline Survey (41%) and the Shrimp Survey (70%) have the largest percent overlap with the draft Gulf of Maine Wind Energy Areas. The remaining surveys range from 1-16% overlap.  Individual survey strata have significant interaction with wind areas, including the sea scallop survey (up to 96% of individual strata) and the bottom trawl survey (up to 60% strata overlap). Additionally, up to 50% of the southern New England North Atlantic right whale survey’s area overlaps with proposed project areas and a region-wide survey mitigation program is underway

Social vulnerabilities of communities are priority concerns with offshore wind development and fisheries impacts in the Northeast, and the impacts of offshore wind development are expected to differentially [impact specific coastal communities](https://noaa-edab.github.io/catalog/wind_port.html) (Fig. \ref{fig:wea-port-rev}). Additionally, impacts of offshore wind development may unevenly affect individual operators, with some permit holders deriving a much higher proportion of revenue from wind areas than the port-based mean.   

```{r wea-port-rev, fig.cap="Percent of Mid-Atlantic port revenue from Wind Energy Areas (WEA) in descending order from most to least port revenue from WEA. EJ = Environmental Justice.",  fig.width=6.5, fig.asp=.85, fig.align='left'}
#code=readLines("https://raw.githubusercontent.com/NOAA-EDAB/ecodata/master/chunk-scripts/human_dimensions_MAB.Rmd-wea-port-rev.R"),

# df.symbol <- ecodata::wind_port |>
#   dplyr::filter(EPU == "MAB",
#                 !Var %in% c("WEA_MAX", "TOT_MAX", 
#                                                      "perc_MIN", "perc_MAX")) |>
#   tidyr::pivot_wider( names_from = Var, values_from = Value) |>
#   dplyr::mutate(City = paste0(City,State)) |>
#   dplyr::select(City, EJ, Gent) |>
#   pivot_longer(cols = c(EJ, Gent), names_to = "Variable") |>
#   filter(!value == "NA") |>
#   dplyr::mutate(symbol = recode(Variable, EJ = -7, Gent = -3), 
#                 Variable = recode(Variable,"EJ"= "Mid-High to High EJ Concerns" , 
#                                   "Gent" ="Mid-High to High Gentrificaiton Concerns"))
# 
#   

ecodata::plot_wind_port(report = "MidAtlantic") #+
  # ggplot2::geom_point(data = df.symbol, aes(x = symbol,y = City, shape = Variable)) +
  # ggplot2::scale_shape_manual(values = c(17, 16))

```

For example, Atlantic City, NJ had the highest potential revenue loss (minimum of 10\% and maximum of 23\%) from potential wind development areas based on 2008-2022 total port fisheries revenue. BOEM reports that cumulative offshore wind development (if all proposed projects are developed) could have moderate impacts on low-income members of communities with environmental justice concerns who work in the commercial fishing and for-hire fishing industry due to disruptions to fish populations, restrictions on navigation and increased vessel traffic, as well as existing vulnerabilities of low-income workers to economic impacts.

THIS NEEDS AN UPDATE WE DON'T KNOW WHICH COMMUNITIES TO INCLUDE HERE

Some ports in New England land Mid-Atlantic managed species from wind areas as well. For the maximum percent value reported in each New England port, the majority (at least 50% based on both value and pounds) of those landings were Mid-Atlantic managed species within wind areas for Barnstable, MA, and Point Judith, RI. 

```{r wind-rev-NE-MAFMC, fig.cap="Percent of New England port revenue with majority MAFMC landings from Wind Energy Areas (WEA) in descending order from most to least port revenue from WEA. EJ = Environmental Justice.", fig.width=6.5, fig.asp=.3, fig.align='left'}
df.symbol <- ecodata::wind_port %>% filter(EPU == "NE", 
                                           !Var %in% c("MaxVal", "TOT_MAX", 
                                                     "perc_MIN", "perc_MAX")) %>% 
   tidyr::pivot_wider( names_from = Var, values_from = Value) %>%
  dplyr::mutate(City = paste0(City,",", State)) %>% 
  dplyr::filter(City %in% c("Point Judith, RI",
                          #  "North Kingstown RI",
                          #  "Boston MA",
                          #  "Woods Hole MA",
                          #  "Hyannis MA",
                            "Barnstable, MA")) %>%
  dplyr::select(City, EJ, Gentrification) %>% 
  pivot_longer(cols = c(EJ, Gentrification), names_to = "Variable") %>% 
  filter(!value == "NA") %>% 
  dplyr::mutate(symbol = recode(Variable, EJ = -7, Gentrification = -3), 
                Variable = recode(Variable,"EJ"= "Mid-High to High EJ Concerns" , 
                                  "Gentrification" ="Mid-High to High Gentrificaiton Concerns"))

# Percentage plot

df.all.perc<- ecodata::wind_port %>% filter(EPU == "NE") %>%
  tidyr::pivot_wider( names_from = Var, values_from = Value) %>%
  dplyr::mutate(ordering = MaxVal,
                   City = paste0(City, ",", State),
                   perc_dif =  c(perc_MAX - perc_MIN),
                   TOT_MAX = c(100 - perc_dif - perc_MIN)) %>% 
  tidyr::pivot_longer(cols = c(perc_MIN,  perc_dif, TOT_MAX), names_to="Var", values_to = "Value") %>% 
  dplyr::arrange(desc(ordering)) %>%
  dplyr::filter(City %in% c("Point Judith, RI",
                            #"North Kingstown RI",
                            #"Boston MA",
                            #"Woods Hole MA",
                            #"Hyannis MA",
                            "Barnstable, MA")) %>%
  dplyr::mutate(City = factor(City, levels = unique(City))) %>% 
  dplyr::filter(!Var %in% c("MaxVal", "EJ", "Gentrification")) %>% 
  dplyr::mutate(Var = recode(Var,"perc_MIN"= "WEA Revenue" , 
                             "perc_dif" ="WEA Revenue Range", 
                             "TOT_MAX" = "Non-WEA Revenue"), 
                Var = factor(Var, levels = c("Non-WEA Revenue", 
                                             "WEA Revenue Range", 
                                             "WEA Revenue")))
p2<-ggplot2::ggplot()+
  ggplot2::geom_bar(data = df.all.perc, aes(fill = Var, y = reorder(City, ordering), x = Value), stat="identity" )+
  ggplot2::scale_fill_brewer()+
  ggplot2::theme(legend.position = "bottom", 
                 legend.title = element_blank(), 
                 legend.box="vertical", legend.margin=margin())+
  ggplot2::geom_point(data = df.symbol, aes(x = symbol,y = City, shape = Variable)) +
  ggplot2::scale_shape_manual(values = c(17, 16)) +
  ggplot2::ggtitle("Port Revenue from WEA, majority MAFMC species")+
  ggplot2::xlab(expression("Port Revenue (%)"))+
  ggplot2::ylab(element_blank())+
  ecodata::theme_ts()

p2
```


Top fishing communities with [environmental justice concerns](https://noaa-edab.github.io/catalog/engagement.html) (i.e., Atlantic City, NJ, Newport News, VA, Hobucken and Beaufort, NC) should be considered in decision making to reduce the social and economic impacts and aid in the resilience and adaptive capacity of underserved communities. These are communities where we need to provide further resources to reach underserved and underrepresented groups and create opportunities for and directly involve these groups in the decision-making process. 

### Implications
Current plans for rapid buildout of offshore wind in a patchwork of areas spreads the impacts differentially throughout the region (Fig. \ref{fig:wind-dev-cumul}).  Up to 17\% of maximum annual fisheries revenue for major Mid-Atlantic commercial species in lease areas and draft call areas could be forgone or reduced and associated effort displaced if all sites are developed. Displaced fishing effort can alter historic fishing area, timing, and method patterns, which can in turn change habitat, species (managed and protected), and fleet interactions. Several factors, including fishery regulations, fishery availability, and user conflicts affect where, when, and how fishing effort may be displaced, along with impacts to and responses of affected fish species.

Planned development [overlaps NARW](https://noaa-edab.github.io/catalog/persistent_hotspots.html) mother and calf migration corridors and a significant foraging habitat that is used throughout the year (Fig. \ref{fig:whales-wind}). Turbine presence and extraction of energy from the system could alter local oceanography and may affect right whale prey availability. For example, persistent foraging hotspots of right whales and seabirds overlap on Nantucket Shoals, where unique hydrography aggregates enhanced prey densities. Wind leases (OCS-A 0521 and OCS-A 0522) currently intersect these hotspots on the southwestern corner of Nantucket Shoals and a prominent tidal front associated with invertebrate prey swarms important to seabirds and possibly right whales. Proposed wind development areas also bring increased vessel strike risk from construction and operation vessels. In addition, there are a number of potential impacts to whales from pile driving and operational noise such as displacement, increased levels of communication masking, and elevated stress hormones.

```{r whales-wind, out.width="55%", fig.cap="Northern Right Whale persistent hotspots (red shading) and Wind Energy Areas (black outlines)."}
#knitr::include_graphics("https://github.com/NOAA-EDAB/ecodata/raw/master/docs/images/right_whales.png")

#magick::image_read("https://github.com/NOAA-EDAB/ecodata/raw/master/docs/images/right_whales.jpg")

# NEW FILE 
magick::image_read("https://github.com/NOAA-EDAB/ecodata/raw/dev/data-raw/workshop/images/NARW_hotspots_final_2024.jpg")


```

Scientific data collection surveys for ocean and ecosystem conditions, fish, and protected species will be altered, potentially increasing uncertainty for stock assessments and associated management decision making.

The increase of offshore wind development can have both positive (e.g., employment opportunities) and negative (e.g., space-use conflicts) effects. Continued increase in coastal development and gentrification pressure has resulted in loss of fishing infrastructure space within ports. Understanding these existing pressures can allow for avoiding and mitigating negative impacts to our shore support industry and communities dependent on fishing. Some of the communities with the highest fisheries revenue overlap with offshore wind development areas that are also vulnerable to gentrification pressure are Point Pleasant and Atlantic City, NJ, Ocean City, MD, and Beaufort, NC.


### 2024 Highlights

This new section is common to the Mid-Atlantic and New England reports. Multiple [anomalous conditions](https://noaa-edab.github.io/catalog/observation_synthesis.html) and extreme events were observed in 2023 that could have brief local effects and/or widespread long-term ecosystem, fishery and management implications. This section intends to provide a record of these observations, the implications they may have for other ecosystem processes, and a reflection on how they fit into our understanding of the ecosystem. Many of these observations are being actively studied but should be noted and considered in future analyses and management decisions.

Globally, 2023 was the warmest year on record with record high sea surface temperatures in the North Atlantic. In contrast, Northeast U.S. shelf surface temperatures were more variable, with near record highs in winter and near average conditions in other seasons.

#### Regional/Coastal phenomena
There was a documented [die-off of scallops](https://noaa-edab.github.io/catalog/observation_synthesis.html) in the Mid-Atlantic Elephant Trunk regions between the 2022 and 2023 surveys. In 2022, Elephant Trunk experienced [stressful temperatures](https://noaa-edab.github.io/catalog/thermal_habitat_persistence.html) for scallops (17 - 19 ℃) for an average of 30 days (Fig. \ref{fig:elephantrunk}), but ongoing research is being conducted to identify contributing factors. A  fish and shellfish mortality event was observed in coastal New Jersey linked to [hypoxia and ocean acidification](https://noaa-edab.github.io/catalog/ocean_acidification.html) (Fig. \ref{fig:NJhypoxia}).

```{r elephantrunk, out.width='50%', fig.cap="The number of days where bottom temperature was between 17 and 19 °C in each GLORYS grid cell for 2022. The gray lines show the sea scallop estimation areas, with the Elephant Trunk region highlighted in black lines."}
magick::image_read("https://github.com/NOAA-EDAB/ecodata/raw/dev/data-raw/workshop/images/bottom_temp_threshold_17_19_Elephant_Trunk.png")

```

```{r NJhypoxia, out.width='50%', fig.cap = "Mission tracks of three gliders (left) deployed off the coast of New Jersey in August and September of 2023. Locations of hypoxic levels of dissolved oxygen (magenta; < 3 mg/liter) and low aragonite saturation state (cyan; < 1) measured along the glider mission tracks and locations of reported fish, lobster, and/or crab mortalities (red X)."}
magick::image_read("https://github.com/NOAA-EDAB/ecodata/raw/dev/data-raw/workshop/images/Figure2-GraceSaba_2024.png")
```

Summer [bottom temperatures](https://noaa-edab.github.io/catalog/bottom_temp_comp.html) in the Gulf of Maine were the warmest on record (since 1959) resulting in the second largest [bottom marine heatwave](https://noaa-edab.github.io/catalog/heatwave_year.html). The heatwave started in February, peaked in May and likely continued beyond August (pending data update). [2023 bottom temperature](https://noaa-edab.github.io/catalog/thermal_habitat_persistence.html) exceeded the 15 ℃ threshold for up to 59 days along the shelf break. 

A wide-spread, long-duration [phytoplankton bloom](https://noaa-edab.github.io/catalog/observation_synthesis.html) of the dinoflagellate *Tripos muelleri* was observed in the GOM and generated chlorophyll concentrations up to ten times greater than average (a record high since 1998) from March to August (Figs. \ref{fig:chl-current-year}, \ref{fig:GOMbloom}). The bloom severely reduced water clarity, impacting harpoon fishing and likely affecting visual predators. Despite *Tripos* being a similar size to typical large phytoplankton (diatoms), this extra production was not grazed nor did it sink to the bottom. The specific drivers of the bloom and implications to the food web are still under investigation. 

```{r chl-current-year, fig.cap="2023 median weekly chlorophyll concentrations (green line) with standard deviation 1998-2023 (gray shading).", fig.width=7, fig.asp=.35}

a <- ecodata::plot_chl_pp(report = "MidAtlantic") + ggplot2::ggtitle("Mid Atlantic 2024")

b <- ecodata::plot_chl_pp(report = "NewEngland") + ggplot2::ggtitle("New England 2024")

a + b + patchwork::plot_layout(widths = c(2.6, 5.2))

```

```{r GOMbloom, out.width='50%', fig.cap="June 2023 chlorophyll anomaly shown as the ratio of the June 2023 average compared to climatological (1998-2023) June average. The black line depicts the 100 m isobath."}
magick::image_read("https://github.com/NOAA-EDAB/ecodata/raw/dev/data-raw/workshop/images/2023-GOMbloom-internalreview.png")
```

In Chesapeake Bay, [hypoxia conditions](https://noaa-edab.github.io/catalog/ches_bay_synthesis.html) were the lowest on record (since 1995), creating more suitable habitat for multiple fin fish and benthic species. Cooler [Chesapeake Bay water temperatures](https://noaa-edab.github.io/catalog/ches_bay_sst.html) paired with low hypoxia in the summer suggest conditions that season were favorable for striped bass. Cooler summer temperatures also support juvenile summer flounder growth.  However, warmer winter and spring water temperatures in the Chesapeake Bay, along with other environmental factors (such as low flow), may have played a role in low production of juvenile striped bass in 2023. Higher-than-average [salinity](https://noaa-edab.github.io/catalog/ch_bay_sal.html) across the Bay was likely driven by low precipitation and increased the area of available habitat for species such as croaker, spot, menhaden, and red drum, while restricting habitat area for invasive blue catfish.

#### Shelf-wide Phenomena
The [Gulf Stream](https://noaa-edab.github.io/catalog/observation_synthesis.html) was highly variable in 2023, with northward shifts intermittently throughout the year and a more notable prolonged shift north along the continental shelf break in the southern Mid-Atlantic in the fall (Fig. \ref{fig:gulfstream23}). This shift severely constricted the Slope Sea (the waters between the Gulf Stream and continental shelf), inhibited warm core ring formation and interactions, resulted in unusually warm and salty surface waters, and strong northeastward currents in the southern Mid-Atlantic. Intermittent warm waters like this can be threats to temperature sensitive species, especially species at the southern end of their range or are not mobile (e.g. scallops), while also providing suitable habitat for more southern species. 
```{r gulfstream23, out.width='50%', fig.cap="October 8-14, 2023 sea surface temperature average derived from the Advanced Clear Sky Processor for Ocean (ACSPO) SST data. The black line depicts the 100 m isobath and the white line is the mean path of the Gulf Stream."}
magick::image_read("https://github.com/NOAA-EDAB/ecodata/raw/dev/data-raw/workshop/images/MAP_SST_Gulf_Stream_2024.png")

```


While the total number of [warm core rings](https://noaa-edab.github.io/catalog/wcr.html) in 2023 (18) was below the decadal average (31), there were a few notable events. A large early season ring pulled continental shelf water into the Slope Sea. Events like these can create biological hotspots, aggregating multiple species in small areas, increasing bycatch risks, and marine mammal shipstrike risks. In spring 2023, concentrations of North Atlantic right whales, humpback whales, basking sharks, and other large baleen whales were observed feeding near the edge of warm core rings near the shelf break. 

Multiple fall 2023 tropical and coastal storms caused several flash flood events, above average coastal water levels, strong winds and high rainfall totals throughout the Northeast. These storms may be related to the shift from 2020-2022 La Niña conditions to strong El Niño conditions in late spring 2023. El Niño winters are associated with more frequent East Coast storms, which can result in increased risk of coastal flooding,  increased freshwater runoff into the coastal ocean, and delayed spring transition from a well mixed water column to stratified. In estuaries, increased freshwater flow decreases salinity, reduces the amount of suitable habitat for juvenile marine fish, and is related to increased hypoxia (low oxygen). However, precipitation is not uniform throughout the Northeast U.S., and [Chesapeake Bay 2023 conditions](https://noaa-edab.github.io/catalog/ch_bay_sal.html)  did not align with El Nino expectations.  The current El Niño is expected to weaken by spring 2024. 

# Contributors

**Editors** (NOAA NMFS Northeast Fisheries Science Center, NEFSC): Sarah Gaichas, Joseph Caracappa, Andy Beet, Brandon Beltz, Geret DePiper, Kimberly Hyde, Scott Large, Sean Lucey, Laurel Smith

**Contributors** (NEFSC unless otherwise noted): Kimberly Bastille, Aaron Beaver (Anchor QEA), Andy Beet, Brandon Beltz, Ruth Boettcher (Virginia Department of Game and Inland Fisheries), Mandy Bromilow and CJ Pellerin (NOAA Chesapeake Bay Office), Joseph Caracappa, Baoshan Chen (Stony Brook University), Doug Christel (GARFO), Patricia Clay, Lisa Colburn, Jennifer Cudney and Tobey Curtis (NMFS Atlantic HMS Management Division), Geret DePiper, Dan Dorfman (NOAA-NOS-NCCOS), Hubert du Pontavice, Emily Farr and Grace Roskar (NMFS Office of Habitat Conservation), Michael Fogarty, Paula Fratantoni, Kevin Friedland, Marjy Friedrichs (VIMS), Sarah Gaichas, Ben Galuardi (GAFRO), Avijit Gangopadhyay (School for Marine Science and Technology, University of Massachusetts Dartmouth), James Gartland (Virginia Institute of Marine Science), Lori Garzio (Rutgers University), Glen Gawarkiewicz (Woods Hole Oceanographic Institution), Sean Hardison, Dvora Hart, Kimberly Hyde, John Kocik, Steve Kress and Don Lyons (National Audubon Society’s Seabird Restoration Program), Young-Oh Kwon and Zhuomin Chen (Woods Hole Oceanographic Institution), Andrew Lipsky, Sean Lucey, Chris Melrose, Shannon Meseck, Ryan Morse, Ray Mroch (SEFSC), Brandon Muffley (MAFMC), Kimberly Murray, David Moe Nelson (NCCOS), Janet Nye (University of North Carolina at Chapel Hill), Chris Orphanides, Richard Pace,Debi Palka, Tom Parham (Maryland DNR), Charles Perretti, Jeffrey Runge (U Maine), Grace Saba and Emily Slesinger (Rutgers University), Vincent Saba, Sarah Salois, Chris Schillaci (GARFO), Amy Schueller (SEFSC), Teresa Schwemmer (Stony Brook University), Dave Secor (CBL), Angela Silva, Adrienne Silver (UMass/SMAST), Laurel Smith, Talya tenBrink (GARFO), Abigail Tyrell, Bruce Vogt (NOAA Chesapeake Bay Office), Ron Vogel (University of Maryland Cooperative Institute for Satellite Earth System Studies and NOAA/NESDIS Center for Satellite Applications and Research), John Walden, Harvey Walsh, Changhua Weng, Dave Wilcox (VIMS), Timothy White (Environmental Studies Program, BOEM), Mark Wuenschel, Qian Zhang (U Maryland)

# Document Orientation

The figure format is illustrated in Fig \ref{fig:docformat}a. Trend lines are shown when slope is significantly different from 0 at the p < 0.05 level. An orange line signifies an overall positive trend, and purple signifies a negative trend. To minimize bias introduced by small sample size, no trend is fit for < 30 year time series. Dashed lines represent mean values of time series unless the indicator is an anomaly, in which case the dashed line is equal to 0. Shaded regions indicate the past ten years. If there are no new data for 2022, the shaded region will still cover this time period. The spatial scale of indicators is either coastwide, Mid-Atlantic states (New York, New Jersey, Delaware, Maryland, Virginia, North Carolina), or at the Mid-Atlantic Bight (MAB) Ecosystem Production Unit (EPU, Fig. \ref{fig:docformat}b) level.

```{r docformat, fig.cap = "Document orientation. a. Key to figures. b.The Northeast Large Marine Ecosystem.",  fig.width = 8, fig.asp = 0.3, ft.arraystretch = 1}
#fig.subcap= c('Key to figures.', 'The Northeast Large Marine Ecosystem.'), out.width = '.49\\linewidth', fig.show="hold" 
# Figure orientation subfigure
m <- 0.1
x <- 1985:2022
y <-  m*x + rnorm(30, sd = 0.35)

data <- data.frame(x = x,
                  y = y)

#Define constants for figure plot
x.shade.max <- max(x)
x.shade.min <- x.shade.max - 9 
hline = mean(y)

setup <- ecodata::plot_setup(shadedRegion = c(2014,2023), report = "MidAtlantic")

#Plot series with trend 
psample <- ggplot2::ggplot(data = data,aes(x = x, y = y)) +
  #Highlight last ten years
  annotate("rect", fill = setup$shade.fill, alpha = setup$shade.alpha,
      xmin = setup$x.shade.min , xmax = setup$x.shade.max,
      ymin = -Inf, ymax = Inf) +
  geom_point(size = setup$pcex) +
  scale_color_manual(aesthetics = "color")+
  guides(color = FALSE) +
  geom_hline(aes(yintercept = hline),
              size = setup$hline.size,
             alpha = setup$hline.alpha,
             linetype = setup$hline.lty)+
  geom_line() +
  geom_gls(aes(x = x, y = y)) +
  scale_y_continuous(labels = function(l){trans = l / 1000})+
  scale_x_continuous(breaks = seq(1985, 2015, by = 5), expand = c(0.01, 0.01)) +
  ylab(expression("Invented Index, 10"^3*"widgets")) +
  xlab("Time") +
  labs(tag = "a")  +
  theme_ts()


# EPU map subfigure
# Set lat/lon window for map
xmin = -78
xmax = -65
ymin = 36
ymax = 45
xlims <- c(xmin, xmax)
ylims <- c(ymin, ymax)

#CRS
crs <- "+proj=longlat +lat_1=35 +lat_2=45 +lat_0=40 +lon_0=-77 +x_0=0 +y_0=0 +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"

#Specify data frame with lat/lon locations for labels
epu_labels <- data.frame(EPU = c("Mid-Atlantic\n Bight",
                                 "Gulf of Maine",
                                 "Georges Bank"),
                         latitude = c(40,42.85,41),
                         longitude = c(-72.7,-69,-68.5))

#Map of NE LME
epumap <- ggplot() +
  geom_sf(data = coast, size = setup$map.lwd) +
  geom_sf(data = epu_sf, fill = "transparent", size = setup$map.lwd) +
  coord_sf(crs = setup$crs, xlim = setup$xlims, ylim = setup$ylims) +
  geom_text(data = epu_labels, aes(x = longitude,
                                    y = latitude,
                                    label = EPU),
            size = 1.7)+
  theme_map() +
  scale_x_continuous(breaks = seq(-78, -65, by = 4), expand = c(0.01, 0.01)) +
  xlab("Longitude") +
  ylab("Latitude") +
  labs(tag = "b")

#epumap <- magick::image_read("images/EPU_Designations_Map.jpg")

psample + plot_spacer() + epumap + plot_layout(ncol = 3, widths=c(3,.5,4.5))

# psample
# 
# knitr::include_graphics("images/EPU_Designations_Map.jpg")
# magick::image_read("images/EPU_Designations_Map.jpg")
# epumap <- magick::image_read("images/EPU_Designations_Map.jpg") %>%
#  magick::image_scale("10%")
# epumap <- magick::image_scale(epumap, "50%")

# epumap

```

Fish and invertebrates are aggregated into similar feeding categories (Table \ref{tab:species-groupings}) to evaluate ecosystem level trends in predators and prey. 

```{r species-groupings, warning=F}

# new table with all species listed by management entity
df <- ecodata::species_groupings %>%
  dplyr::select(SOE.24, COMNAME, Fed.Managed) %>%
  filter(SOE.24 != "Other") %>%
  distinct() %>%
  group_by(SOE.24, Fed.Managed) %>%
  summarize_all(funs(paste(na.omit(.), collapse = ", "))) %>%
  spread(Fed.Managed, COMNAME) %>%
  arrange(factor(SOE.24, levels = c("Apex Predator", "Piscivore", "Planktivore", "Benthivore", "Benthos")))
df<-df[c(1,3,2,4,5)] %>%
  mutate_all(tolower)


# knitr::kable(df, booktabs = TRUE, caption = 'Feeding guilds and management bodies.', 
#              col.names = c("Guild", "MAFMC", "Joint", "NEFMC", "State or Other")) %>%
#   kable_styling(font_size=10, latex_options=c("repeat_header", "scale_down", "hold_position")) %>%
#   row_spec(0,bold=TRUE) %>%
#   column_spec(1, width="2cm") %>%
#   column_spec(2, width="4cm") %>%
#   column_spec(3, width="2cm") %>%
#   column_spec(4, width="5cm") %>%
#   column_spec(5, width="6cm") %>%
#   #column_spec(3, width="7.5cm") #%>%
#   collapse_rows(columns = 1, latex_hline = "major", valign = "middle")

flextable::flextable(df) %>%
  flextable::set_caption('Feeding guilds and management bodies.') %>%
  flextable::fontsize(size=8, part = "all") %>%
  flextable::set_header_labels(SOE.24 = "Guild",
                               MAFMC = "MAFMC",
                               JOINT = "Joint",
                               NEFMC = "NEFMC",
                               "<NA>" = "State or Other") %>%
  flextable::width(width = c(1,1,1,1,3))

  #flextable::fit_to_width(max_width = 20)
  #flextable::autofit()
  #flextable::fontsize(i = NULL, j = NULL, size = 10, part = "body") %>%


```