OHI Science | Citation policy
This script prepares the final B/Bmsy data: 1. Calculates the 5 year running average of B/Bmsy data generated using the CMSY method 2. Obtains a B/Bmsy value for each catch record (each FAO/OHI/year/species combination), prioritizing RAM data
None.
B/Bmsy values from the RAM Legacy Stock Assessment data are generated in RAM_data_prep.Rmd
B/Bmsy values from the CMSY method are generated in calculate_bbmsy.Rmd
Mean catch data created in catch_data_prep.Rmd
For the CMSY generated B/Bmsy values we use the five year running mean of the values to smooth the data and to account for model uncertainty.
cmsy <- read.csv('output/cmsy_bbmsy.csv') %>%
filter(!is.na(bbmsy_mean)) %>%
dplyr::select(stock_id, year, bbmsy_q2.5,bbmsy_q97.5,bbmsy_sd, bbmsy_mean, model) %>%
arrange(stock_id, year) %>%
group_by(stock_id) %>%
mutate(mean_5year = rollmean(bbmsy_mean, 5, align="right", fill=NA))
write.csv(cmsy, "int/cmsy_b_bmsy_mean5yrs.csv", row.names=FALSE)A few regions have multiple RAM stocks for the same species (see scientific name). In these cases, we will average the B/Bmsy values of the species, weighted by the area of the RAM stock.
Read in the three data tables:
cmsy <- read.csv('int/cmsy_b_bmsy_mean5yrs.csv') %>%
dplyr::select(stock_id, year, cmsy_bbmsy=mean_5year)
ram <- read.csv("int/ram_bmsy.csv") %>% # final output from RAM_data_prep
rename(stock_id = stockid) # to match other two data tables
mean_catch <- read.csv("output/mean_catch.csv", stringsAsFactors = FALSE) %>% # final output from Watson catch
mutate(taxon_key = str_extract(stock_id_taxonkey, "(\\d)+$")) %>% # extract ending consecutive digits
mutate(stock_id = str_extract(stock_id_taxonkey, "^(\\w+).(\\d){1,2}")) Check number of Watson stocks that have CMSY or RAM B/Bmsy values:
## Watson v RAM BBMSY
setdiff(ram$stock_id, mean_catch$stock_id)
setdiff(mean_catch$stock_id, ram$stock_id)
intersect(ram$stock_id, mean_catch$stock_id) #342 stocks with RAM-B/Bmsy data (although RAM is matched by fao and rgn ids) - v2019
#332 stocks with RAM-B/Bmsy data (although RAM is matched by fao and rgn ids) - v2020
## Watson v CMSY
setdiff(cmsy$stock_id, mean_catch$stock_id)
setdiff(mean_catch$stock_id, cmsy$stock_id)
intersect(mean_catch$stock_id, cmsy$stock_id) #738 stocks with CMSY-B/Bmsy data - v2019
#703 stocks with CMSY-B/Bmsy data - v2020Combine Watson to RAM-B/Bmsy:
data <- mean_catch %>%
left_join(ram, by=c('rgn_id', 'stock_id', "year"))# 608006 catch records (catch from specific fao and ohi regions) when joined with ram increases to 612599 because there are multiple stocks for some species - v2019
# 704173 catch records (catch from specific fao and ohi regions) when joined with ram increases to 709274 because there are multiple stocks for some species - v2020
sum(!is.na(data$ram_bmsy))/nrow(data) # about 6% of catch records have RAM data - v2020
sum(data$mean_catch[!is.na(data$ram_bmsy)])/sum(data$mean_catch) # about 54% of tons of catch have RAM data - v2020 Save & view duplicate stocks:
sum(duplicated(paste(data$rgn_id, data$stock_id, data$year, sep="_")))
# 7773 regions with multiple RAM stocks (stockid_ram) for the same species (see scientific name in stockid)
## save the duplicate stock values to take a look at an example
tmp <- data[duplicated(paste(data$rgn_id, data$stock_id, data$year, sep="_")), ]
## Examples of a region with multiple RAM stocks of the same species
filter(data, rgn_id == 9, year == 2001, stock_id == "Thunnus_alalunga-71") # stocks ALBANPAC and ALBASPAC are the same species; mean catch 1153.528 tonnes
filter(data, rgn_id == 9, year == 2001, stock_id == "Tetrapturus_audax-71") # stocks STMARLINWCNPAC and STMARLINSWPO are the same species; mean catch 86 tonnesRegions with multiple stocks of the same species will have B/Bmsy values averaged, weighted by the area of the RAM stock within the region
## Group by location, year, and species before taking a weighted mean of the catch
data <- data %>%
group_by(rgn_id, taxon_key, stock_id, year, mean_catch) %>%
summarize(ram_bmsy = ifelse(all(!is.na(RAM_area_m2)), weighted.mean(ram_bmsy, RAM_area_m2, na.rm=TRUE), mean(ram_bmsy, na.rm = TRUE)),
gapfilled = ifelse(all(is.na(gapfilled)), NA, max(gapfilled, na.rm=TRUE)),
method = paste(method, collapse = ", ")) %>%
ungroup()
## check that averaging went ok - compare with mean catch values earlier (1153 and 86)
filter(data, rgn_id == 9, year == 2001, stock_id == "Thunnus_alalunga-71")
filter(data, rgn_id == 9, year == 2001, stock_id == "Tetrapturus_audax-71")
## check example of duplicate stock catch with ram_bmsy but no RAM_area_m2 value, ram_bmsy should not be NA (should be 29256 tonnes, checked)
filter(data, rgn_id == 20, year == 2001, stock_id == "Todarodes_pacificus-61")
# add in the B/Bmsy values from the CMSY approach
data <- data %>%
left_join(cmsy, by=c("stock_id", "year"))A lot of the NAs for both RAM-bmsy and CMSY are due to unmatched Watson-RAM stocks
B/Bmsy values for each catch record are generated (for the species where this is possible) and saved. A corresponding gapfilling dataset is also saved.
data_gf <- data %>%
mutate(bmsy_data_source = ifelse(!is.na(ram_bmsy), "RAM", NA)) %>%
mutate(bmsy_data_source = ifelse(is.na(bmsy_data_source) & !is.na(cmsy_bbmsy), "CMSY", bmsy_data_source)) %>%
mutate(bbmsy = ifelse(is.na(ram_bmsy), cmsy_bbmsy, ram_bmsy)) %>%
dplyr::select(rgn_id, stock_id, taxon_key, year, bbmsy, bmsy_data_source, RAM_gapfilled=method, mean_catch) %>%
filter(year >= 2001)
#old <- read_csv(file.path("../v2018/output/fis_bbmsy_gf.csv"))
write.csv(data_gf, "output/fis_bbmsy_gf.csv", row.names=FALSE)
data_gf <- read.csv("output/fis_bbmsy_gf.csv")
bbmsy <- data_gf %>%
dplyr::select(rgn_id, stock_id, year, bbmsy) %>%
dplyr::filter(!is.na(bbmsy))
bbmsy_dups_fixed <- bbmsy %>%
group_by(rgn_id, stock_id, year) %>%
summarise(bbmsy = mean(bbmsy)) %>% ## account for the duplicated TaxonName/CommonName noted in "catch_data_prep.Rmd"...
ungroup()
#old <- read_csv(file.path("../v2019/output/fis_bbmsy.csv"))
write.csv(bbmsy_dups_fixed, "output/fis_bbmsy.csv", row.names=FALSE)