Skip to content

Commit

Permalink
first commit ref #6
Browse files Browse the repository at this point in the history
  • Loading branch information
zoyafuso-NOAA committed Nov 1, 2023
1 parent e42e9e0 commit 180e7d4
Show file tree
Hide file tree
Showing 2 changed files with 218 additions and 0 deletions.
162 changes: 162 additions & 0 deletions code/goa_split_fractions/goa_split_fractions.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,162 @@
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Project: Reproduce the GOA.SPLIT_FRACTIONS tables using the gapindex
## R package
## Author: Zack Oyafuso ([email protected])
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## Restart R Session before running
rm(list = ls())

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Setup
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Import Libraries
library(gapindex)

## Connect to Oracle
sql_channel <- gapindex::get_connected()

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Pull Data from the gapindex R Package
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Pull GOA Management Groups from the GOA schema
goa_mgmt_spp <- RODBC::sqlQuery(channel = sql_channel,
query = "SELECT MANAGEMENT_GROUP, SPECIES_CODE
FROM GOA.GOA_MANAGEMENT_GROUPS")
names(x = goa_mgmt_spp)[1] <- "GROUP"

## Pull data from gapindex based on the species groupings from
## GOA.GOA_MANAGEMENT_GROUPS for only the past few survey years.
gp_data <- gapindex::get_data(year_set = c(2019, 2021, 2023),
survey_set = "GOA",
spp_codes = goa_mgmt_spp,
sql_channel = sql_channel,
pull_lengths = F)

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Post-stratify hauls based on the NMFS Areas
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Calculate midpoint of the observed start and end longitudes for each haul
gp_data$haul$MID_LONGITUDE <-
(gp_data$haul$START_LONGITUDE + gp_data$haul$END_LONGITUDE) / 2

## Following ... reclassify hauls
gp_data$haul$STRATUM[gp_data$haul$STRATUM == 40 &
gp_data$haul$MID_LONGITUDE >= -140] <- 50
gp_data$haul$STRATUM[gp_data$haul$STRATUM == 142 &
gp_data$haul$MID_LONGITUDE < -140] <- 141
gp_data$haul$STRATUM[gp_data$haul$STRATUM == 142 &
gp_data$haul$MID_LONGITUDE >= -140] <- 152

gp_data$haul$STRATUM[gp_data$haul$STRATUM == 143] <- 152 ## Note: new stratum

gp_data$haul$STRATUM[gp_data$haul$STRATUM == 240 &
gp_data$haul$MID_LONGITUDE >= -140] <- 250
gp_data$haul$STRATUM[gp_data$haul$STRATUM == 241 &
gp_data$haul$MID_LONGITUDE >= -140] <- 250

gp_data$haul$STRATUM[gp_data$haul$STRATUM == 340 &
gp_data$haul$MID_LONGITUDE >= -140] <- 350
gp_data$haul$STRATUM[gp_data$haul$STRATUM == 341 &
gp_data$haul$MID_LONGITUDE >= -140] <- 351

gp_data$haul$STRATUM[gp_data$haul$STRATUM == 440 &
gp_data$haul$MID_LONGITUDE >= -140] <- 450
gp_data$haul$STRATUM[gp_data$haul$STRATUM == 540 &
gp_data$haul$MID_LONGITUDE >= -140] <- 550

## This is a hard-copied version of GOA.SPLIT_something
updated_stratum_area <-
data.frame(
SURVEY_DEFINITION_ID = 47, SURVEY = "GOA", DESIGN_YEAR = 1984,
STRATUM = c(40, 41, 50, 140, 141, 152, 150, 151, 240, 241, 250, 251,
340, 341, 350, 351, 440, 450, 540, 550),
AREA_NAME = NA, DESCRIPTION = NA,
AREA_KM2 = c(4980.005501, 6714.745, 11514.5325, 7346.035, 9993.915778,
12043.32322, 4196.599, 6888.172, 2286.139849, 1503.635678,
1882.079151, 4551.089322, 751.2781795, 1296.716466,
2700.424821, 996.6395344, 1252.954196, 1249.897804,
1609.551032, 1484.372968) )

## Update the stratum information in the gp_data list
gp_data$strata <-
rbind(subset(x = gp_data$strata,
subset = !(STRATUM %in% updated_stratum_area$STRATUM)),
updated_stratum_area)

## Update the stratum groupings in the gp_data list to include the new stratum
## 152. Note I'm only modifying the stratum grouping for AREA_ID 959 because
## it is the only AREA_ID that is affected.
gp_data$stratum_groups <- rbind(
gp_data$stratum_groups,
data.frame(SURVEY_DEFINITION_ID = 47,
SURVEY = "GOA",
AREA_ID = c(959),
DESIGN_YEAR = 1984,
STRATUM = 152))

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Calculate CPUE, stratum biomass, then aggregate biomass across subareas
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
gp_cpue <-
gapindex::calc_cpue(racebase_tables = gp_data)
gp_stratum_biomass <-
gapindex::calc_biomass_stratum(racebase_tables = gp_data,
cpue = gp_cpue)
gp_subarea_biomass <-
gapindex::calc_biomass_subarea(racebase_tables = gp_data,
biomass_strata = gp_stratum_biomass)

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Reformat gp_subarea_biomass to something that looks like
## GOA.GOA_SPLIT_FRACTION
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
goa_split_fractions <-
gp_subarea_biomass[gp_subarea_biomass$AREA_ID %in% c(949, 959),
c("YEAR", "SPECIES_CODE", "AREA_ID", "BIOMASS_MT")]
goa_split_fractions <-
reshape(data = goa_split_fractions,
idvar = c("SPECIES_CODE", "YEAR"),
timevar = "AREA_ID",
direction = "wide")
names(x = goa_split_fractions) <- c("YEAR", "MANAGEMENT_GROUP",
"WEST_BIOMASS", "EAST_BIOMASS")
goa_split_fractions$WEST_BIOMASS <-
round(x = goa_split_fractions$WEST_BIOMASS, digits = 1)
goa_split_fractions$EAST_BIOMASS <-
round(x = goa_split_fractions$EAST_BIOMASS, digits = 1)

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Compare output from gapindex to what is currently in the GOA schema
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## Pull current split fractions table
historical_split <-
RODBC::sqlQuery(channel = sql_channel,
query = "SELECT YEAR, MANAGEMENT_GROUP,
WEST_BIOMASS, EAST_BIOMASS
FROM GOA.SPLIT_FRACTIONS
WHERE YEAR >= 2019")

## Merge the current split fraction table with the gapindex-produced table
test_split <-
merge(x = historical_split,
y = goa_split_fractions,
by = c("YEAR", "MANAGEMENT_GROUP"),
all.x = TRUE,
suffixes = c("_HIST", "_GP"))

## Calculate absolute differences between the west and east biomass estimates
test_split$WEST_BIOMASS_DIFF <-
test_split$WEST_BIOMASS_GP - test_split$WEST_BIOMASS_HIST
test_split$EAST_BIOMASS_DIFF <-
test_split$EAST_BIOMASS_GP - test_split$EAST_BIOMASS_HIST

write.csv(x = test_split[order(test_split$MANAGEMENT_GROUP),
c("YEAR", "MANAGEMENT_GROUP",
"WEST_BIOMASS_HIST", "WEST_BIOMASS_GP",
"WEST_BIOMASS_DIFF", "EAST_BIOMASS_HIST",
"EAST_BIOMASS_GP", "EAST_BIOMASS_DIFF")],
file = "code/goa_split_fractions/goa_split_fractions.csv",
row.names = F)

56 changes: 56 additions & 0 deletions code/goa_split_fractions/goa_split_fractions.csv
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
"YEAR","MANAGEMENT_GROUP","WEST_BIOMASS_HIST","WEST_BIOMASS_GP","WEST_BIOMASS_DIFF","EAST_BIOMASS_HIST","EAST_BIOMASS_GP","EAST_BIOMASS_DIFF"
2019,"arrowtooth flounder",68988.4,68988.2,-0.19999999999709,142295.9,142296,0.100000000005821
2021,"arrowtooth flounder",53981.6,53981.8,0.200000000004366,95332.3,95332.2,-0.100000000005821
2023,"arrowtooth flounder",92762.8,92762.8,0,191769.6,191769.6,0
2019,"deep water flatfish",16948.8,16948.8,0,13929.1,13929.1,0
2021,"deep water flatfish",11099.8,11100.8,1,16151.1,16150.9,-0.200000000000728
2023,"deep water flatfish",13823.5,13823.6,0.100000000000364,17234.1,17234.2,0.100000000002183
2019,"dusky rockfish",5776.1,5776,-0.100000000000364,2227.8,2227.8,0
2021,"dusky rockfish",659.5,659.6,0.100000000000023,1016.3,1016.3,0
2023,"dusky rockfish",762.6,762.6,0,2335.9,2336,0.0999999999999091
2019,"flathead sole",14118,14117.8,-0.200000000000728,11139.7,11139.7,0
2021,"flathead sole",12631.5,12631.4,-0.100000000000364,15680.2,15680.2,0
2023,"flathead sole",15742.1,15742.1,0,8314.4,8314.3,-0.100000000000364
2019,"harlequin rockfish",342.7,342.7,0,209.5,209.5,0
2021,"harlequin rockfish",50.2,50.2,0,72.7,72.6,-0.100000000000009
2023,"harlequin rockfish",0,0,0,81.7,81.8,0.0999999999999943
2019,"OROXM006",70.6,70.6,0,464.7,464.7,0
2021,"OROXM006",0,0,0,2257.1,2257.2,0.0999999999999091
2023,"OROXM006",1,1,0,4740.1,4740,-0.100000000000364
2019,"OROXM007",547.3,547.3,0,2964.2,2964.2,0
2021,"OROXM007",392.9,392.9,0,1381.4,1381.6,0.199999999999818
2023,"OROXM007",603.9,603.9,0,2066.2,2066.2,0
2019,"Pacific cod",2978,2978,0,6688.1,6688.1,0
2021,"Pacific cod",8837.3,8837.3,0,15218.2,15218,-0.200000000000728
2023,"Pacific cod",13318.6,13318.7,0.100000000000364,12710.1,12710.1,0
2019,"Pacific ocean perch",47901.8,47901.8,0,200557.5,200557.3,-0.200000000011642
2021,"Pacific ocean perch",30463.6,30463.5,-0.0999999999985448,127524.4,127524.4,0
2023,"Pacific ocean perch",39406,39406,0,834319.7,834319.8,0.100000000093132
2019,"redbanded rockfish",1077,1077,0,2587.6,2587.7,0.0999999999999091
2021,"redbanded rockfish",283.2,283.2,0,4818.2,4818,-0.199999999999818
2023,"redbanded rockfish",654.5,654.5,0,4345.1,4344.8,-0.300000000000182
2019,"redstripe rockfish",1631.9,1631.8,-0.100000000000136,10058,10058,0
2021,"redstripe rockfish",76.7,76.7,0,2346.5,2346.5,0
2023,"redstripe rockfish",61.3,61.3,0,6671.9,6672,0.100000000000364
2019,"rex sole",8308.7,8308.7,0,15737.9,15737.9,0
2021,"rex sole",8703.5,8703.3,-0.200000000000728,17086,17086,0
2023,"rex sole",6539.4,6539.4,0,10917.4,10917.3,-0.100000000000364
2023,"rougheye blackspotted complex",4128.6,4128.8,0.199999999999818,17018.8,17019.6,0.799999999999272
2019,"shallow water flatfish",5167.8,5167.7,-0.100000000000364,11641.5,11641.4,-0.100000000000364
2021,"shallow water flatfish",12915.9,12916.1,0.200000000000728,12347.4,12347.3,-0.100000000000364
2023,"shallow water flatfish",16311,16311,0,17227.4,17227.6,0.19999999999709
2019,"sharpchin rockfish",2870.8,2870.7,-0.100000000000364,5669.8,5669.7,-0.100000000000364
2021,"sharpchin rockfish",1394,1394,0,6332.1,6332.2,0.0999999999994543
2023,"sharpchin rockfish",306.8,306.9,0.0999999999999659,5319.2,5319.2,0
2023,"shortraker rockfish",8812.1,8812.1,0,13239.4,13239.4,0
2019,"shortraker/rougheye",19044.1,NA,NA,22088.9,NA,NA
2021,"shortraker/rougheye",8459,NA,NA,8885.1,NA,NA
2019,"shortspine thornyhead",13414,13413.9,-0.100000000000364,10905.3,10905.3,0
2021,"shortspine thornyhead",11383.3,11383.5,0.200000000000728,12355.6,12355.6,0
2023,"shortspine thornyhead",11976.3,11976.2,-0.0999999999985448,12641.2,12641.1,-0.100000000000364
2019,"silvergray rockfish",5181.5,5181.6,0.100000000000364,23577,23577,0
2021,"silvergray rockfish",1480.4,1480.4,0,42528.1,42528,-0.0999999999985448
2023,"silvergray rockfish",12039.1,12039.1,0,30745.8,30745.9,0.100000000002183
2019,"walleye pollock",10952.4,10952.5,0.100000000000364,49695.1,49695.1,0
2021,"walleye pollock",19466.2,19466.3,0.0999999999985448,34273.5,34273.5,0
2023,"walleye pollock",19655.4,19655.5,0.0999999999985448,30767,30767,0

0 comments on commit 180e7d4

Please sign in to comment.