In HaikunXu/BSE: IATTC's BSE code

knitr::opts_chunk$set(warning = FALSE, message = FALSE)

This example code demonstrates how to compile the purse-seine catch and length composition data for the stock assessment of bigeye tuna in the eastern Pacific Ocean.

library(tidyverse)
library(BSE)

save_dir <- "D:/OneDrive - IATTC/IATTC/2024/SAC15/Assessment/Data/PS/"
yr.end <- 2023

BET.OBJ.Catch.20002023 <- read.csv(paste0(save_dir,"BET.OBJ.Catch.20002023.csv"))
BET.OBJ.Catch.19751999 <- read.csv(paste0(save_dir,"BET.OBJ.Catch.19751999.csv"))

BET.NOA.Catch.20002023 <- read.csv(paste0(save_dir,"BET.NOA.Catch.20002023.csv"))
BET.NOA.Catch.19751999 <- read.csv(paste0(save_dir,"BET.NOA.Catch.19751999.csv"))

BET.DEL.Catch.20002023 <- read.csv(paste0(save_dir,"BET.DEL.Catch.20002023.csv"))
BET.DEL.Catch.19751999 <- read.csv(paste0(save_dir,"BET.DEL.Catch.19751999.csv"))

Year_OBJ <- data.frame(Year = seq(1,(yr.end-1974)*4),
                       Area = rep(c("A1","A2","A3","A4","A5"), each = (yr.end-1974)*4))

BET_OBJ_Catch <- rbind(BET.OBJ.Catch.19751999,BET.OBJ.Catch.20002023) %>% 
  mutate(Year=(year-1975)*4+quarter) %>%
  gather(3:7,key="Area",value="Catch") %>% 
  select(Year,Area,Catch)
BET_OBJ_Catch <- left_join(Year_OBJ,BET_OBJ_Catch) %>%
  mutate(catch=ifelse(is.na(Catch),0,Catch),
         Type="OBJ")

# add discard to the OBJ catches
BET_OBJ_Catch <- left_join(BET_OBJ_Catch, OBJ_discard_rate_BET)
BET_OBJ_Catch <- BET_OBJ_Catch %>%
  mutate(Catch = catch * scaler) %>%
  select(Year, Area, Catch, Type)

# COVID year's catch adjustment based on the CART model
BET_OBJ_Catch$Catch[which(BET_OBJ_Catch$Year %in% 181:184)] <- BET_OBJ_Catch$Catch[which(BET_OBJ_Catch$Year %in% 181:184)] * 0.88
BET_OBJ_Catch$Catch[which(BET_OBJ_Catch$Year %in% 185:188)] <- BET_OBJ_Catch$Catch[which(BET_OBJ_Catch$Year %in% 185:188)] * 0.82

Year_NOA <- data.frame(Year = seq(1,(yr.end-1974)*4),
                       Area = rep(c("A1","A2"), each = (yr.end-1974)*4))
BET_NOA_Catch <- rbind(BET.NOA.Catch.19751999,BET.NOA.Catch.20002023) %>% 
  mutate(Year=(year-1975)*4+quarter) %>%
  gather(3:4,key="Area",value="Catch") %>% 
  select(Year,Area,Catch)
BET_NOA_Catch <- left_join(Year_NOA,BET_NOA_Catch) %>%
  mutate(Catch=ifelse(is.na(Catch),0,Catch),
         Type="NOA")

Year_DEL <- data.frame(Year = seq(1,(yr.end-1974)*4),
                       Area = rep(c("A1","A2"), each = (yr.end-1974)*4))
BET_DEL_Catch <- rbind(BET.DEL.Catch.19751999,BET.DEL.Catch.20002023) %>% 
  mutate(Year=(year-1975)*4+quarter) %>%
  gather(3:4,key="Area",value="Catch") %>% 
  select(Year,Area,Catch)
BET_DEL_Catch <- left_join(Year_DEL,BET_DEL_Catch) %>%
  mutate(Catch=ifelse(is.na(Catch),0,Catch),
         Type="DEL")

BET_NOA_Catch$Catch <- BET_NOA_Catch$Catch + BET_DEL_Catch$Catch

# prepare to the SS format
BET_PS_Catch <- rbind(BET_OBJ_Catch, BET_NOA_Catch) %>%
  mutate(Seas = 1,
         CV = 0.01,
         Fleet = paste0(Type, "-", Area)) %>%
  select(Year, Seas, Fleet, Catch, CV, Type, Area) %>%
  filter(Year>16)

