Skip to content

Commit

Permalink
copied from optimal_allocation_goa repo, modified
Browse files Browse the repository at this point in the history
  • Loading branch information
@zoyafuso-NOAA
zoyafuso-NOAA committed Dec 31, 2023
1 parent 3d6c7d9 commit bf00dc4
Show file tree
Hide file tree
Showing 9 changed files with 22,102 additions and 499 deletions.
6 changes: 4 additions & 2 deletions R/deprecated functions/goa_allocate_stations.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,13 +106,15 @@ goa_allocate_stations <-
FUN = function(x)
sapply(X = split(x = x,
f = optim_df$STRATUM),
FUN = function(xx) sd(as.vector(xx))))[strata_names,]
FUN = function(xx)
stats::sd(x = as.vector(xx))))[strata_names,]
strs_mean <- apply(X = dens[, spp_idx, ],
MARGIN = 2,
FUN = function(x)
sapply(X = split(x = x,
f = optim_df$STRATUM),
FUN = function(xx) mean(as.vector(xx))))[strata_names,]
FUN = function(xx)
stats::mean(x = as.vector(xx))))[strata_names,]

strs_stats <- cbind(data.frame(
STRATO = 1:length(strata_names),
Expand Down
252 changes: 52 additions & 200 deletions analysis_scripts/GOA/1A_operating_model_data.R
View file
Original file line number Diff line number Diff line change
@@ -1,207 +1,59 @@
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Project: GOA Groundfish CPUE Data Synthesis
## Author: Zack Oyafuso ([email protected])
## Contributors: Lewis Barnett ([email protected])
## Description: Create CPUE dataset used for VAST for species of interest
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
###############################################################################
## Project: GOA Groundfish CPUE Data Synthesis
## Author: Zack Oyafuso ([email protected])
## Contributors: Lewis Barnett ([email protected])
## Description: Create CPUE dataset used for VAST for species of interest
###############################################################################
rm(list = ls())

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Import libraries
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
library(dplyr)
library(terra)
library(sumfish)

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Constants used
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## crs used
lonlat_crs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
utm_crs <- "+proj=utm +zone=5N +units=km"

## years used
start_year <- 1996
current_year <- 2021

##################################################
#### Import CPUE survey data
#### Import Haul-level data
#### Import Species codes
#### Import gapindex, Connect to Oracle
##################################################
sumfish::getSQL()
GOA <- sumfish::getRacebase(year = c(start_year, current_year),
survey = 'GOA')

data <- sumfish::sumHaul(GOA) %>% dplyr::mutate(REGION = "GOA")
species_codes <- GOA$species

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Merge bathy rasters ----
## These rasters are really dense, so they are stored in parts and then
## "puzzled" togethered using terra::merge()
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
split_bathy <- list()

n_split_rasters <- length(grep(x = dir("data/GOA/processed_rasters/"),
pattern = "aigoa"))
for (i in 1:n_split_rasters) {
split_bathy[[i]] <- terra::rast(x = paste0("data/GOA/processed_rasters/",
"aigoa_", i, ".tif"))
}

bathy <- do.call(what = terra::merge, args = split_bathy)
rm(split_bathy, i, n_split_rasters)

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Join species names
## Select and rename columns, dropping rows with mising depths
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
species_codes <- dplyr::select(species_codes, -YEAR_ADDED)
data <- dplyr::inner_join(data, species_codes)
data <- data %>% dplyr::select(YEAR, SURVEY = REGION, HAULJOIN = HAULJOIN,
SURFACE_TEMPERATURE, GEAR_TEMPERATURE,
BOTTOM_DEPTH,
EFFORT, WEIGHT,
LATITUDE = START_LATITUDE,
LONGITUDE = START_LONGITUDE,
SPECIES_NAME, COMMON_NAME) %>%
tidyr::drop_na(BOTTOM_DEPTH, LATITUDE, LONGITUDE)

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Filter to GOA survey, remove tows with 0 bottom depth, and drop 2001,
## the year when the survey was incomplete and years before 1996 when a
## different net/soak time was used
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
data <- data %>% filter(SURVEY == "GOA",
BOTTOM_DEPTH > 0,
YEAR != 2001 & YEAR >= 1996)

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Sum catches of northern and southern rock sole with rock sole unid.
## (not distinguished until 1996), rename species complex
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
rock_soles <- data %>%
dplyr::filter(COMMON_NAME %in% c("rock sole unid.",
"southern rock sole",
"northern rock sole")) %>%
dplyr::group_by_at(vars(-WEIGHT, -COMMON_NAME, -SPECIES_NAME)) %>%
dplyr::summarise(WEIGHT = sum(WEIGHT)) %>%
dplyr::ungroup() %>%
dplyr::mutate(SPECIES_NAME = "Lepidopsetta spp.",
COMMON_NAME = "rock soles")

data <- as.data.frame(rbind(data, rock_soles))
library(gapindex)
sql_channel <- gapindex::get_connected()

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Sum catches of blackspooted and rougheye rocks with rougheye and
## blackspotted rockfish unid., rename species complex
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
B_R_rockfishes <- data %>% dplyr::filter(
COMMON_NAME %in% c("blackspotted rockfish",
"rougheye rockfish",
"rougheye and blackspotted rockfish unid.")) %>%
dplyr::group_by_at(vars(-WEIGHT, -COMMON_NAME, -SPECIES_NAME)) %>%
dplyr::summarise(WEIGHT = sum(WEIGHT)) %>%
dplyr::ungroup() %>%
dplyr::mutate(SPECIES_NAME = "Sebastes B_R",
COMMON_NAME = "BS and RE rockfishes")
data <- as.data.frame(rbind(data, B_R_rockfishes))

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Sum catches of sculpins
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
sculpins <- data %>% dplyr::filter(
COMMON_NAME %in% c("bigeye sculpin", "great sculpin",
"plain sculpin", "yellow Irish lord")) %>%
dplyr::group_by_at(vars(-WEIGHT, -COMMON_NAME, -SPECIES_NAME)) %>%
dplyr::summarise(WEIGHT = sum(WEIGHT)) %>%
dplyr::ungroup() %>%
dplyr::mutate(SPECIES_NAME = "sculpins",
COMMON_NAME = "sculpins")
data <- as.data.frame(rbind(data, sculpins))

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Change name of spiny dogfish to Pacific spiny dogifsh
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
data$COMMON_NAME[data$COMMON_NAME == "spiny dogfish"] <-
"Pacific spiny dogfish"

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Subset data to only the species of interest
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
data <- subset(data,
COMMON_NAME %in% c(
## Species included in the survey optimization
"arrowtooth flounder", ## Atherestes stomias
"Pacific cod", ## Gadus macrocephalus
"walleye pollock", ## Gadus chalcogrammus
"rex sole", ## Glyptocephalus zachirus
"flathead sole", ## Hippoglossoides elassodon
"Pacific halibut", ## Hippoglossus stenolepis
"southern rock sole", ## Lepidopsetta bilineata
"northern rock sole", ## Lepidopsetta polyxystra
"Pacific ocean perch", ## Sebastes alutus
"silvergray rockfish", ## Sebastes brevispinis
"northern rockfish", ## Sebastes polyspinis
"dusky rockfish", ## Sebastes variabilis
"BS and RE rockfishes", ## Sebastes aleutianus and S. melanostictus
"Dover sole", ## Microstomus pacificus
"shortspine thornyhead", ## Sebastolobus alascanus

## Species not included in the survey optimization, but
## included when simulating surveys
"sablefish",
"Atka mackerel",
"shortraker rockfish",
"Pacific spiny dogfish",
"yelloweye rockfish",
"giant octopus",
"longnose skate",
"big skate",
"harlequin rockfish",
"giant grenadier",
"sculpins"
))

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Assign station depths from EFH layer
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
cpue_shape <- terra::vect(x = data[, c("LONGITUDE", "LATITUDE")],
geom = c("LONGITUDE", "LATITUDE"),
keepgeom = TRUE,
crs = lonlat_crs)
cpue_shape_aea <- terra::project(x = cpue_shape,
y = terra::crs(bathy))
cpue_shape_aea$depth <- terra::extract(x = bathy,
y = cpue_shape_aea)$AIGOA_ba
cpue_shape_aea[, c("COMMON_NAME", "BOTTOM_DEPTH")] <-
data[, c("COMMON_NAME", "BOTTOM_DEPTH")]

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Plot correlation between EFH depths and reported depth from BTS
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
plot(depth ~ BOTTOM_DEPTH,
data = as.data.frame(cpue_shape_aea)[cpue_shape_aea$COMMON_NAME == "Pacific cod", ],
xlab = "Depth recorded by the BTS",
ylab = "Depth extracted from EFH layer")
cor(cpue_shape_aea$depth, cpue_shape_aea$BOTTOM_DEPTH, use = "complete.obs")
abline(a = 0, b = 1)

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Attach depths to dataset, scaled
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
data$DEPTH_EFH <- cpue_shape_aea$depth
data$DEPTH_EFH[is.na(data$DEPTH_EFH)] <-
data$BOTTOM_DEPTH[is.na(data$DEPTH_EFH)]
data$LOG10_DEPTH_EFH <- log10(data$DEPTH_EFH)
data$LOG10_BOTTOM_DEPTH <- log10(data$BOTTOM_DEPTH)
##################################################
#### Set up haul-level CPUE survey
##################################################
spp_set <- read.csv(file = "data/GOA/species_list/species_list.csv")

## Pull Data from RACEBASE
gapindex_data <- gapindex::get_data(survey_set = "GOA",
year_set = c(1996, 1999,
seq(from = 2003,
to = format(x = Sys.Date(),
format = "%Y"),
by = 2)),
spp_codes = spp_set[,-3],
pull_lengths = FALSE,
sql_channel = sql_channel)

## Calculate and zero-fill CPUE
gapindex_cpue <- gapindex::calc_cpue(racebase_tables = gapindex_data)

## Format catch and effort data for VAST. Note: The `AreaSwept_km2` field set
## to 1 when using CPUE in the `Catch_KG` field.
data_geostat <-
with(gapindex_cpue,
data.frame(Hauljoin = HAULJOIN,
Species = SPECIES_CODE,
Year = YEAR,
Vessel = "missing",
AreaSwept_km2 = 1,
Catch_KG = CPUE_KGKM2,
Lat = LATITUDE_DD_START,
Lon = LONGITUDE_DD_START,
Depth_m = DEPTH_M,
LOG10_DEPTH_M = log10(x = DEPTH_M),
LOG10_DEPTH_M_CEN = scale(x = log10(x = DEPTH_M)),
Pass = 0) )

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Save
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
data <- data[order(data$YEAR, data$SPECIES_NAME),]
write.csv(x = data,
file = "data/GOA/goa_vast_data_input.csv",
##################################################
#### Save
##################################################
if (!dir.exists(paths = "data/GOA/vast_data/"))
dir.create(path = "data/GOA/vast_data/")
write.csv(x = data_geostat,
file = "data/GOA/vast_data/goa_data_geostat.csv",
row.names = F)
saveRDS(object = data_geostat, file = "data/GOA/vast_data/goa_data_geostat.RDS")
Loading

0 comments on commit bf00dc4

Please sign in to comment.