OHI Science | Citation policy
The habitat pressure data layer is created from fish catch by gear type of industrial and non industrial fishing provided by Watson (2019). Total tonnes of catch are calculated for each OHI region and then standardized by area (km^2) soft-bottom habitat in each region. These values are then rescaled to be between 0 and 1.
Instead of using the maximum density, used the 95th quantile (line ~255).
Reference: Watson, R. A. and Tidd, A. 2018. Mapping nearly a century and a half of global marine fishing: 1869–2015. Marine Policy, 93, pp. 171-177. (Paper URL)
Downloaded: December 11, 2019 from IMAS portal
Description: Global fisheries landings data per cell separated by Industrial versus Non-Industrial catch, IUU, and discards.
Native data resolution:
Time range: 1950 - 2017
Format: csv format
First create one raster a year with all catches (industrial and non industrial) using gear that harm soft-bottom habitat. Then extract catches per OHI regions and standardize by soft-bottom habitat of each region. We used a ln(x + 1) transformation because density data were extremely skewed. Finally, we rescale the transformed data dividing by the maximum value across all years. A second rescale was done using the media values across years.
This script provides pressure scores, health indicators, and trend values for soft-bottom habitats.
library(tidyverse)
library(parallel)
library(foreach)
library(doParallel)
library(raster)
library(rasterVis)
library(seaaroundus) # devtools::install_github("ropensci/seaaroundus") if install.packages("seaaroundus") doesn't work.
library(RColorBrewer)
library(cowplot)
library(colorspace)
library(sp)
library(rgdal)
library(sf)
library(readxl)
library(here)
registerDoParallel(4) # Registering cores for parallel processing
source('http://ohi-science.org/ohiprep_v2020/workflow/R/common.R')
## Color palette
cols = rev(colorRampPalette(brewer.pal(11, 'Spectral'))(255)) # Rainbow color scheme
mytheme = rasterTheme(region = cols)
options(scipen = 999)
IMAS_d2020 <- file.path(dir_M, "git-annex/globalprep/_raw_data/IMAS_GlobalFisheriesLandings/d2020")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.
## Create a raster of Cell numbers
# This Codes.xlsx was downloaded from the same place as the raw Watson data.
cells <- read_excel(file.path(IMAS_d2020, "Codes.xlsx"), sheet = "Cell") %>%
dplyr::rename(x = LonCentre, y = LatCentre, z = Cell) %>% # renamed these xyz to use in the rasterFromXYZ() function below
dplyr::select(x,y,z)
## Turn the lat/long points into a raster
cells_raster <- rasterFromXYZ(cells)
crs(cells_raster) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
crs(cells_raster)
# Check out the cells raster to make sure it looks good
plot(cells_raster)Then create one raster per year for industrial and non industrial catches using gear that harm soft-bottom habitat.
years <- c(2003:2017)
## Data frame with each gear in the dataset with the corresponding multiplier according to how much they harm the ocean bottom.
multipliers <- data.frame(FleetGearName = c("Dredge", "Gillnet", "Lines Non-tuna", "Longline Non-tuna", "LonglineTuna", "Other", "Pole and Line Tuna", "Purse seine Non-tuna", "Purse seine Tuna", "Seine", "Trap", "Trawl", "Trawl midwater"), multiplier = c(1, 0.5, 0, 0, 0, 0.25, 0, 0, 0, 0, 0.5, 1, 0))
## Loop to create rasters of catches using gear that harm soft-bottom habitat.
foreach(yr = years) %dopar%{
# yr = 2015 if you want to just run one year without looping through all of them
## Read in raw data for industrial catch per the year
raw_ind <- readRDS(paste0(file.path(dir_M,'git-annex/globalprep/prs_fish/v2020/int/annual_catch/IndCatch_'),yr,'.rds'))
## Wrangle data, discount catch by multiplier and add catches per cell
catch_ind <- raw_ind %>%
dplyr::rowwise() %>%
dplyr::mutate(catch = sum(Reported, IUU, Discards, na.rm=TRUE)) %>%
dplyr::left_join(multipliers, by = "FleetGearName") %>%
dplyr::mutate(catch_discount = catch*multiplier) %>%
dplyr::group_by(Cell) %>%
dplyr::summarise(cell_catch_ind = sum(catch_discount))
## Read in raw data for non-industrial catch per the year
raw_Nind <- readRDS(paste0(file.path(dir_M,'git-annex/globalprep/prs_fish/v2020/int/annual_catch/NIndCatch_'),yr,'.rds'))
## Wrangle data, discount catch by multiplier and add catches per cell
catch_Nind <- raw_Nind %>%
rowwise() %>%
dplyr::mutate(catch = sum(Reported, IUU, Discards, na.rm=TRUE)) %>%
left_join(multipliers, by = "FleetGearName") %>%
dplyr::mutate(catch_discount = catch*multiplier) %>%
dplyr::group_by(Cell) %>%
dplyr::summarise(cell_catch_Nind = sum(catch_discount))
## Sum industrial and non-industrial catches
catch_total <- catch_ind %>%
dplyr::full_join(catch_Nind, by= "Cell") %>%
dplyr::rowwise() %>%
dplyr::mutate(total = sum(cell_catch_ind, cell_catch_Nind, na.rm=TRUE)) %>%
dplyr::select(Cell, total)
## Rasterize catch by swapping cell IDs with (calls the cells raster we made above)
raster::subs(cells_raster, catch_total, by = 'Cell', which = 'total', subsWithNA=TRUE, filename = file.path(dir_M, paste0('git-annex/globalprep/hab_prs_hd_subtidal_soft_bottom/v2020/int/catch_count_',yr,'.tif')), overwrite=TRUE) # Saves raster directly to Mazu
}
## Reading in the raster to make sure it looks as expected
test <- raster(file.path(dir_M, ('git-annex/globalprep/hab_prs_hd_subtidal_soft_bottom/v2020/int/catch_count_2014.tif')))
plot(test) ## Looks like all information is there!## Read tif as rasters for all catches using gear that harm soft-bottom habitat
## Reference raster
catch_crs <- raster(file.path (dir_M,'git-annex/globalprep/hab_prs_hd_subtidal_soft_bottom/v2020/int/catch_count_2003.tif'))
## OHI region polygons
ohi_regions <- st_read(dsn = file.path(dir_M, "git-annex/globalprep/spatial/v2017"), layer="regions_2017_update")
## Make the region coordinate reference system the same as the catch data
ohi_regions_wgs <- st_transform(ohi_regions, proj4string(catch_crs))
## List of all tifs with catches that need to be read in the loop
catch_rasts <- list.files(file.path(dir_M, "git-annex/globalprep/hab_prs_hd_subtidal_soft_bottom/v2020/int"), pattern = 'catch_count_.*.tif', full.names = TRUE)
## Registering cores for parallel processing
registerDoParallel(4)
## Loop to extract values for each OHI region
foreach(sb_catch = catch_rasts) %dopar% {
# sb_catch = catch_rasts[12] if you just want to run one year's raster without looping through all of them
cat(sb_catch)
catch_year <- str_extract(basename(sb_catch), "(\\d)+") # Extract the corresponding year for each file
catch <- raster(sb_catch)
data <- raster::extract(catch, ohi_regions_wgs, weights = TRUE, normalizeWeights = FALSE, progress = 'text')
names(data) <- paste(ohi_regions$type_w_ant, ohi_regions$rgn_ant_id, sep="_")
sb_catch_rgn <- plyr::ldply(data, rbind)
# The following code keeps the raster information when the raster lands partially on the land polygons
sb_catch_rgn <- sb_catch_rgn %>%
tidyr::separate(.id, c("rgn_type", "rgn_id"), sep="_") %>%
dplyr::mutate(tonnes = value*weight) %>%
dplyr::group_by(rgn_type, rgn_id) %>%
dplyr::summarize(tonnes = sum(tonnes, na.rm=TRUE)) %>%
dplyr::group_by(rgn_id) %>%
dplyr::summarize(tonnes = sum(tonnes, na.rm=TRUE)) %>%
dplyr::ungroup() %>%
dplyr::arrange(as.numeric(as.character(rgn_id))) %>%
dplyr::mutate(year = catch_year)
write.csv(sb_catch_rgn, file = file.path(here(), "globalprep/hab_prs_hd_subtidal_soft_bottom/v2020/int", sprintf("sb_catch_rgn_%s.csv", catch_year)), row.names=FALSE)
}area <- read.csv(file.path(here(), "globalprep/hab_prs_hd_subtidal_soft_bottom/v2016/output/habitat_extent_softbottom.csv"))
## Read in csv for each year and combine them in one data frame
files <- list.files("int", pattern = "sb_catch_rgn_.*.csv", full.names = TRUE)
sb_catch_list <- lapply(files, function(x){read.csv(file=x)})
sb_catch_all <- bind_rows(sb_catch_list)
# sum(sb_catch_all$tonnes) # v2018: 707463134
# v2019: 5899354
# v2020: 682566393
data_density <- sb_catch_all %>%
left_join(area, by="rgn_id") %>%
filter(rgn_id <= 250) %>%
mutate(density = tonnes/km2)
## Find 99th quantile density for each year
quant_99 <- data_density %>%
group_by(year) %>%
summarize(quantDensity = quantile(density, c(0.99), na.rm=TRUE)) %>%
data.frame()
write.csv(quant_99, "int/reference_point_quant_99.csv", row.names=FALSE) # Make sure you are working in the correct working directory
# # Scatterplot of 2014 and 2015
# catch_2014 <- read.csv("int/sb_catch_rgn_2014.csv") %>%
# mutate(tonnes_2014 = tonnes) %>%
# dplyr::select(-year)
#
# catch_2015 <- read.csv("int/sb_catch_rgn_2015.csv") %>%
# mutate(tonnes_2015 = tonnes) %>%
# dplyr::select(-year)
#
# catch_combined <- catch_2014 %>%
# left_join(catch_2015, "rgn_id")
#
# catch_plot <- catch_combined %>%
# ggplot(aes(x=tonnes_2014, y=tonnes_2015)) +
# geom_point()
#
# catch_plot
## Reference point = 99th quantile across all years
ref_point <- read.csv("int/reference_point_quant_99.csv")
ref_point_median <- median(ref_point$quantDensity) # v2020 95th: 20.5649376263763
# v2020 99th: 62.23165
# ref_point (max) 2018 = 7052.172 (year 2015)
# ref_point (max) 2019 = 45.76575 (year 2015)
# ref_point (max) 2020 = 4842.57385 (year 2015)
# ref_point (95th quantile) 2020 = 19.62075 (year 2015)
## Checking for distribution of data. In the past density data happens to be very skewed
hist(data_density$density) # Yes, very skewed data
## Rescale the density:
data_rescale <- data_density %>%
mutate(density_rescaled_median = density/ref_point_median) %>%
mutate(density_rescaled_median_capped = ifelse(density_rescaled_median > 1,
1,
density_rescaled_median))
hist(data_rescale$density_rescaled_median_capped)
## code from v2019:
## Rescale the density:
## Density is rescaled twice to reduce skew
# data_rescaled <- data_density %>%
# mutate(density_rescaled_max = log(density + 1)/log(ref_point_max + 1)) %>% # pressure-type measure
# mutate(inv_dens_rescaled_max = 1 - density_rescaled_max) # health-type measure
#
# hist(data_rescaled$density_rescaled_max)
# hist(data_rescaled$inv_dens_rescaled_max)
#
#
# ## Find the second rescaling point - median value of healthy-type measure for each year
# ref_median <- data_rescaled %>%
# group_by(year) %>%
# summarize(ref_median = median(inv_dens_rescaled_max, na.rm=TRUE)) %>%
# data.frame()
#
# write.csv(ref_median, "int/reference_point_median.csv", row.names=FALSE)
#
#
# ref_point_median <- min(ref_median$ref_median) # 0.8590695 (v2018)
# # 0.99485138963124 (v2019)
# # 0.87582237227035 (v2020)
#
#
# ## Rescale for second time:
# data_rescaled_2 <- data_rescaled %>%
# mutate(density_rescaled_median = inv_dens_rescaled_max/ref_point_median) %>%
# mutate(density_rescaled_median_capped = ifelse(density_rescaled_median > 1,
# 1,
# density_rescaled_median))
#
# hist(data_rescaled_2$density_rescaled_median_capped)
#
#
# hist(data_rescaled_2$density_rescaled_median_capped[data_rescaled_2$year==2014])
# filter(data_rescaled_2, rgn_id==163)
# filter(data_rescaled_2, year==2010)## Check against old data
old <- read.csv("../v2018/output/habitat_health_softbottom.csv") %>%
dplyr::filter(year == 2014) %>%
dplyr::select(rgn_id, old_health=health)
compare <- data_rescale %>%
filter(year==2014) %>%
dplyr::select(rgn_id, health = density_rescaled_median_capped) %>%
dplyr::mutate(health = 1 - health) %>%
left_join(old)
plot(health ~old_health, data=compare)
abline(0,1, col = 'red')
# Not a great correlation, but different years data and fisheries data is very different## Get relevant data
save_dir <- "output"
condition_pressure <- data_rescale %>%
mutate(habitat = "soft_bottom") %>%
mutate(pressure = density_rescaled_median_capped) %>%
mutate(health = 1 - density_rescaled_median_capped) %>%
dplyr::select(rgn_id, year, habitat, health, pressure) %>%
filter(!is.na(pressure))
hist(condition_pressure$pressure)
hist(condition_pressure$health)
## Get habitat trends
stop_year <- max(condition_pressure$year)
trend <- data.frame()
for (status_year in (stop_year-4):stop_year){ # status_year = 2011
trend_years <- status_year:(status_year - 4)
first_trend_year <- min(trend_years)
trend_new <- condition_pressure %>%
filter(year %in% trend_years) %>%
group_by(rgn_id) %>%
do(mdl = lm(health ~ year, data=.),
adjust_trend = .$health[.$year == first_trend_year]) %>%
summarize(rgn_id = rgn_id,
trend = round(coef(mdl)['year']/adjust_trend * 5, 4)) %>%
ungroup() %>%
mutate(trend = ifelse(trend > 1, 1, trend)) %>%
mutate(trend = ifelse(trend < (-1), (-1), trend))
trend_new <- trend_new %>%
dplyr::mutate(habitat = "soft_bottom") %>%
dplyr::mutate(year = status_year) %>%
dplyr::select(rgn_id, year, habitat, trend)
trend <- rbind(trend, trend_new)
}
write.csv(trend, file.path(save_dir, "habitat_trend_softbottom.csv"), row.names=FALSE)
health <- condition_pressure %>%
dplyr::mutate(habitat = "soft_bottom") %>%
dplyr::select(rgn_id, year, habitat, health)
write.csv(health, file.path(save_dir, "habitat_health_softbottom.csv"), row.names=FALSE)
pressure <- condition_pressure %>%
dplyr::select(rgn_id, year, pressure_score = pressure)
write.csv(pressure, file.path(save_dir, "hd_sb_subtidal.csv"), row.names=FALSE)Compare trend, pressure, pressure output to last year’s values
## Health data
old_health <- read.csv("../v2018/output/habitat_health_softbottom.csv") %>%
dplyr::filter(year == 2014) %>%
dplyr::select(rgn_id, old_health=health)
compare_health <- read.csv("output/habitat_health_softbottom.csv")%>%
filter(year==2017) %>%
dplyr::select(rgn_id, health) %>%
left_join(old_health)
plot(health ~old_health, data=compare_health)
abline(0,1, col = 'red')
hist(compare_health$health)
## Pressure
old_pressure <- read.csv("../v2018/output/hd_sb_subtidal.csv") %>%
dplyr::filter(year == 2014) %>%
dplyr::select(rgn_id, old_pressure= pressure_score)
compare_pressure <- read.csv("output/hd_sb_subtidal.csv")%>%
filter(year==2014) %>%
dplyr::select(rgn_id, pressure=pressure_score) %>%
left_join(old_pressure)
plot(pressure ~old_pressure, data=compare_pressure)
abline(0,1, col = 'red')
compare_pressure_plot <- ggplot(data = compare_pressure, aes(x = old_pressure, y = pressure, label = rgn_id))+
geom_point()+
geom_abline(color = "red")
plot(compare_plot)
library(plotly)
ggplotly(compare_pressure_plot)
## Trend
old_trend <- read.csv("../v2018/output/habitat_trend_softbottom.csv") %>%
dplyr::filter(year == 2014) %>%
dplyr::select(rgn_id, old_trend= trend)
compare_trend <- read.csv("output/habitat_trend_softbottom.csv")%>%
filter(year==2014) %>%
dplyr::select(rgn_id, trend) %>%
left_join(old_trend)
plot(trend ~old_trend, data=compare_trend)
abline(0,1, col = 'red')
## Comparisons in general are not great but not too bad. Consider: compared to v2018 we included more destructive gear, we also included artisanal (non-commercial catches) and Watson updated its methods.There was no gapfilling for these data. Created gapfill files with values of 0. Note: Every data layer has to have a corresponding gapfill file. If there is no gapfilling the gapfill column is filled in with 0 as done below.
hab_health_sb_gf <- read.csv("output/habitat_health_softbottom.csv")%>%
mutate(gapfilled = 0) %>%
dplyr::select(rgn_id, year, gapfilled)
write.csv(hab_health_sb_gf, "output/habitat_health_softbottom_gf.csv", row.names=FALSE)
hab_trend_sb_gf <- read.csv("output/habitat_trend_softbottom.csv")%>%
mutate(gapfilled = 0) %>%
dplyr::select(rgn_id, year, gapfilled)
write.csv(hab_trend_sb_gf, "output/habitat_trend_softbottom_gf.csv", row.names=FALSE)
hd_sb_subtidal_gf <- read.csv("output/hd_sb_subtidal.csv")%>%
mutate(gapfilled = 0) %>%
dplyr::select(rgn_id, year, gapfilled)
write.csv(hd_sb_subtidal_gf, "output/hd_sb_subtidal_gf.csv", row.names=FALSE)