# add initial catch for each fishery
BET_PS_Catch_int <- BET_PS_Catch %>%
  filter(Year == 100) %>%
  mutate(Year = -999,
         Catch = ifelse(Fleet == "OBJ-A5", 999, 0))

BET_PS_Catch <- rbind(BET_PS_Catch, BET_PS_Catch_int) %>%
  arrange(Fleet, Year)

write.csv(BET_PS_Catch,file=paste0(save_dir,"BET_PS_Catch_1975-",yr.end,".csv"),row.names = FALSE)

ggplot(data = BET_PS_Catch %>% filter(Year > 0)) +
  geom_line(aes(x = Year, y = Catch, color = Type)) +
  facet_wrap( ~ Area, nrow = 5, scales = "free") +
  geom_vline(xintercept = 100.5, linetype = "dashed") +
  theme_bw(16)

BET.OBJ.Comp.20002023 <- read.csv(paste0(save_dir,"BET.OBJ.Comp.20002023.csv"))
BET.OBJ.Comp.19751999 <- read.csv(paste0(save_dir,"BET.OBJ.Comp.19751999.csv"))

BET.NOA.Comp.20002023 <- read.csv(paste0(save_dir,"BET.NOA.Comp.20002023.csv"))
BET.NOA.Comp.19751999 <- read.csv(paste0(save_dir,"BET.NOA.Comp.19751999.csv"))

# BET.DEL.Comp.20002023 <- read.csv(paste0(save_dir,"BET.DEL.Comp.20002023.csv"))
# BET.DEL.Comp.19751999 <- read.csv(paste0(save_dir,"BET.DEL.Comp.19751999.csv"))

BET_OBJ_Comp <- rbind(BET.OBJ.Comp.19751999,BET.OBJ.Comp.20002023) %>%
    mutate(Year=(year-1975)*4+quarter, Type="OBJ") %>%
  arrange(area,Year)
BET_OBJ_Comp <- BET_OBJ_Comp[c(207,206,3:205)]

BET_NOA_Comp <- rbind(BET.NOA.Comp.19751999,BET.NOA.Comp.20002023) %>%
    mutate(Year=(year-1975)*4+quarter, Type="NOA") %>%
  arrange(area,Year)
BET_NOA_Comp <- BET_NOA_Comp[c(207,206,3:205)]

BET_PS_Comp <- rbind(BET_OBJ_Comp,BET_NOA_Comp)

# process PS LF data
names(BET_PS_Comp)[5:205] <- 1:201

Nwells <- BET_PS_Comp %>% select(Type,Year,area,nwells)

BET_PS_Comp_Long <- BET_PS_Comp %>% gather(5:205,key="l",value="lf") %>%
  mutate(L=cut(as.numeric(l), breaks = c(0,seq(22,198,2), Inf), right=F, labels = seq(20,198,2))) %>% na.omit() %>%
  select(Type,Year,area,L,lf) %>%
  group_by(Type,area,Year,L) %>%
  summarise(LF=sum(lf))

BET_PS_Comp_Short <- BET_PS_Comp_Long %>%
  spread(L,LF)

BET_PS_Comp_Final <- left_join(BET_PS_Comp_Short,Nwells) %>%
  filter(Year > 16, nwells > 4)

BET_PS_Comp_SS <- data.frame("Type" = BET_PS_Comp_Final$Type,
                             "Year" = BET_PS_Comp_Final$Year,
                             "Month" = 1,
                             "Fleet" = BET_PS_Comp_Final$area,
                             "Sex" = 0,
                             "Part" = 0,
                             "Nsamp" = BET_PS_Comp_Final$nwells)

BET_PS_Comp_SS <- cbind(BET_PS_Comp_SS, BET_PS_Comp_Final[,4:93], BET_PS_Comp_Final[,4:93])

write.csv(BET_PS_Comp_SS,file=paste0(save_dir,"BET_PS_Comp_1975-",yr.end,".csv"),row.names = FALSE)

names(BET_PS_Comp)[5:205] <- 1:201
BET_PS_Comp_mean <- BET_PS_Comp %>%
  gather(5:205,key="Length",value=comp) %>%
  group_by(Type,area,Length) %>%
  summarise(Comp=sum(comp*nwells)) %>%
  group_by(Type,area) %>%
  mutate(Length=as.numeric(Length),Comp=Comp/sum(Comp))

ggplot(data=BET_PS_Comp_mean) +
  geom_line(aes(x=Length,y=Comp,color=Type)) +
  facet_wrap(~area,nrow = 3) +
  theme_bw(16)