In hmorzaria/pssalmonsurvival: Sensitivity analysis of salmon survival in Puget Sound

Package for the article Effect of multiple pressures on early marine survival of juvenile salmon in Puget Sound

Hem Nalini Morzaria-Luna"," I.C. Kaplan, C.J. Harvey, M. Schmidt, E.A. Fulton, R. Girardin, and P. MacCready

install.packages("rnaturalearthhires", repos = "http://packages.ropensci.org", type = "source")
library(pssalmonsurvival)

Map of model extent in Puget Sound

#https://www.ngdc.noaa.gov/mgg/shorelines/ A Global Self-consistent, Hierarchical, High-resolution Geography Database
#system("wget https://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/gshhg-shp-2.3.7.zip")

file.name <- "amps_model_map.png"

make_map(file.name)

Plot food web

#Plot model food web

data("ppreymatrix")
plot.name <- "ps_foodweb.png"

plot_foodweb(ppreymatrix, plot.name)

Plot salmon time series

salmon.abundance <- readxl::read_xlsx(here::here("modelfiles","Salmon_timeline_Losee.xlsx"), sheet = 1) %>% 
 # dplyr::rename("Salmon" = `Salmon Sum`) %>% 
  dplyr::select(-`Salmon Sum`,-Steelhead) %>% 
  tidyr::pivot_longer(cols=Chinook:Sockeye,names_to = "salmon_group",values_to = "biomass_mt")

plot_salmon_abundance(salmon.abundance)

Model ensemble was developed with alternate parametrizations. Plot biomass trajectories based on simulations that run the model ensemble 30 years under constant fishing rates. We only use 6 model variants, 2:5 and 7:8 Model biomass comes from 30-year runs of the base models. There are 8 model variants that differ in invertebrate growth rate and vertebrate density dependence

#biomass comparison between model variants, base runs
data("ensemblebiomass")

plot_ensemblebiomass(ensemblebiomass)

Plot salmon survival for model ensemble for base scenarios, survival is defined as the proportion of age 1 salmon that survive to age 5, the cohort is lagged over time. Also summarises base scenario which get used later to estimate relative survival.

#numbers at age from base runs
data("ensemblenumbersage")
data("salmongroups")
#plotmodels <- c(1,6) # eliminated model versions 1 & 6

plot_ensemble_survival(ensemblenumbersage, salmongroups)

Turn eval=True to make forcings for model simulations, by scenario.

#groups for each scenario
sim.groups <- list(c("Pinniped", "California_sea_lion", "Harbor_seal"), 
                c("PinkSY_Fish", "ChumFSY_Fish"), 
                c("Herring_PS", "Herring_Cherry"),
                "Porpoise",
                "Spinydog_Fish", 
                c("Pisc_Seabird","NonPisc_Seabird"),
                "Gel_Zoo", 
                c("ChinookY_Fish","ChinookSY_Fish","CohoHY_Fish","ChumHSY_Fish"),
                #chinook hatchery are repeated twice for a scenario with just Chinook
                c("ChinookY_Fish_2","ChinookSY_Fish_2"))

data("functionalgroups")
data("salmongroups")

ncfile = "AMPS.nc"
num.forcing.years = 50
bash.file = "amps_cal.sh"
force.prm = "PugetSound_force.prm"
func.groups = functionalgroups

#scenario 20% decrease

rate.multiplier = 0.8

lapply(sim.groups, make_forcing, ncfile, func.groups, num.forcing.years, rate.multiplier, bash.file, force.prm, salmongroups)

#scenario 20% increase

rate.multiplier = 1

lapply(sim.groups, make_forcing, ncfile, func.groups, num.forcing.years, rate.multiplier, bash.file, force.prm, salmongroups)

rate.multiplier = 1.2

lapply(sim.groups, make_forcing, ncfile, func.groups, num.forcing.years, rate.multiplier, bash.file, force.prm, salmongroups)

Plot salmon survival for scenario results, survival is defined as the proportion of age 1 salmon that survive to age 5, the cohort is lagged over time.

#base case survival for each model variant created above

base.survival <- readr::read_csv(here::here("modelfiles","base_survival.csv")) %>% 
  tidyr::drop_na() %>%
  dplyr::mutate(max_year = max(year_no)) %>% 
  dplyr::filter(year_no<=(max_year-3)) %>% 
  dplyr::filter(year_no==max(year_no)) %>% 
  dplyr::mutate(scenario_name = dplyr::if_else(scenario_name=="salmon competition","wild pink and chum salmon competition",
                                                 dplyr::if_else(scenario_name=="mammal predation","pinniped predation",
dplyr::if_else(scenario_name=="seabirds predation","seabird predation", scenario_name)))) 



