updated script that creates the benthic time series

fitting/benthic_time.R

@@ -1,72 +1,72 @@
 # Benthic time series
+# Author: M.T.Grezlik
+# [email protected]
 # Following Sarah Weisberg's approach for the GOM
 # https://github.com/SarahJWeisberg/GOM-Rpath/blob/main/fitting/benthic_time.R
 
-#pulling in mega and macrobenthos time series for forcing
+# Load packages and data --------------------------------
 
 #load packages
 library(dplyr)
 library(tidyr)
 library(here)
+library(ggplot2)
 
 #Load balanced model
 load(here("data/alternate.GB.bal.rda"))
 load(here("data/alternate.GB.params.bal.rda"))
 
 #macrobenthos 
-url<-"https://github.com/NOAA-EDAB/Rpathdata/blob/dd034d1573f79ce011c01054bdc017e241e7857e/data-raw/fallmacrobenthosindex.rds?raw=true"
-download.file(url,destfile = "data/fallmacrobenthosindex.rds")
-macro <- readRDS("data/fallmacrobenthosindex.rds")
+macro <- readRDS(url('https://github.com/NOAA-EDAB/Rpathdata/blob/dd034d1573f79ce011c01054bdc017e241e7857e/data-raw/fallmacrobenthosindex.rds?raw=true'))
 
-load(url(macro_url))
 
 #megabenthos
-url<-"https://github.com/NOAA-EDAB/Rpathdata/blob/main/data-raw/fallmegabenthosindex.rds?raw=true"
-download.file(url,destfile = "data/fallmegabenthosindex.rds")
-mega <- readRDS("data/fallmegabenthosindex.rds")
+mega <-readRDS(url("https://github.com/NOAA-EDAB/Rpathdata/blob/main/data-raw/fallmegabenthosindex.rds?raw=true"))
 
-#filter for GOM only, remove unneeded column, relabel biomass, SE for simplicity
-macro_gom <- macro %>% 
-  filter(EPU == "GOM") %>% 
-  select(-Units) %>%
+# Data processing --------------------------------
+
+#filter for GB only, remove unneeded column, relabel biomass, SE for simplicity
+macro_gb <- macro |>  
+  filter(EPU == "GB") |> 
+  select(-Units)  |> 
   mutate(Var = ifelse(Var == "Fall Macrobenthos Biomass Index Estimate","biomass","SE"))
-mega_gom <- mega %>% 
-  filter(EPU == "GOM") %>% 
-  select(-Units) %>%
+mega_gb <- mega  |>  
+  filter(EPU == "GB") |> 
+  select(-Units) |> 
   mutate(Var = ifelse(Var == "Fall Megabenthos Biomass Index Estimate","biomass","SE"))
 
 #calculate 1980-85 mean
-macro_start <- macro_gom %>% 
-  filter(Time < 1986) %>% 
-  group_by(Var) %>%
+macro_start <- macro_gb |>  
+  filter(Time < 1986) |> 
+  group_by(Var) |> 
   summarise(mean_start = mean(Value))
-mega_start <- mega_gom %>% 
-  filter(Time < 1986) %>% 
-  group_by(Var) %>%
+mega_start <- mega_gb |>  
+  filter(Time < 1986) |>  
+  group_by(Var) |> 
   summarise(mean_start = mean(Value))
 
 #merge with time series, calculate anomalies
-macro_gom <- left_join(macro_gom,macro_start,by="Var")
-macro_gom <- macro_gom %>% mutate(anom = (Value-mean_start)/mean_start)
-mega_gom <- left_join(mega_gom,mega_start,by="Var")
-mega_gom <- mega_gom %>% mutate(anom = (Value-mean_start)/mean_start)
+macro_gb <- left_join(macro_gb,macro_start,by="Var")
+macro_gb <- macro_gb |>  mutate(anom = (Value-mean_start)/mean_start)
+mega_gb <- left_join(mega_gb,mega_start,by="Var")
+mega_gb <- mega_gb |>  mutate(anom = (Value-mean_start)/mean_start)
 
 #generate forcing time series
-macro<-GOM.params$model[Group == "Macrobenthos",Biomass]
-mega<-GOM.params$model[Group == "Megabenthos",Biomass]
+macro<-alternate.GB.params.bal$model[Group == "Macrobenthos",Biomass]
+mega<-alternate.GB.params.bal$model[Group == "Megabenthos",Biomass]
 
-macro_time<-macro_gom %>% 
-  filter(Var == "biomass" & Time > 1984) %>% 
-  select(Time,anom) %>%
+macro_time<-macro_gb |>  
+  filter(Var == "biomass" & Time > 1984) |>  
+  select(Time,anom) |> 
   mutate(biomass = macro*(anom+1))
 
-mega_time<-mega_gom %>% 
-  filter(Var == "biomass" & Time > 1984) %>% 
-  select(Time,anom) %>%
+mega_time<-mega_gb |>  
+  filter(Var == "biomass" & Time > 1984)  |>  
+  select(Time,anom)  |> 
   mutate(biomass = macro*(anom+1))
 
 #visualize
-# mega_gom %>% filter(Var == "biomass" & Time > 1984) %>%
+# mega_gb |>  filter(Var == "biomass" & Time > 1984)  |> 
 #   ggplot(aes(x=Time,y=anom))+
 #   geom_line()+
 #   geom_point()+