NewCatchWorkflow.Rmd

---
title: "New catch functions from CATCH.nc file"
author: "Sarah Gaichas"
date: "`r format(Sys.time(), '%d %B %Y')`"
output:
  html_document:
    code_fold: hide
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r message=FALSE, warning=FALSE}
library(here)
library(stringr)
library(magrittr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(atlantisom)
```


## Intro

Total catch in tons has been pulled from the Atlantis `Catch.txt` output file from 2019-2023 due to difficulties scaling catch in numbers from the `CATCH.nc` file to total weight. See comparisons [here](https://sgaichas.github.io/poseidon-dev/TestCatchNOBA.html). 

For testing multispecies models, there is a need for subannual catch information, which is not provided in the annual `Catch.txt` output. In addition, there is no polygon information in the `Catch.txt` output, so spatial subsetting is not possible.

In order to get catch in tons subannually and by area, we can use the fleet-specific catch information in the `CATCH.nc` file. These outputs are reported in tons by species aggregated over all age classes for each defined fleet. Variable names in the `CATCH.nc` file for these outputs follow the format "XXX_Catch_FCN" and "XXX_Discard_FCN" where XXX is the 3 letter group Code defined in the Atlantis `groups.csv` file and N is the integer Index for the fishery defined in the Atlantis `fisheries.csv` file. 

## Methods

Workflow

1.  Modify `atlantisom::load_nc` to get catch in tons from different `CATCH.nc` outputs. New function because current function already gets catch in numbers, would be complicated to use the same "Catch" variable to get tons. Name new function `atlantisom::load_nc_catchtons`

2.  Test new `atlantisom::load_nc_catchtons` with NOBA, NEUS, CC by summing output over polygons and fleets to each year and comparing with annual output of current `atlantisom::load_catch` based on `Catch.txt`

3.  Incorporate `atlantisom::load_nc_catchtons` into `atlantisom::run_truth`. Add a new catchtons object keeping polygon, timestep, and fleet information to the output list. Consider removing the `Catch.txt` output in this function, currently called `catch_all` but not used? Or keep it, since catchtons can be summed to catch_all this could be used as an internal test as was started in the run_truth function. 

4.  Adjust `atlantisom::om_init` for the new catchtons output. Should we keep reading in the `Catch.txt` here? It may be ok to keep it just to keep existing functions for aggregate annual catch working.

5.  Adjust `atlantisom::om_species` for the new catchtons output. 

6.  Adjust `atlantisom::om_index` for the new catchtons output. Create a new index that will be consistent across the subannual and annual datasets: this means applying the cv at the level of reporting? So possibly by fleet and timestep, then summing to annual? This involves including and updating two functions: `atlantisom::create_fishery_subset` parallel to the `atlantisom::create_survey` function to subset areas, currently used only on fishery comps but can now be applied to the index as well, and updating another `atlantisom::sample_fishery_totalcatch` for the cv application. This will save perhaps one new output and update one existing output, or just one output? Will we have the same current `censusfishCatch.rds` output that is annually aggregated plus a disaggregated version? Or have a single output with the same name `censusfishCatch.rds` that is disaggregated, to be summed separately in mskeyrun?

7. Adjust `mskeyrun` functions creating aggregated and subannual catch datasets. Do we need a fleet specific index?


### 1. New `atlantisom::load_nc_catchtons` function

Here is the new function, tested with snippets in mskeyrun sarah_wgsamsim branch and implemented in atlantisom dev branch

```{r, code = readLines("https://raw.githubusercontent.com/r4atlantis/atlantisom/dev/R/load_nc_catchtons.R"), eval=F}
```


### 2. Test `load_nc_catchtons` with the NOBA, NEUS, CC model outputs

For all plots, the total catch derived from Catch.txt is considered correct, and is represented by the black line. The total catch summed across polygon, subannual timestep and fleet is in blue. 

Only NOBA first, this one is model sacc30 with climate, different from sacc38 in `mskeyrun` without climate: 

```{r, message=FALSE, warning=FALSE}

# define for each model

## NOBA (different from mskeyrun version)

atlmod <- here("config/NOBAsacc30Config.R")

select_groups <- c("Long_rough_dab",
                   "Green_halibut",
                   "Mackerel",
                   "Haddock",
                   "Saithe",
                   "Redfish",
                   "Blue_whiting",
                   "Norwegian_ssh",
                   "North_atl_cod",
                   "Polar_cod",
                   "Capelin")

stepperyr <- 5



## NEUS



##### aggregate over polygons to fleets #######
# make this a function

comparecatch <- function(atlmod, select_groups, stepperyr) {
  source(atlmod, local = TRUE)
  
  dir <- d.name
  nc_catch <- paste0(scenario.name, 'CATCH.nc')
  file_fgs <- functional.groups.file
  file_init <- initial.conditions.file
  file_fish <- fisheries.file
  
  # Get the boundary boxes
  allboxes <- atlantisom::load_box(dir = d.name, file_bgm = box.file)
  boxes <- atlantisom::get_boundary(allboxes)
  
  
  # Read in information
  # Read in the functional groups csv since that is used by many functions
  fgs <- load_fgs(dir = dir, file_fgs = file_fgs)
  # Read in the biomass pools
  bps <- load_bps(dir = dir,
                  fgs = file_fgs,
                  file_init = file_init)
  
  catchtons <- load_nc_catchtons(
    dir = dir,
    file_nc = nc_catch,
    file_fish = file_fish,
    bps = bps,
    fgs = fgs,
    select_groups = select_groups,
    select_variable = "Catch",
    check_acronyms = TRUE,
    bboxes = boxes
  )
  #if(verbose) message("Catch tons read in.")
  
  # result is output of load_nc_catchtons
  
  aggcatchtons <- catchtons %>%
    dplyr::select(-agecl) %>%
    dplyr::group_by(species, fleet, time) %>%
    dplyr::summarise(totcatch = sum(atoutput)) %>%
    dplyr::mutate(year = ceiling(time / stepperyr))
  
  # compare annual catch in tons to Catch.txt output
  yearcatchtons <- aggcatchtons %>%
    dplyr::group_by(species, year) %>%
    dplyr::summarise(yeartot = sum(totcatch))
  
  txtCatchtons <- atlantisom::load_catch(d.name, catch.file, fgs) %>%
    dplyr::filter(species %in% select_groups) %>%
    dplyr::mutate(year = time / 365) %>%
    dplyr::left_join(dplyr::select(fgs, Code, Name), by = c("species" = "Name"))
  
  comparetons <- txtCatchtons %>%
    dplyr::left_join(yearcatchtons, by = c("Code" = "species", "year" = "year")) %>%
    ggplot2::ggplot() +
    geom_line(aes(year, atoutput)) +
    geom_point(aes(year, yeartot), color = "blue", alpha = 0.2) +
    facet_wrap( ~ species, scales = "free_y")
  
  # looks like a match!
  comparetons
  
}

comparecatch(atlmod = atlmod, 
             select_groups = select_groups,
             stepperyr = stepperyr)

```

Now try the Isaac's recent CC model:

```{r, message=FALSE, warning=FALSE}
## CC

atlmod <- here("config/CCConfigSep22.R")

select_groups <- c("Arrowtooth_flounder",
                   "Petrale_sole",
                   "Pacific_sardine",
                   "Anchovy",
                   "Herring",
                   "Pacific_Ocean_Perch",
                   "Bocaccio_rockfish",
                   "Yelloweye_rockfish",
                   "Mesopel_M_Fish", #Pacific hake
                   "Mesopel_N_Fish", #Sablefish
                   "Demersal_P_Fish") #Dover sole

stepperyr <- 5

comparecatch(atlmod = atlmod, 
             select_groups = select_groups,
             stepperyr = stepperyr)
```

And finally the NEUS model, thanks to Joe for output files: 

```{r, message=FALSE, warning=FALSE}
## NEUS
atlmod <- here("config/NEUStestConfigMay23.R")

select_groups <- c("Summerflounder",
                   "Winterflounder",
                   "Mackerel",
                   "Herring",
                   "Cod",
                   "Monkfish",
                   "Haddock",
                   "Spiny_Dogfish",
                   "Winter_Skate",
                   "Yellowtail_Flounder",
                   "Silver_Hake")

stepperyr <- 5

comparecatch(atlmod = atlmod, 
             select_groups = select_groups,
             stepperyr = stepperyr)

```

Looks like we have a match to Catch.txt output for all systems using the same function and inputs from the CATCH.nc file. However, it looks like there is a final value in Catch.txt that is not present, or not being summed correctly in the CATCH.nc output in the CC and NEUS models. Worth investigating...

### 3. Incorporate new `load_nc_catchtons` function into `atlantisom::run_truth`

Added catchtons and disctons (same function with `select_variable = "Discards"`) to `run_truth` in the atlantisom dev branch. At present, catchtons and disctons are in the run_truth output and catch_all which was originally from the `Catch.txt` file have been removed from the output. Sections of the function that read in both the `Catch.txt` and the `CatchPerFishery.txt` files have been commented out. They could be reinstated if needed, but `om_init` reads in `Catch.txt` already and seems like a good place to read in text files. 

Here is the new function, tested by re-running the code in `SimData.Rmd` in the mskeyrun repo using the NOBA sacc38 outputs. The new runtruth object includes the new catchtons and disctons (0) objects. 

```{r, code = readLines("https://raw.githubusercontent.com/r4atlantis/atlantisom/dev/R/run_truth.R"), eval=F}
```

### 4.  Adjust `atlantisom::om_init` for the new catchtons output

Actually, no adjustment needed for `om_init`. This was run as is to produce the new truth output. I think we keep reading in `Catch.txt` here.

### 5.  Adjust `atlantisom::om_species` for the new catchtons output

Added the catchtons and disctons truth outputs here, subset by code rather than species name, update documentation. Note that code will need to be changed to species prior to applying the fishery sampling functions. 

New function, tested by re-running code in `SimData.Rmd` in mskeyrun as above. The most recent `omlist_ss.rds` file in the mskeyrun simulated data atlantisoutput NOBAsacc38 directory reflects the updated catch in tons at full resolution.

```{r, code = readLines("https://raw.githubusercontent.com/r4atlantis/atlantisom/dev/R/om_species.R"), eval=F}
```

### 6. Adjust `atlantisom::om_index` for the new catchtons output

The complex part. This requires adjusting the `atlantisom::create_fishery_subset` to take catchtons as an input, then feeding this into the index using `atlantisom::sample_fishery_totalcatch` and applying cv at a more disaggregated point, then summing to the annual index so that subannual and annual indices are consistent. Then `om_index` will output the subannual index only, or both? It should probably just make one and users can decide how to aggregate to annual (ie. mskeyrun will provide both a subannual and annual index based on the same `om_index` output).

#### 6a. Update `om_index`

The fishery portion of this wrapper can now have parallel structure to the survey portion. It allows multiple fishery.R config files so the user can get a different index for each fleet if that is desired. Specifications for which polygons are kept for which fleet number would be in the `fishery.R` config file, with an appropriate fleet.name. 

I think we want to aggregate across all fleets for a species for the given set of polygons: for each `fishery.R` file the first step is to select fleets identified in `fishery.R` and aggregate across those fleets for the polygons specified in `fishery.R`, keeping species separated out--this is done in `create_fishery_subset` below.

Note that allowing for multiple `fishery.R` files for each fleet means we need to make a parallel change to `om_comps`.

Note that species in here are the codes not the names, so conversion back to names needs to happen somewhere. The fishery.R file now specifies the species codes in fishspp instead of species names (survspp still does that).

Time in this output is *timestep* not days.

Here is the new `om_index`

```{r, code = readLines("https://raw.githubusercontent.com/r4atlantis/atlantisom/dev/R/om_index.R"), eval=F}

```

And an example `fishery.R` file (modified `msfishery.R` file used in ms-keyrun project):

```{r, code = readLines("https://raw.githubusercontent.com/NOAA-EDAB/ms-keyrun/sarah_wgsamsim/data-raw/simulated-data/config/msfishery.R"), eval=F}

```


#### 6b. Update `create_fishery_subset`

Left alone, this function will aggregate across all fleets and polygons given in the input data (since it is not expecting a fleet column). This needs modification to use the fleet column found in catchtons. That way the user can control which fleets are included in the aggregated output. Specifying fleets=NULL means all fleets are aggreagated.

```{r, code = readLines("https://raw.githubusercontent.com/r4atlantis/atlantisom/dev/R/create_fishery_subset.R"), eval=F}

```


#### 6c. Update `sample_fishery_totcatch`

Now this function determines if there are polygons in the input, and aggregates over polygon if there are. It then applies cv at the subannual level for the specified fishery. 

```{r, code = readLines("https://raw.githubusercontent.com/r4atlantis/atlantisom/dev/R/sample_fishery_totcatch.R"), eval=F}

```

#### Testing

Compare the previous total catch time series (annual) with the new one summed to annual and hope the mostly match (both have a cv applied but it is 0.01 for all species so shouldn't matter.)

Read in saved files created previously in ms-keyrun and the one just created with the new functions. This is now comparing the "sampled" fishery data. 

```{r}

mskeyrundir <- "/Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38"

scenario.name <- "nordic_runresults_01"

stepperyr <- 5

fisherycomp <- atlantisom::read_savedfisheries(mskeyrundir, "Catch")

from_nc <- fisherycomp$allfleet[[1]]

from_txt <- fisherycomp$census[[1]]

fgs <- atlantisom::load_fgs(mskeyrundir, "nordic_groups_v04.csv")

  aggcatchtons <- from_nc %>%
    dplyr::select(species, time, atoutput) %>%
    dplyr::mutate(year = ceiling(time / stepperyr)) %>%
    dplyr::group_by(species, year) %>%
    dplyr::summarise(yeartot = sum(atoutput))
  
  txtCatchtons <- from_txt %>%
    dplyr::select(species, time, atoutput) %>%
    dplyr::mutate(year = time / 365) %>%
    dplyr::left_join(dplyr::select(fgs, Code, Name), by = c("species" = "Name"))
  
  comparetons <- txtCatchtons %>%
    dplyr::left_join(aggcatchtons, by = c("Code" = "species", "year" = "year")) %>%
    ggplot2::ggplot() +
    geom_line(aes(year, atoutput)) +
    geom_point(aes(year, yeartot), color = "blue", alpha = 0.2) +
    facet_wrap( ~ species, scales = "free_y")
  
  # looks like a match!
  comparetons


```

We have a match so it seems the workflow gives the same results in aggregate. Lets call the fleet specific and subannual catch good then. It can go into `mskeyrun`.

### 7. Adjust `mskeyrun` functions creating aggregated and subannual catch datasets

See mskeyrun repo. I would like to keep life simple and not have fleet specific subannual catch indices, but since we just did fleet specific indices with the real Georges Bank data I suppose I will have to test this with simulated data.

But maybe not for the WGSAM skill assessment. See below for why I don't quite have a full set of fleet specific outputs.

### 8. Adjust `atlantisom::om_comps` for the new catchtons output

Well this sounded easier than it is. We can have multiple fleets in `om_comps` using a similar structure to surveys. Sadly, the length comps from fisheries cannot be by fleet because they come from the CATCH.nc numbers output, which is by age class and not by fleet. Lengths are estimated using the age class information by numbers, structn, and resn. I don't think we can apply the `calcage2length` function to the true age data because we don't have the growth information from structn and resn at that level.

I can rewrite `om_comps` to use fleet information on all but the lengths. Placeholder code below.

```{r, code = readLines("https://raw.githubusercontent.com/r4atlantis/atlantisom/dev/R/om_comps.R"), eval=F}

```