OHI Science | Citation policy
The commercial fishing pressure layers are created from spatialized catch by gear data provided by Watson (2017), and net primary production data from the Vertically Generalized Production Model (VGPM) as described in Behrenfeld and Falkowski (1997)
Three layers are created in this analysis, commercial fishing pressure from high bycatch gear and low bycatch gear as well as artisanal fishing pressure.
Previously, catch data from the Sea Around Us Project was used in conjunction with spatial data from Halpern 2008 to assign proportional catch to each gear type. Watson (2017) published updated sptial catch data by gear type from 1950-2014 and that is used in this assessment as it more accurately represents the type of information needed for these layers.
Reference: Watson (2017)
Downloaded: April 21, 2017
Description: Catch per half degree cell (raw values are in tons per km2)
Native data resolution: 0.5 degree
Time range: 2003 - 2014 (raw data goes back to 1950 but NPP data limits time series to 2003)
Format: Tabular
Load all relevant libraries including parallel processing packages.
#set options for all chunks in code
knitr::opts_chunk$set(warning=FALSE, message=FALSE,fig.width = 6, fig.height = 4, fig.path = 'figs/')
# comment out when knitting
# setwd("globalprep/prs_fish/v2017")
library(tidyverse)
library(parallel)
library(foreach)
library(doParallel)
library(raster)
library(rasterVis)
library(seaaroundus)
library(RColorBrewer)
library(cowplot)
library(stringr)
source('../../../src/R/common.R')
ocean <- raster::raster(file.path(dir_M, 'model/GL-NCEAS-Halpern2008/tmp/ocean.tif'))
registerDoParallel(14) #registering cores for parallel processing
options(scipen=999)
cols = rev(colorRampPalette(brewer.pal(11, 'Spectral'))(255)) # rainbow color scheme
mytheme=rasterTheme(region=cols)
#set mollweide projection CRS
mollCRS=crs('+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs')
#paths
dir_git <- file.path(dir_M,'git-annex/globalprep/prs_fish/v2017')First get the template raster with a resolution of 0.5 degree cells. The getcells() function comes from the seaaroundus R package.
The values of these cells are the Cell ID numbers. In the fisheries catch dataset, these Cell ID numbers match up with the “Seq” numbers.
saup_cells <- getcells("POLYGON ((-180 90, 180 90, 180 -90, -180 -90, -180 90))")
saup_rast <- raster(ncol=720, nrow=360)
saup_rast[] <- saup_cells
plot(saup_rast,col=cols,main = "SAUP Cell IDs")
For each year read in the raw data, filter for appropriate type (high, low and artisanal) and then summarize total catch per cell (in tons/km2)
years = c(2003:2014)
foreach(yr = years) %dopar%{ #yr = 2014
#read in raw data for the year
raw <- readRDS(paste0(file.path(dir_M,'marine_threats/impact_acceleration/stressors/comm_fish/int/catch_data_'),yr,'.rds'))
#high bycatch
high <- raw %>%
dplyr::filter(bycatch == "high") %>%
rowwise() %>%
dplyr::mutate(catch = sum(LSF_CR, IUU_CR, Discards_CR, na.rm=TRUE) * 1000) %>% #there shouldnt be NAs but just in case, converting from tonnes to kg to get rid of low values
dplyr::group_by(Seq) %>%
dplyr::summarise(cell_catch = sum(catch))
#rasterize catch by swapping cell ids with
raster::subs(saup_rast, high, by = 1, which = 2, subsWithNA=TRUE, filename = paste0('int/high_bycatch/annual_catch/high_bc_',yr,'.tif'), overwrite=TRUE)
#low bycatch
low <- raw%>%
filter(bycatch == "low")%>%
rowwise()%>%
mutate(catch = sum(LSF_CR, IUU_CR, Discards_CR, na.rm=TRUE) * 1000)%>% #there shouldnt be NAs but just in case
group_by(Seq)%>%
summarise(cell_catch = sum(catch))
raster::subs(saup_rast, low, by = 1, which = 2, subsWithNA=TRUE, filename = paste0('int/low_bycatch/annual_catch/low_bc_',yr,'.tif'),overwrite=T)
#artisanal bycatch
artisanal <- raw%>%
group_by(Seq)%>%
summarise(cell_catch = sum(SSF_CR * 1000)) #small scale fisheries catch rate
raster::subs(saup_rast, artisanal, by = 1, which = 2, subsWithNA=TRUE, filename = paste0('int/artisanal/annual_catch/art_',yr,'.tif'),overwrite=T)
}I wanted to look at how well our high and low bycatch gear designations align with the reported discard rates per cell. While not perfect, the majority of cells show that catch caught with high bycatch gear leads to a higher proportion of discards than when caught with low bycatch gear.
#read in raw data for the year
raw <- readRDS(paste0(file.path(dir_M,'marine_threats/impact_acceleration/stressors/comm_fish/int/catch_data_2014.rds')))
df <- raw %>%
dplyr::select(Seq, OceanArea, LSF_CR, SSF_CR, IUU_CR, Discards_CR, bycatch) %>%
rowwise() %>%
mutate(catch = sum(SSF_CR, LSF_CR, IUU_CR, Discards_CR)*OceanArea,
discards = Discards_CR * OceanArea) %>%
group_by(Seq, bycatch) %>%
mutate(gear_cell_catch = sum(catch)) %>%
ungroup() %>%
dplyr::select(-LSF_CR, -OceanArea, -SSF_CR, -IUU_CR, -Discards_CR) %>%
group_by(Seq, bycatch) %>%
mutate(disc_cell_catch = sum(discards)) %>%
dplyr::select(-catch, -discards) %>%
distinct() %>%
ungroup() %>%
mutate(disc_prop = disc_cell_catch/gear_cell_catch) %>%
dplyr::select(Seq, bycatch, disc_prop) %>%
spread(bycatch, disc_prop) %>%
mutate(ratio = high/low)
plot(df$high ~ df$low, ylab = "Prop. of catch discarded using high bycatch gears", xlab = "Prop. of catch discarded using low bycatch gears")
abline(0,1,col="red")
The Net Primary Production data was prepared in npp.Rmd.
npp <- list.files(file.path(dir_M,'git-annex/globalprep/prs_fish/v2016/VGPM_primary_productivity/int/annual_npp'), pattern = 'npp_2', full.names=TRUE)
plot(raster(npp[12]), col=cols, axes=FALSE, main = 'Net Primary Production (mg C/m2/day) \n 2014')
Total catch per cell is standardized by the NPP values. This is done because the same fishing pressure can have different impacts depending on the productivity in the region.
Before standardizing, the NPP data is aggregated to the same spatial resolution as the catch data, 0.5 degree cells, to accurately calculate catch in tons/km2 standardized by mg C/m2/day.
npp_stand <- function(file,fname){
yr <- substr(file,nchar(file)-7,nchar(file)-4)
catch <- raster(file)
#get net primary production for that year
n <- npp[substr(npp,111,114)==yr] %>%
raster() %>%
projectRaster(saup_rast) %>%
resample(.,saup_rast)
#3. Divide catch by npp and save to git-annex
overlay(catch,n,fun=function(x,y){x/y},filename = paste0('int/', fname,'/annual_catch_npp/annual_catch_npp_', yr, '.tif'), overwrite=TRUE, progress="text")
}
#apply the function
hb_files <- list.files('int/high_bycatch/annual_catch', full.names=TRUE)
lb_files <- list.files('int/low_bycatch/annual_catch', full.names=TRUE)
art_files <- list.files('int/artisanal/annual_catch', full.names=TRUE)
lapply(hb_files,npp_stand,fname = "high_bycatch")
lapply(lb_files,npp_stand,fname = "low_bycatch")
lapply(art_files,npp_stand,fname = "artisanal")Mean catch per cell is calculated over a rolling window of 5 years to account for interannual variability. The data is then log transformed.
lb_npp <- list.files('int/low_bycatch/annual_catch_npp', full.names=TRUE)
hb_npp <- list.files('int/high_bycatch/annual_catch_npp', full.names=TRUE)
art_npp <- list.files('int/artisanal/annual_catch_npp', full.names=TRUE)
registerDoParallel(6)
foreach (i = 2003:2010) %dopar%{ # i = 2010
yrs <- c(i:(i+4))
out_lb <- lb_npp[which(substr(lb_npp, 51, 54) %in% yrs)] %>%
stack() %>%
calc(fun=function(x){mean(x, na.rm=TRUE)}) %>%
calc(fun=function(x){log(x+1)}, filename = paste0('int/low_bycatch/five_year_means/mean_catch_',yrs[1],'_',yrs[5],'.tif'),overwrite=TRUE)
out_hb <- hb_npp[which(substr(hb_npp, 52, 55) %in% yrs)] %>%
stack() %>%
calc(fun=function(x){mean(x,na.rm=TRUE)}) %>%
calc(fun=function(x){log(x+1)}, filename = paste0('int/high_bycatch/five_year_means/mean_catch_',yrs[1],'_',yrs[5],'.tif'), overwrite=TRUE)
out_art <- art_npp[which(substr(art_npp,49,52) %in% yrs)]%>%
stack()%>%
calc(fun=function(x){mean(x,na.rm=T)})%>%
calc(fun=function(x){log(x+1)}, filename = paste0('int/artisanal/five_year_means/mean_catch_',yrs[1],'_',yrs[5],'.tif'),overwrite=T)
}Look at all mean catch data standardized by NPP and use the 99.99th quantile as the reference point.
mean_hb <- list.files('int/high_bycatch/five_year_means', full.names=TRUE)
mean_lb <- list.files('int/low_bycatch/five_year_means', full.names=TRUE)
mean_art <- list.files('int/artisanal/five_year_means', full.names = TRUE)
# HIGH BYCATCH
#get data across all years
vals <- c()
for(i in 2007:2014){ # i = 2007
#print(i)
m <- mean_hb[which(str_sub(mean_hb, -8, -5) == i)] %>%
raster()%>%
getValues()
n <- m[!is.na(m)]
vals <- c(vals, n)
}
ref_hb <- quantile(vals, prob = 0.9999, na.rm=T) #9.48
#LOW BYCATCH
#get data across all years
lb_vals <- c()
for(i in 2007:2014){
m <- mean_lb[which(str_sub(mean_lb, -8, -5) == i)] %>%
raster() %>%
getValues()
n <- m[!is.na(m)]
lb_vals <- c(lb_vals, n)
}
ref_lb <- quantile(lb_vals, prob = 0.9999, na.rm=TRUE) #10.24
#Artisanal
#LOW BYCATCH
#get data across all years
art_vals <- c()
for(i in 2007:2014){
m <- mean_art[which(str_sub(mean_art, -8, -5) == i)]%>%
raster()%>%
getValues()
n <- m[!is.na(m)]
art_vals <- c(art_vals,n)
}
ref_art <- quantile(art_vals, prob = 0.9999, na.rm=TRUE) #6.89The reference point is 8.6706183 for high bycatch, 9.4098829 for low bycatch and 5.9811736 for artisanal.
mean_hb <- list.files('int/high_bycatch/five_year_means',full.names=T)
mean_lb <- list.files('int/low_bycatch/five_year_means',full.names=T)
mean_art <- list.files('int/artisanal/five_year_means',full.names = T)
registerDoParallel(8)
foreach (i = 2003:2010) %dopar%{ # i = 2003
yrs <- c(i:(i+4))
out_hb <- mean_hb[which(substr(mean_hb, 45, 48) == i)] %>%
raster()%>%
calc(fun=function(x){ifelse(x>ref_hb, 1, x/ref_hb)}) %>%
calc(fun=function(x){ifelse(x<0, 0, x)}) %>%
projectRaster(crs = mollCRS, over=TRUE, method = 'ngb') %>%
resample(ocean, method = 'ngb', filename = paste0(file.path(dir_git), '/output/high_bycatch/hb_fish_pressure_', yrs[1], '-', yrs[5], '.tif'), overwrite = TRUE)
out_lb <- mean_lb[which(substr(mean_lb, 44, 47) == i)] %>%
raster() %>%
calc(fun=function(x){ifelse(x>ref_lb, 1, x/ref_lb)}) %>%
calc(fun=function(x){ifelse(x<0, 0, x)}) %>%
projectRaster(crs = mollCRS, over=TRUE, method = "ngb") %>%
resample(ocean, method = 'ngb', filename = paste0(file.path(dir_git),'/output/low_bycatch/lb_fish_pressure_',yrs[1],'-',yrs[5],'.tif'),overwrite=T)
out_art <- mean_art[which(substr(mean_art, 42, 45) == i)] %>%
raster() %>%
calc(fun = function(x){ifelse(x>ref_art, 1, x/ref_art)}) %>%
calc(fun = function(x){ifelse(x<0, 0, x)}) %>%
projectRaster(crs = mollCRS, over=TRUE, method = "ngb") %>%
resample(ocean, method = 'ngb', filename = paste0(file.path(dir_git), '/output/artisanal/art_fish_pressure_', yrs[1], '-', yrs[5], '.tif'), overwrite=TRUE)
}hb <- raster(file.path(dir_git, 'output/high_bycatch/hb_fish_pressure_2010-2014.tif'))
lb <- raster(file.path(dir_git, 'output/low_bycatch/lb_fish_pressure_2010-2014.tif'))
art <- raster(file.path(dir_git, 'output/artisanal/art_fish_pressure_2010-2014.tif'))
s = stack(hb, lb, art)
plot(s, col=cols, axes=FALSE, box=FALSE)
setwd('globalprep/prs_fish/v2017')
# raster/zonal data
zones <- raster(file.path(dir_M, "git-annex/globalprep/spatial/v2017/regions_eez_with_fao_ant.tif")) # raster data
rgn_data <- read.csv(file.path(dir_M, "git-annex/globalprep/spatial/d2014/data/rgn_mol_raster_1km/regionData.csv")) # data for sp_id's used in raster
### Low bycatch data first
# get raster data:
rasts <- list.files(file.path(dir_git,'output/low_bycatch'), full =TRUE)
pressure_stack <- stack(rasts)
# extract data for each region:
regions_stats <- zonal(pressure_stack, zones, fun="mean", na.rm=TRUE, progress="text")
regions_stats2 <- data.frame(regions_stats)
write.csv(regions_stats2, "int/low_bycatch.csv", row.names=FALSE)
setdiff(regions_stats2$zone, rgn_data$ant_id)
setdiff(rgn_data$ant_id, regions_stats2$zone)
data <- merge(rgn_data, regions_stats2, all.y=TRUE, by.x="rgn_id", by.y="zone") %>%
tidyr::gather("year", "pressure_score", starts_with("lb"))
lb_data <- data %>%
dplyr::mutate(year = stringr::str_sub(year, -4, -1)) %>%
dplyr::mutate(year = as.numeric(year)) %>%
dplyr::filter(rgn_id <= 250) %>%
dplyr::select(rgn_id, year, pressure_score)
write.csv(lb_data, "output/lb.csv", row.names=FALSE)
summary(lb_data)
filter(lb_data, is.na(pressure_score))
rasts <- list.files(file.path(dir_git,'output/high_bycatch'), full = TRUE)
pressure_stack <- stack(rasts)
# extract data for each region:
regions_stats <- zonal(pressure_stack, zones, fun="mean", na.rm=TRUE, progress="text")
regions_stats2 <- data.frame(regions_stats)
write.csv(regions_stats2, "int/high_bycatch.csv", row.names = FALSE)
setdiff(regions_stats2$zone, rgn_data$ant_id)
setdiff(rgn_data$ant_id, regions_stats2$zone)
data <- merge(rgn_data, regions_stats, all.y=TRUE, by.x="rgn_id", by.y="zone") %>%
tidyr::gather("year", "pressure_score", starts_with("hb"))
hb_data <- data %>%
dplyr::mutate(year = stringr::str_sub(year, -4, -1)) %>%
dplyr::mutate(year = as.numeric(year)) %>%
dplyr::filter(rgn_id <= 250) %>%
dplyr::select(rgn_id, year, pressure_score) # na values 71, 72, 74, 75, 188, 215
write.csv(hb_data, "output/hb.csv", row.names=FALSE)
rasts <- list.files(file.path(dir_git,'output/artisanal'), full = TRUE)
pressure_stack <- stack(rasts)
# extract data for each region:
regions_stats <- zonal(pressure_stack, zones, fun="mean", na.rm=TRUE, progress="text")
regions_stats2 <- data.frame(regions_stats)
write.csv(regions_stats2, "int/artisanal.csv", row.names = FALSE)
setdiff(regions_stats2$zone, rgn_data$ant_id) # antarctica regions are in there, makes sense....no land
setdiff(rgn_data$ant_id, regions_stats2$zone) # 213 is in there, that makes sense (Antarctica)
data <- merge(rgn_data, regions_stats, all.y=TRUE, by.x="rgn_id", by.y="zone") %>%
tidyr::gather("year", "pressure_score", starts_with("art"))
art_data <- data %>%
dplyr::mutate(year = stringr::str_sub(year, -4, -1)) %>%
dplyr::mutate(year = as.numeric(year)) %>%
dplyr::filter(rgn_id <= 250) %>%
dplyr::select(rgn_id, year, pressure_score)
write.csv(art_data, "output/art.csv", row.names=FALSE)
## compare with previous year's data
new <- read.csv("output/lb.csv") %>%
filter(year == 2010) %>%
dplyr::select(rgn_id, new_pressure_score = pressure_score)
old <- read.csv("../v2016/output/comm_fish_lb_2016.csv") %>%
left_join(new, by = 'rgn_id')
plot(old$new_pressure_score, old$pressure_score)
abline(0, 1, col="red")
new <- read.csv("output/hb.csv") %>%
filter(year == 2010) %>%
dplyr::select(rgn_id, new_pressure_score = pressure_score)
old <- read.csv("../v2016/output/comm_fish_hb_2016.csv") %>%
left_join(new, by = 'rgn_id')
plot(old$new_pressure_score, old$pressure_score)
abline(0, 1, col="red")
new <- read.csv("output/art.csv") %>%
filter(year == 2010) %>%
dplyr::select(rgn_id, new_pressure_score = pressure_score)
old <- read.csv("../v2016/output/artisanal_fish_lb_2016.csv") %>%
left_join(new, by = 'rgn_id')
plot(old$new_pressure_score, old$pressure_score)
abline(0, 1, col="red")Watson, R. A. (2017). A database of global marine commercial, small-scale, illegal and unreported fisheries catch 1950–2014. Scientific Data, 4.