#numbers at age from base runs
data("ensemblenumbersagescenarios")
data("salmongroups")
data("salmonbybasin")
data("indsalmoneffect")


#plotmodels <- c(1,6) # eliminated model versions 1 & 6

plot_ensemble_survival_scenarios(ensemblenumbersagescenarios, salmongroups, plotmodels, base.survival, salmonbybasin, indsalmoneffect)

Plot survival as time series

data("ensemblenumbersagescenarios")
data("salmongroups")
data("indsalmoneffect")
#data("salmoneffect")

survival.rates <- readr::read_csv(here::here("modelfiles","survival_rates.csv"))


plot_ensemble_survival_scenarios_timeseries(ensemblenumbersagescenarios, salmongroups, indsalmoneffect)

#Not run
#Using trajectory analysis, but it was impossible to interpret
data("ensemblenumbersagescenarios")
data("salmongroups")

make_trajectories(ensemblenumbersagescenarios, salmongroups)

Extract base survival for cumulative impact scenarios

#numbers at age from base runs
data("ensemblenumberscum")
data("salmongroups")
#plotmodels <- c(1,6) # eliminated model versions 1 & 6

extract_base_cum_survival(ensemblenumberscum, salmongroups)

Plot survival in cumulative impact scenarios

base.cum.survival <- readr::read_csv(here::here("modelfiles","base_cum_survival.csv")) %>% 
  tidyr::drop_na() %>%
  dplyr::mutate(max_year = max(year_no)) %>% 
  dplyr::filter(year_no<=(max_year-3)) %>% 
  dplyr::filter(year_no==max(year_no)) 


#numbers at age from base runs
data("ensemblenumberscum")
data("salmongroups")
data("salmonbybasin")
data("salmoneffect")

#plotmodels <- c(1,6) # eliminated model versions 1 & 6

plot_ensemble_cum_survival_scenarios(ensemblenumberscum, salmongroups, base.cum.survival, salmonbybasin, salmoneffect)

Plot interaction effect

ensemble.survival <- readr::read_csv(here::here("modelfiles","ensemble_survival.csv")) 
ensemble.cum.survival <- readr::read_csv(here::here("modelfiles","ensemble_cum_survival.csv")) 

data("scenariocategories")
data("salmonbybasin")


plot_interaction_effect(ensemble.survival, ensemble.cum.survival, scenariocategories, salmonbybasin, thiscutoff = 125)

Plot biomass of higher trophic level guilds in individual scenarios

#only run once
#data("ppreymatrix")
#data("functionalgroups")

#get salmon predators from food web matrix
#get_salmonpredators(ppreymatrix, functionalgroups)

data("ensemblebiomass")
data("predgroups")

plot_predators(ensemblebiomass, predgroups, thiscutoff = 50)

plot biomass of higher trophic level guilds in cumulative scenarios

data("ensemblebiomasscum")
data("predgroups")

plot_predatorscum(ensemblebiomasscum, predgroups, thiscutoff = 20)

Explore predation using dietcheck for individual scenarios

data("scenariocategories")
data("salmonbybasin")
data("dietcheck")

plot_dietcheck(dietcheck, indsalmoneffect, predgroups, thiscutoff=30)

Explore predation using dietcheck for cumulative scenarios

data("salmoneffect")
data("salmonbybasin")
data("dietcheckcum")

plot_dietcheckcum(dietcheckcum, salmoneffect, predgroups, thiscutoff = 100)

Plot weight-at-age for individual scenarios

plot_watage(ensemblewageind, salmonbybasin, salmongroups, thiscutoff=100)

Plot weight-at-age for cumulative scenarios

plot_watage_cum(ensemblewagecum, salmonbybasin, salmongroups, salmoneffect, thiscutoff=20)

Cite R and packages

devtools::session_info()

c("broom","ggplot2","ggthemes","ggspatial","ggsci","rgdal","rnaturalearth","sf","sp","magrittr","ggforce","dplyr","here","ncdf4","readr",
  "stringr","RNetCDF","scales","GGally","network","sna","RColorBrewer","grDevices","colorRamps","colorspace","tidyr","Hmisc","forcats",
  "ecotraj","ggnewscale","paletteer","wesanderson","purrr") %>%
  purrr::map(citation) %>%
  print(style = "text")

hmorzaria/pssalmonsurvival documentation built on Aug. 20, 2023, 12:27 a.m.