OHI Science | Citation policy
[REFERENCE RMD FILE: http://ohi-science.org/ohiprep_v2023/globalprep/prs_coral_harvest/v2023/prs_coral_harvest.html]
This analysis converts FAO commodities data into data layers used to calculate OHI 2023 global coral harvest pressure.
Reference: https://www.fao.org/fishery/statistics-query/en/trade App release date: July 2023 FAO raw commodities quantity 1976_2021 FAO raw commodities value 1976_2021 FAO raw commodities metadata found here and overall FAO data directory found here
Downloaded: July 19, 2023
Description: Quantity (tonnes) and value (USD) of raw commodities (Exports only) for each country, taxa, year. The FAO data is subset to include commodities in these categories: shells, corals, ornamental fish, fish oil, seaweed and plants, sponges (see: raw/commodities2products.csv for details).
Time range: 1976-2021
knitr::opts_chunk$set(eval=FALSE)
## load libraries, set directories
library(ohicore) #devtools::install_github('ohi-science/ohicore@dev')
library(dplyr)
library(stringr)
library(tidyr)
library(zoo)
library(ggplot2)
library(here)
library(tidyverse)
library(plotly)
library(tictoc)
scen_year <- "2023"
prev_scen_year <- as.character(as.numeric(scen_year) - 1)
common_year <- as.character(as.numeric(scen_year) - 4) # this depends on what the latest data year is; there's often a 3 year lag so the common year between the current and previous version year would be 4 years prior to the current version year; this time there is a 2 year lag but the common year is still 4 years prior at the moment
## Load FAO-specific user-defined functions
source(here('workflow/R/fao_fxn.R')) # function for cleaning FAO files
source(here('workflow/R/common.R')) # directory locations## This file makes it easier to process data for the OHI global assessment
## by creating the following objects:
##
## * dir_M = identifies correct file path to Mazu (internal server) based on your operating system
## * mollCRS = the crs code for the mollweide coordinate reference system we use in the global assessment
## * regions_shape() = function to load global shapefile for land/eez/high seas/antarctica regions
## * ohi_rasters() = function to load two rasters: global eez regions and ocean region
## * region_data() = function to load 2 dataframes describing global regions
## * rgn_syns() = function to load dataframe of region synonyms (used to convert country names to OHI regions)
## * low_pop() = function to load dataframe of regions with low and no human population
## * UNgeorgn = function to load dataframe of UN geopolitical designations used to gapfill missing dataSimultaneously read and process FAO commodities value and quantity data.
## NOTE: This can be run as a loop, but the "value" and "quant" datasets need to be run individually to make sure
## there are no problems (after this check, they can be looped for efficiency)
## describe where the raw data are located:
dir_fao_data <- file.path(dir_M, paste0('git-annex/globalprep/_raw_data/FAO_commodities/d', scen_year))
## list files included in d2023 folder (value and quant datasets)
files <- list.files(dir_fao_data, pattern=glob2rx('*.csv'), full.names=TRUE)
## To compare to old data:
# dir_fao_data <- file.path(dir_M, 'git-annex/globalprep/_raw_data/FAO_commodities/d2020')
# files <- list.files(dir_fao_data, pattern=glob2rx('*.csv'), full.names=T)
## loop
tic()
for (f in files){ # f = files[2]
cat(sprintf('\n\n\n====\nfile: %s\n', basename(f)))
d <- read.csv(f, check.names = FALSE, strip.white = TRUE, stringsAsFactors = FALSE) # stringsAsFactors=T
# checks names syntactically, strips leading and trailing whitespace, prevents conversion of characters to factors
## Specifies that units are tonnes if we are reading in the Commodities Quantity data csv, and usd if we are reading in the Commodities Value data csv
units <- c('tonnes','usd')[str_detect(f, c('quant','value'))] # detect unit name using lowercase American English
## gather into long format and clean up FAO-specific data foibles
## warning: attributes are not identical across measure variables; they will be dropped: this is fine
m <- d %>%
dplyr::select(-`Unit (Name)`) %>%
rename(country = `Reporting country (Name)`,
commodity = `Commodity (Name)`,
trade = `Trade flow (Name)`) %>%
rename_with(~ gsub("\\[", "", .)) %>%
rename_with(~ gsub("\\]", "", .)) %>%
pivot_longer(cols = -c(country, commodity, trade, Unit),
names_to = "year", values_to = "value")
## Include only the "Exports" data:
## 2022 - changed from "Export" to "Exports"
m <- m %>%
filter(trade == "Exports")
m <- m %>%
fao_clean_data_new() %>% # swaps out FAO-specific codes. NOTE: optional parameter 'sub_N' can be passed to control how an 'N' code is interpreted.
select(-trade, -Unit) %>% # eliminate 'trade' column
arrange(country, commodity, is.na(value), year)
## Products join: attach product categories from com2prod, and
## filter out all entries that do not match a product category.
## Note: commodity_lookup is user-defined function to compare
## commodities in data vs commodities in lookup table
## load lookup for converting commodities to products
com2prod <- read.csv(here(paste0('globalprep/prs_coral_harvest/v', scen_year, '/raw/commodities2products_weighting.csv')), na.strings='')
## version used in 2019:
## read.csv(here('globalprep/np/v2019/raw/commodities2products.csv'), na.strings='')
## version used in 2018:
## com2prod <- read.csv('raw/commodities2products.csv', na.strings='')
## version used in 2015: use when testing....
## com2prod <- read.csv('../v2014_test/commodities2products.csv', na.strings='')
## Check the current commodity-to-product lookup table. If necessary, make changes to "raw/commodities2products_weighting.csv"
## v2021: Missing a couple in each category.
# Make sure to examine the fish oil category thoroughly though. Some of the instances that are caught here are NOT fish oil. For example: "Whole lobsters Homarus spp, in shell or not, dried, salted/in brine, smoked, cooked by steaming/boiling in water" is including because "boiling" has "oil" in it. Add the appropriate missing commodities to commodities2products_weighting.csv
# Don't worry about any that show up under corals, shells, or sponges, since we don't include them in the assessment anyways. We are just keeping them in for now, because we use them as an example in the methods document (see below).
np_commodity_lookup(m, com2prod)
## inner_join will attach product names to matching commodities according to
## lookup table 'com2prod', and eliminate all commodities that do not appear in the lookup table.
m <- m %>%
inner_join(com2prod, by='commodity')
## Special case: user-defined function deals with
## breaking up Antilles into separate reported rgns
m <- split_regions(m) # used to use np_split_antilles(); this incorporates more than just antilles
## Some changes to region names that aren't working in name_2_rgn()
# m <- m %>%
# mutate(country = ifelse(country == "Côte d'Ivoire", "Ivory Coast", country)) %>%
# mutate(country = ifelse(country == "C<f4>te d'Ivoire ", "Ivory Coast", country)) %>%
# mutate(country = ifelse(country == "C\xf4te d'Ivoire", "Ivory Coast", country)) %>%
# mutate(country = ifelse(country == "Cura<e7>ao","Curacao", country)) %>%
# mutate(country = ifelse(country == "Curaçao","Curacao", country)) %>%
# mutate(country = ifelse(country == "Cura\xe7ao","Curacao", country)) %>%
# mutate(country = ifelse(country == "R\xe9union", "Reunion", country)) %>%
# mutate(country = ifelse(country == "Réunion", "Reunion", country)) %>%
# mutate(country = ifelse(country == "R<e9>union", "Reunion", country)) %>%
# mutate(country = ifelse(country == "Micronesia (Fed. States)", "Micronesia", country)) %>%
# mutate(country = ifelse(country == "Saint Helena/Asc./Trist.", "Saint Helena", country)) %>%
# mutate(country = ifelse(country == "Türkiye", "Turkey", country)) %>%
# mutate(country = ifelse(country == "Venezuela (Boliv Rep of)", "Venezuela", country)) %>%
# mutate(country = ifelse(country == "Netherlands (Kingdom of the)", "Netherlands", country)) %>%
# filter(country != "Azerbaijan") # landlocked, but not being removed by name_2_rgn?
# the above shouldn't be needed anymore in v2023 onwards due to updates to name_2_rgn
m_rgn <- name_2_rgn(df_in = m,
fld_name='country',
flds_unique=c('commodity', 'product', 'year'))
# v2023 removed countries:
# These data were removed for not having any match in the lookup tables:
#
# yugoslavia sfr
# 1 (not used and only has data in 1984-1988)
# v2023 duplicates: [1] "China" "China, Hong Kong SAR" "China, Macao SAR" "Guadeloupe"
# [5] "Martinique" "Montenegro" "Russian Federation" "Serbia and Montenegro"
# [9] "Sudan" "Sudan (former)" "Un. Sov. Soc. Rep." - these are all fixed in the group by and summarise below
## combine composite regions
## When summarizing the dataset, this function provides a modified way to sum the value column while maintaining NA values when both variables are NA (rather than turning to zero values). The function will sum non-NA values normally.
sum_function <- function(x) {
if (sum(is.na(x)) == length(x))
return(NA)
return(sum(x, na.rm = T))}
m_rgn <- m_rgn %>%
group_by(rgn_id, rgn_name, commodity, product, year) %>%
summarize(value = sum_function(value)) %>%
ungroup()
## units: rename value field to units based on filename
names(m_rgn)[names(m_rgn) == 'value'] <- units
## output to .csv - should create two csvs (tonnes.csv and usd.csv)
harvest_out <- sprintf(here(paste0('globalprep/prs_coral_harvest/v', scen_year, '/int/%s.csv')), units)
write.csv(m_rgn, harvest_out, row.names = FALSE, na = '')
}
toc()
# v2023: 152.934 sec elapsedRead in the tonnes and usd data that was completed in the Natural Products ornamentals dataprep (and above). Combining the quantity and value data and a bit of cleaning to remove data prior to first reporting year for coral commodities and regions.
## Read in quant dataset from intermediate folder
h_tonnes <- read.csv(here(paste0('globalprep/prs_coral_harvest/v', scen_year, '/int/tonnes.csv')))
## Read in value dataset from intermediate folder
h_usd <- read.csv(here(paste0('globalprep/prs_coral_harvest/v', scen_year, '/int/usd.csv')))
## concatenates h_tonnes and h_usd data
## h includes rgn_name, rgn_id, commodity, product, year, tonnes, usd.
h <- h_usd %>%
full_join(h_tonnes, by=c('rgn_name', 'rgn_id', 'commodity', 'product', 'year')) %>%
mutate(commodity = as.character(commodity)) %>%
arrange(rgn_id, product, commodity, year) # %>%
# dplyr::filter(product == "corals") ##filter for our commodities of interest
## clips out years prior to first reporting year, for each commodity per region
h <- h %>% np_harvest_preclip()Summary of gapfilling that is performed:
h <- h %>% np_harvest_gapflag()
## Adds flag for required gap-filling, based upon NAs in data.
## NOTE: Does not perform any gap-filling.
## At this point, h includes:
## rgn_name rgn_id commodity product year tonnes usd gapfill
## 'gapfill' will be in (zerofill, endfill, tbd, none)
data_check <- h %>% np_datacheck()
## for each commodity within each region, creates (but doesn't save...) summary info:
## num_years: the length of the data series for this commodity in this region
## usd_unique_nz: (or 'tns') number of unique non-zero values for usd or tonnes
## usd_na & tns_na: number of NA occurrences
## paired_obs: number of non-zero paired observations
## usd_unique_pairs: (or 'tns') within set of paired observations, count of unique usd and tonnes
## unique_pairs: lesser of usd_unique_pairs and tns_unique_pairs
## count_no_data: number of paired NAs - years with no value reported
h <- h %>% np_zerofill()
## for post-reporting years with NA for both tonnes and USD, fill zero -
## assumes that non-reporting indicates zero harvest to report.
## Also cross-fills zeros where one side is 0, other is NA (not flagged as gapfill)
h <- h %>% np_lowdata_filter()
## Exclude commodities (within a region) that have few non-zero data points.
## Optional parameter with default: nonzero_h_yr_min = 4
## NOTE: This filter has consequences for the regression, but also has meaning in terms of
## not inflicting a penalty on regions trying, and then stopping, an experimental harvest.
## Melanie's script to add a georegional ID tag based on country keys and IDs.
h <- h %>%
add_georegion_id()
h <- h %>% np_regr_fill(years_back = 10, vars = 'td', scope = 'rgn_id')
h <- h %>% np_regr_fill(vars = 'tdy', scope = 'georgn_id')
h <- h %>% np_regr_fill(vars = 'tdy', scope = 'global')
## np_regr_fill() is a generalized regression gapfill function. Parameters (with defaults):
## * years_back=50 (int): This determines how far back in the time series to include within the regression.
## * min_paired_obs=4 (int): This determines how many paired observations are required to attempt a regression.
## * scope = 'rgn_id' (str): ('rgn_id', 'georgn_id', 'global') Determines grouping scale for regression.
## * vars = 'tdy' (str): ('td', 'tdy') Determines model: (tonnes ~ usd) or (tonnes ~ usd + year) [and vice versa]
h <- h %>% np_end_fill()
## For final year of data, if both usd and tonnes originally reported as NA, pull forward
## values for usd and tonnes from the previous year. This should happen after regression fill.
h_comm <- h
## Store commodity-level data, before moving on to the product-level smoothing.
## Output gapfilling report to .csv files.
## Very few usd gapfilling, and none in recent years (data used to weight contributions), so will ignore this: gapfill=="r2_u_gr"
h_gap <- h %>%
mutate(gapfill = ifelse(gapfill == "r2_u_gr", "none", gapfill)) %>% # focusing only on tonnes gapfilling
select(rgn_id, commodity, product, year, gapfill) %>%
filter(product == "corals")
write.csv(h_gap, file.path(here(), paste0('globalprep/prs_coral_harvest/v', scen_year, '/output/prs_coral_gf.csv')), row.names = FALSE, na = '')Summarize each product per country per year, e.g., all corals in Albania in 2011. And, do some error checking.
h_prod <- h_comm %>%
filter(product == "corals") %>%
group_by(rgn_name, rgn_id, product, year) %>%
summarize(tonnes = sum(tonnes, na.rm = TRUE),
usd = sum(usd, na.rm = TRUE))
## Error-checking and table exports to see if there are duplicates
stopifnot(sum(duplicated(h_prod[ , c('rgn_id', 'product', 'year')])) == 0)Look at wide format with all commmodities and product subtotal (where commodity column value is “Z_TOTAL”), compared with the input data prior to summing.
h_x_tonnes <- h_comm %>%
bind_rows(mutate(h_prod, commodity='Z_TOTAL')) %>%
select(rgn_name, rgn_id, commodity, product, year, tonnes) %>%
arrange(rgn_name, product, commodity, year) %>%
pivot_wider(names_from = year, values_from = tonnes)
h_x_usd <- h_comm %>%
bind_rows(mutate(h_prod, commodity='Z_TOTAL')) %>%
select(rgn_name, rgn_id, commodity, product, year, usd) %>%
arrange(rgn_name, product, commodity, year) %>%
pivot_wider(names_from = year, values_from = usd)
## Check a random country and commodity
australia <- h_x_usd %>% filter(product == "corals", rgn_name == "Australia")
australia ## perfect
## Can open up in Excel to compare subtotals per country-product-year
write.csv(h_x_tonnes, paste0('globalprep/prs_coral_harvest/v', scen_year, '/int/coral_harvest_tonnes_wide.csv'), row.names = FALSE, na = 'NA')
write.csv(h_x_usd, paste0('globalprep/prs_coral_harvest/v', scen_year, '/int/coral_harvest_usd_wide.csv'), row.names = FALSE, na = 'NA')Determine rolling averages for tonnes and USD in order to determine peak values. This is based upon total harvests by product group, not individual commodity.
# Find max year in the summarized data table
year_max <- max(h_prod$year)
roll_prod <- h_prod %>%
arrange(rgn_id, product, year) %>%
group_by(rgn_id, product) %>%
mutate(
tonnes_rollmean = rollapply(tonnes, width=4, FUN=mean, align='right', partial=TRUE, na.rm=FALSE),
usd_rollmean = rollapply( usd, width=4, FUN=mean, align='right', partial=TRUE, na.rm=FALSE)) %>%
rename(
tonnes_orig = tonnes, # prevent overwriting of reported and gapfilled values
usd_orig = usd) %>% # prevent overwriting of reported and gapfilled values
mutate(
tonnes = ifelse(!is.na(tonnes_rollmean), tonnes_rollmean, tonnes_orig),
usd = ifelse(!is.na(usd_rollmean), usd_rollmean, usd_orig)) %>%
select(rgn_id, rgn_name, product, year, tonnes, usd, tonnes_orig, usd_orig)
write.csv(roll_prod, paste0("globalprep/prs_coral_harvest/v", scen_year, "/int/tonnes_coral_harvest.csv"), row.names = FALSE)## read in production harvest data
roll_prod <- read_csv(paste0("globalprep/prs_coral_harvest/v", scen_year, "/int/tonnes_coral_harvest.csv"))
## read in coral extent
coral_ext <- read_csv(paste0("globalprep/hab_coral/v", (as.numeric(scen_year) - 2), "/data/habitat_extent_coral_updated.csv")) %>% # UPDATE TO MOST RECENT VERSION
dplyr::select(-habitat, -year) %>%
filter(km2 != 0)
## read in coral health data
coral_health <- read_csv(paste0("globalprep/hab_coral/v", (as.numeric(scen_year) - 2), "/data/habitat_health_coral_updated.csv")) %>% # UPDATE TO MOST RECENT VERSION
dplyr::select(-habitat, -year)
# join together with the coral harvest data
coral_harvest <- roll_prod %>%
left_join(coral_ext, by = "rgn_id") %>%
left_join(coral_health, by = "rgn_id")
coral_harvest <- coral_harvest %>%
mutate(intensity = tonnes/km2)
ref = quantile(coral_harvest$intensity, probs = 0.95, na.rm = TRUE) ## find the 95th quantile for a reference point
coral_harvest <- coral_harvest %>%
mutate(pressure_no_health = ifelse(intensity > ref, 1, intensity/ref)) %>%
mutate(pressure_health = pressure_no_health*health) %>% ## calculate the pressure score
filter(!is.na(pressure_health)) %>%
dplyr::select(rgn_id, year, pressure_score = pressure_health)
# test <- coral_harvest %>%
# filter(is.na(km2))
#
# unique(test$rgn_id)write.csv(coral_harvest, paste0("globalprep/prs_coral_harvest/v", scen_year, "/output/prs_coral_harvest.csv"), row.names = FALSE)Datacheck
# get current version and previous version data and then combine for comparison
current_version_coral <- read_csv(paste0("globalprep/prs_coral_harvest/v", scen_year, "/output/prs_coral_harvest.csv")) %>%
rename("new_prs" = "pressure_score")
previous_version_coral <- read_csv(paste0("globalprep/prs_coral_harvest/v", prev_scen_year, "/output/prs_coral_harvest.csv")) %>%
rename("old_prs" = "pressure_score")
combined_coral <- previous_version_coral %>%
left_join(current_version_coral, by = c("rgn_id", "year"))
# if wanting to look at highest differences
changes <- combined_coral %>%
mutate(diff = abs(old_prs-new_prs))
# see what rows have NAs and take note if anything is odd
na_regions <- combined_coral[!complete.cases(combined_coral), ]
na_regions %>% select(rgn_id) %>% unique()
# v2023: no NAs
# check all data from previous version vs. all data from current version
plot(combined_coral$old_prs, combined_coral$new_prs,
xlab = paste0("Old Pressure All Years v", prev_scen_year),
ylab = paste0("New Pressure All Years v", scen_year))
abline(0,1, col = "red")
# explore some deviations v2023
previous_version_coral_latest_yr <- previous_version_coral %>%
filter(year == common_year)
current_version_coral_latest_yr <- current_version_coral %>%
filter(year == as.character(as.numeric(scen_year) - 2)) # set to latest year of data
combined_latest_yr <- previous_version_coral_latest_yr %>%
left_join(current_version_coral_latest_yr, by = "rgn_id") %>%
select(-year.x, -year.y)
plot(combined_latest_yr$old_prs, combined_latest_yr$new_prs,
xlab = paste0("Old Pressure Data Year ", common_year, " v", prev_scen_year),
ylab = paste0("New Pressure Data Year ", (as.numeric(scen_year) - 2), " v", scen_year))
abline(0,1, col = "red")
# compare within the new data itself 2020 vs. 2021 since both years are new v2023
current_version_coral_prior_yr <- current_version_coral %>%
filter(year == as.character(as.numeric(scen_year) - 3)) # set to 1 prior to latest year of data
combined_two_recent_yrs <- current_version_coral_prior_yr %>%
rename(prior_prs = new_prs) %>%
left_join(current_version_coral_latest_yr, by = "rgn_id") %>%
select(-year.x, -year.y)
plot(combined_two_recent_yrs$prior_prs, combined_two_recent_yrs$new_prs,
xlab = paste0("New Pressure Data Year ", (as.numeric(scen_year) - 3), " v", scen_year),
ylab = paste0("New Pressure Data Year ", (as.numeric(scen_year) - 2), " v", scen_year))
abline(0,1, col = "red")
# check the most recent year in common in both versions
combined_coral_common_yr <- combined_coral %>%
filter(year == common_year)
plot(combined_coral_common_yr$old_prs, combined_coral_common_yr$new_prs,
xlab = paste0("Old Pressure Data Year ", (as.numeric(prev_scen_year) - 3), " v", prev_scen_year),
ylab = paste0("New Pressure Data Year ", (as.numeric(scen_year) - 4), " v", scen_year))
abline(0,1, col = "red")
# glance at rolling production for one region
roll_prod <- read_csv(paste0("globalprep/prs_coral_harvest/v", scen_year, "/int/tonnes_coral_harvest.csv")) %>%
filter(rgn_id == 207)