OHI Science | Citation policy
This script takes the Sea Around Us Project (SAUP) catch data, provided at a resolution of half-degree cells globally, and aggregates catch to stock levels. For the Ocean Health Index, we assume a stock is represented by the FAO region in which the species is caught.
In order to aggregate to FAO regions, we associate each cell to the FAO region and the OHI region in which it is located.
An example of our aggregation proces: New Zealand is located entirely in FAO region 81. All catch reported by New Zealand will be aggregated by species to the FAO region. If a species was reported as caught in both New Zealand waters and in the High Seas of area 81, these two records will be combined into one by summing the catch.
Catch data now goes from 1950 - 2014. Previously catch data was only available through 2010.
The Sea Around Us Project shared the spatialized catch data with OHI on joined to a lookup table that links SAUP region names and ids to the FAO region they are located in. The proportional area of each EEZ within the FAO region was also calculated for overlapping EEZs.
Reference:
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',message = FALSE, warning = FALSE)
## Libraries
library(readr)
library(dplyr)
library(raster)
library(parallel)
library(purrr)
source('~/github/ohiprep/src/R/common.R')
## Paths for data
path_data = file.path(dir_M, "git-annex/globalprep/_raw_data/SAUP/d2017/annual_data")These files are large so using the data.table package is recommended due to R memory limitations.
## information about half degree cells
cells <- read_csv('cells.csv') %>%
select(-X1)Aggregate catch per OHI region and FAO area. This catch will be used twice.
The catch is used to weight scores per region. For this we need to use all catch records, including those not reported at the species level.
The catch data at species level is used to calculate stock status (BBmsy) per stock (remember that our definition of a stock is a species caught within a single FAO area).
## list all data files
data_files <- list.files(file.path(path_data), full.names = T)
## function to wrangle data into what we need (total catch per OHI region per stock)
stock_rgn_catch <- function(df) {
output_df <- readRDS(df) %>%
select(year, taxon_scientific_name, taxon_common_name, cell_id, catch_sum, taxon_key) %>%
rename(CellID = cell_id) %>% #change cellid name to match what is in cells.csv
left_join(cells) %>%
mutate(catch_prop = catch_sum * area) %>%
group_by(year, rgn_id,fao_id, taxon_scientific_name, taxon_common_name, taxon_key) %>%
summarise(catch = sum(catch_prop)) %>%
ungroup() %>%
mutate(stock_id = gsub(" ", "_", paste(taxon_scientific_name, fao_id, sep='-'), fixed=TRUE))%>%
rename(fao_rgn = fao_id,
tons = catch)
return(output_df)
}
## use purrr::map to apply the function to all yearly datasets
a <- map_df(data_files, stock_rgn_catch)
#this is a large file (150 MB) so it is saved on the server
write.csv(a, file = file.path(dir_M,'git-annex/globalprep/fis/v2017/int/stock_catch_by_rgn.csv'),row.names=FALSE)Catch-MSY is the model we use to estimate stock status for all global stocks. This model requires information about the resilience of each species in addition to the catch data.
Load taxonomic resilience information
# add the taxon_resilence data to catch for b/bmsy calculations
taxon_res = read_csv('../v2017/data/taxon_resilience_lookup.csv') %>%
#mutate(common = ifelse(common %in% "Silver croaker", paste(common, sciname, sep=" "), common)) %>%
dplyr::select(taxon_common_name=common, Resilience)Filter out all stocks that don’t meet our conditions:
#set variables to filter by
min_yrs = 20
min_tons = 1000
#read in catch data created above
df <- read_csv(file.path(dir_M,'git-annex/globalprep/fis/v2017/int/stock_catch_by_rgn.csv'))
#create dataset ready to run through catch only models
stks <- df%>%
filter(taxon_key >= 600000, #remove all records of catch reported at higher taxonomic levels than species
tons > 0)%>% #remove records of 0 catch
select(-rgn_id)%>% #remove rgn_id since we aggregate stocks to the FAO level
group_by(stock_id,year,fao_rgn,taxon_scientific_name,taxon_common_name,taxon_key)%>%
summarise(tons = sum(tons))%>% #calculate total tons per stock
ungroup()%>%
group_by(stock_id)%>%
mutate(nyrs = n(), #get the total number of years the stock has records for
avg_ann_catch = mean(tons))%>% #calculate the mean catch over all catch years
ungroup()%>%
filter(avg_ann_catch >= min_tons, #keep only those stocks that meet our conditions
nyrs >= min_yrs)%>%
left_join(taxon_res)%>% #add resilience information
dplyr::select(year,taxon_scientific_name,taxon_common_name,fao_rgn,stock_id,taxon_key,Resilience,tons)
# check on stocks that don't have a resilience
no_res <- filter(stks, is.na(Resilience)) %>%
select(taxon_scientific_name, taxon_common_name) %>%
distinct()
nrow(no_res) #148 species do not have a Resilience. These will get assigned a Medium Resilience by default by the CMSY model.
write.csv(stks, file = 'data/stock_catch.csv')stks = read.csv('data/stock_catch.csv')
DT::datatable(head(stks,n=100))Pauly D. and Zeller D. (Editors), 2015. Sea Around Us Concepts, Design and Data (seaaroundus.org)