knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'Figs/',
                      echo = TRUE, message = FALSE, warning = FALSE)

library(dplyr)
# library(rgdal)
library(raster)
library(here)
library(ggridges)
library(ggplot2)
library(tidyverse)
library(terra)

source(here('workflow/R/common.R'))

region_data()
rgns_global <- rgns_eez  %>%
  dplyr::select(rgn_id) %>%
  dplyr::filter(rgn_id != 213)

1 Source

We use the methods and code from Casey O’Hara: https://github.com/oharac/spp_risk_dists to prepare the data for the species subgoal.

Reference:

2 Summary

Spatial species range data and extinction risk data from IUCN is used to generate regional scores for the Species subgoal (part of the Biodiversity goal) and resilience layers.

Mean risk status per cell: Species ranges are converted to a global spatial raster of 10 km resolution.

The mean extinction risk for each cell, $\bar{R}_{cell}$, is calculated by averaging the IUCN extinction risk of the species with ranges overlapping the cell.

Risk is a scaled value representing the species extinction risk category: * ‘LC’ = 0.0, ‘NT’ = 0.2, ‘VU’ = 0.4, ‘EN’ = 0.6, ‘CR’ = 0.8, ‘EX’ = 1.0

\[\bar{R}_{cell} = \frac{\displaystyle\sum_{species}(Risk)}{n_{spp}}\]

Mean risk status per region: The mean extinction risk for a region, $\bar{R}_{SPP}$, is estimated by averaging the risk values of the raster cells falling within each OHI region, with each cell’s contribution weighted by the number of species in the cell.

Species goal model

The regional risk values are converted to species status scores by subtracting the risk values from 1 and rescaling so a risk value of $ $ receives a score of zero.

From Halpern et al (2012):

The target for the Species sub-goal is to have all species at a risk status of Least Concern. We scaled the lower end of the biodiversity goal to be 0 when 75% species are extinct, a level comparable to the five documented mass extinctions and would constitute a catastrophic loss of biodiversity.

\[X_{SPP} = \frac{((1 - \bar{R}_{SPP}) - 0.25)}{(1 - 0.25)} * 100%\]

where:

$X_{SPP}$ is Species goal status
$\bar{R}_{SPP}$ is mean extinction risk for a region

Trend is calculated using time series of risk categories based on current and past assessments.

3 Updates from previous assessment

Additional year of data

4 Data Sources

Reference:
- IUCN 2024. The IUCN Red List of Threatened Species. Version 2024-1. http://www.iucnredlist.org.
  - Shapefiles available from: https://www.iucnredlist.org/resources/spatial-data-download
  - Downloaded: August, 2024 *BirdLife International and Handbook of the Birds of the World (2020) Bird species distribution maps of the world. Version 2020.1. BirdLife International, Cambridge, UK and NatureServe, Arlington, USA. http://datazone.birdlife.org/species/requestdis.
  - Zipped shapefile available from BirdLife International.
  - Downloaded: Feb 26, 2021
Description: Shapefiles containing polygons of assessed species ranges; each shapefile represents all assessed species within a comprehensively-assessed (i.e. >90% assessed) taxonomic group.
Native data resolution: NA
Time range: NA
Format: Shapefile

Reference: IUCN and UNEP-WCMC (2021), The World Database on Protected Areas (WDPA) [On-line], February 2021. Cambridge, UK: UNEP-WCMC. Available at: www.protectedplanet.net.

Downloaded: August 15, 2022

Description: Shapefile of World Database on Protected Areas

Time range: 1819 - 2022; some protected areas do not have an associated “status year” and are reported as year 0.

Format: Shapefile

File location: Mazu:git-annex/globalprep/_raw_data/wdpa_mpa/d2021/WDPA_WDOECM_wdpa_shp/

5 Methods

There are several steps that need to be taken to get to this point in the data prep.

Here is an overview of the organization of files and data that are run prior to this:

5.0.1 Code

5.0.1.1 Run this first! Setup directory: `spp/v20XX/_setup`

In this directory are a sequence of files used to generate the bits and pieces that are later assembled into the rasters of biodiversity risk.

.Rmd files are sequenced with a prefix number (and letter) to indicate the order of operations. Briefly:

Pull information from the IUCN Red List API to determine an overall species list, habitat information, and current risk (conservation status).
- 1_set_up_iucn_habs_and_risk.Rmd
Pull information from API on risk from regional assessments; also recode the regions according to Marine Ecoregions (Spalding et al, 2007) for later spatialization.
- 2_set_up_iucn_risk_regional.Rmd
- 2a_set_up_regional_codes.Rmd
Pull historical assessment information from API for possible trend analysis.

3_set_up_iucn_trends.Rmd

Set up spatial layers in Gall-Peters, 100 km² cells. Layers include:
- cell ID (cells are sequentially numbered for combining with tabular data)
- ocean area
- marine protected area (classification, year of protection, proportion of protection)
- Exclusive Economic Zones (EEZ) and FAO fishing regions
- Marine Ecoregions of the World
- bathymetry
- NOTE: these layers are all saved in the spp_risk_dists/_spatial directory.
  - 4a_set_up_ocean_area_and_mpa_pct.Rmd
  - 4b_set_up_eez_meow_rasts.Rmd: You can skip this one.
Convert species range maps to rasters.
- For maps provided directly by IUCN, aggregate into multispecies files based on family. There is some cleaning done at this stage to fix problematic extents and attributes.
- From the list of all available maps, generate a master list of all mapped, assessed species for inclusion in the study.
- Rasterize each species to a .csv that includes cell ID and presence. A .csv format was used for file size and ease of reading and binding into dataframes.
  - 5b_generate_spp_map_list.Rmd
  - 5c_rasterize_spp_shps.Rmd
Aggregate individual species ranges into larger taxonomic groups, and summarize key variables (mean risk, variance of risk, number of species, etc) by group.
- Technically this is not necessary but makes it easier to quality check the process along the way, and supports mapping at the level of taxonomic group rather than the entire species list level.
- This process is done twice: once for uniform weighting and once for range-rarity weighting. Resulting files are saved separately.
  - 6a_aggregate_spp_ranges.Rmd
  - 6b_aggregate_rr_spp_ranges.Rmd

5.0.1.2 Then run this! Root directory: `v2024`

At this level there are several scripts, prefixed 1x_biodiversity_maps, that collate the various taxonomic group level files (generated in setup part 6) and summarize to the global level. These need to be run before spp_data_prep.Rmd! * 1a_biodiversity_maps_comp_assessed.Rmd * 1c_biodiversity_maps_all_spp.Rmd * spp_data_prep.Rmd

Note each creates a specific aggregation - comprehensively assessed species vs all available species; uniform vs range-rarity weighting.
The rasters generated in these scripts are saved in the _output folder.

5.0.2 Data and output files

The spp_risk_dists/_data folder contains tabular data about IUCN species used throughout the processing of this analysis. These files are generated by scripts in the setup directory.

The spp_risk_dists/_spatial folder contains general spatial data generated and/or used in the setup scripts. These include:

rasters for cell ID, EEZ ID, marine ecoregion ID, ocean area, and bathymetry masks.
tabular data of region names and lookups for IUCN regional assessment to marine ecoregion.
tabular data of marine protected area level/year/coverage to cell ID.
shapefiles used for map plotting from Natural Earth.

The spp_risk_dists/_output folder contains the rasters of biodiversity risk, species richness, variance of risk, etc generated from the scripts in the base directory.

5.1 Compare all vs. comprehensively assessed species

In the past, we have used all species with IUCN risk assessments to calculate the species subgoal. However, some of Casey’s work suggests it is better to use the taxa groups that have been comprehensively assessed by IUCN (> 90% of species assessed). The general concern is that IUCN tends to oversample species in the Atlantic portion of the ocean, relative to other regions. This is indicated by the larger number of species with IUCN status in this region. However, the Atlantic falls in line with the other regions when looking at the comprehensively assessed species.

# devtools::install_github("dill/beyonce")
library(beyonce)
cols <- beyonce_palette(129, 100, type = "continuous")

#2024
n_comp <- terra::rast(here("globalprep/spp/v2024/_output/n_spp_risk_raster_comp.tif"))
log_n_comp <- log(n_comp)
plot(log_n_comp, col=cols)

#2022
n_comp <- terra::rast(here("globalprep/spp/v2022/_output/n_spp_risk_raster_comp.tif"))
log_n_comp <- log(n_comp)
plot(log_n_comp, col=cols)

#2021
n_comp <- terra::rast(here("globalprep/spp/v2021/_output/n_spp_risk_raster_comp.tif"))
log_n_comp <- log(n_comp)
plot(log_n_comp, col=cols)


#2024
n_all <- terra::rast(here("globalprep/spp/v2024/_output/n_spp_risk_raster_all.tif"))
log_n_all <- log(n_all)
plot(log_n_all, col=cols)

#2022
n_all <- terra::rast(here("globalprep/spp/v2022/_output/n_spp_risk_raster_all.tif"))
log_n_all <- log(n_all)
plot(log_n_all, col=cols)

#2021
n_all <- terra::rast(here("globalprep/spp/v2021/_output/n_spp_risk_raster_all.tif"))
log_n_all <- log(n_all)
plot(log_n_all, col=cols)


prop_comp <- n_comp/n_all
plot(prop_comp, col=rev(cols))
#this makes sense. For some reason Damselfish and a few reptiles were classified as "non comprehensive", but they are actually comprehensive.

5.2 SPP: Status

5.2.1 Status: Extract average species risk for each region

For each cell, we multiply the average species risk by the number of species in order to weight each cells contribution by the number of species. We sum these values for each region and calculate: (average species risk * number species)/number of species

mean_risk_comp <- terra::rast(here("globalprep/spp/v2024/_output/mean_risk_raster_comp.tif"))

n_comp <- terra::rast(here("globalprep/spp/v2024/_output/n_spp_risk_raster_comp.tif"))

regions_ohi <- terra::rast(here("globalprep/spp/v2024/_spatial/eez_rast.tif"))

risk_stack_comp <- c(regions_ohi, mean_risk_comp, n_comp)
plot(risk_stack_comp)

risk_vals_comp <- terra::values(risk_stack_comp) %>%
  data.frame()
risk_vals_comp <- filter(risk_vals_comp, !is.na(eez_rast))
risk_vals_comp <- filter(risk_vals_comp, !is.na(mean_risk_raster_comp))

rgn_risk_comp <- risk_vals_comp %>%
  rowwise() %>%
  dplyr::mutate(risk_weight = mean_risk_raster_comp * n_spp_risk_raster_comp) %>%
  group_by(eez_rast) %>%
  summarize(rgn_risk_weight = sum(risk_weight),
            rgn_n_species = sum(n_spp_risk_raster_comp))


## This is the first time this object comes up in this script and its calling on itself ?? (found the above code from the spp_data_prep.Rmd in the v2021 root directory)
 rgn_risk_comp <- rgn_risk_comp %>%
  dplyr::rowwise() %>%
  dplyr::mutate(mean_risk = rgn_risk_weight/rgn_n_species) %>%
  dplyr::select(rgn_id = eez_rast, mean_risk)

5.2.2 Status: estimate for previous years

We use the trend data to estimate risk values for previous years (vs. using the same values for all assessment years). The change in species status across years is based on a linear model.

Trend is calculated using the same method as the risk calculation. For each cell, we multiply the average species trend by the number of species in order to weight each cell’s contribution by the number of species in the cell. We sum these values for each OHI region and calculate for each region: (average species trend * number species)/number of species

trend_comp <- terra::rast(here("globalprep/spp/v2024/_output/trend_raster_comp.tif"))

n_trend_comp <- terra::rast(here("globalprep/spp/v2024/_output/n_trend_raster_comp.tif"))

regions_ohi <- terra::rast(here("globalprep/spp/v2024/_spatial/eez_rast.tif"))

trend_stack_comp <- c(regions_ohi, trend_comp, n_trend_comp)
trend_vals_comp <- terra::values(trend_stack_comp) %>%
  data.frame()

trend_vals_comp <- filter(trend_vals_comp, !is.na(eez_rast))
trend_vals_comp <- filter(trend_vals_comp, !is.na(trend_raster_comp))

rgn_trend_comp <- trend_vals_comp %>%
  rowwise() %>%
  dplyr::mutate(trend_weight = trend_raster_comp * n_trend_raster_comp) %>%
  group_by(eez_rast) %>%
  summarize(rgn_trend_weight = sum(trend_weight),
            rgn_n_species = sum(n_trend_raster_comp)) %>%
  rename(rgn_id = eez_rast) %>%
  ungroup()

rgn_trend_comp <- rgn_trend_comp %>%
  dplyr::rowwise() %>%
  dplyr::mutate(mean_trend = rgn_trend_weight/rgn_n_species) %>%
  dplyr::select(rgn_id, mean_trend)

5.2.3 Status: Get yearly risk scores based on trend

We estimate previous risk for each region, using the trend data. We assume change in risk is linear.

assess_years <- 2012:2023
years <- expand.grid(rgn_id = unique(rgn_risk_comp$rgn_id), year=assess_years)

# this is what the trend will be multiplied by to get a risk estimate for each year:
year_multiplier <- data.frame(year=assess_years, multiplier = rev(0:(length(assess_years)-1))) 

rgn_risk_comp_yrs <- rgn_risk_comp %>%
  dplyr::left_join(rgn_trend_comp, by = "rgn_id") %>%
  dplyr::left_join(years, by = "rgn_id") %>%
  dplyr::left_join(year_multiplier, by="year") %>%
  dplyr::rowwise() %>%
  dplyr::mutate(mean_risk_per_year = mean_risk - mean_trend*multiplier) %>%
  dplyr::select(rgn_id, year, mean_risk = mean_risk_per_year)

5.2.4 Status: Converting regional mean risk to status

We rescale the data so a risk factor of 0.75 is equal to zero.

rgn_status <- rgn_risk_comp_yrs %>%
  mutate(spp_status = (0.75 - mean_risk)/0.75)

5.2.5 Status: Gapfill missing regions

Region 232 (Bosnia) does not have a value, which is not surprising because their coast is super small and results are erratic for this region. We gapfill with surrounding regions.

status_gf <- rgns_global %>%
  left_join(rgn_status) %>%
    dplyr::select(-mean_risk)
summary(status_gf)
filter(status_gf, is.na(spp_status))

croatia <- filter(status_gf, rgn_id == 187)
mont <- filter(status_gf, rgn_id == 186) 

bosnia <- bind_rows(croatia, mont) %>%
  group_by(year) %>%
  summarize(spp_status = mean(spp_status)) %>%
  mutate(rgn_id = 232)

status_gf <- status_gf %>%
  filter(rgn_id !=232) %>%
  bind_rows(bosnia)

5.2.6 Status: Final formatting for ohi-global

status <- status_gf %>%
  dplyr::select(rgn_id, year, score = spp_status)
dim(status)   #220*length(assess_years)
summary(status) # should be no NA values

write.csv(status, here("globalprep/spp/v2024/output/sp_status_global.csv"), row.names=FALSE)

5.2.7 Status: Compare to last year

status <- read_csv(here("globalprep/spp/v2024/output/sp_status_global.csv"))

old_spp <- read.csv(here("globalprep/spp/v2022/output/sp_status_global.csv")) %>%
  filter(year == max(year)) %>%
  rename(old_score = score) %>%
  mutate(year = 2022) %>%
  left_join(status) %>%
  rename(new_score = score)

plot(old_spp$old_score, old_spp$new_score, xlab = "old", ylab = "new")
abline(0,1, col="red")

old_spp_gather <- old_spp %>%
  dplyr::select(rgn_id, old_score, new_score) %>%
  tidyr::gather("assessment", "score", -1) %>%
  dplyr::filter(rgn_id <= 250)

ggplot(old_spp_gather, aes(y=assessment, x=score)) + 
  geom_density_ridges()

ggplot(old_spp, aes(x = old_score, y = new_score)) + 
  geom_point() + 
  geom_abline(intercept = 0)

5.3 SPP: Trend

5.3.1 Trend: calculating

Getting proportional trend requires the status data (trend/status). Proportional trend is multiplied by 5 to get estimated change in five years.

# proportional trend requires status data
status <- read.csv(here("globalprep/spp/v2024/output/sp_status_global.csv")) %>%
  dplyr::filter(year==max(year)) %>%
  dplyr::select(rgn_id, score)

# Calculated in above section: Trend data
rgn_trend_score <- rgn_trend_comp %>%
  mutate(spp_trend_adj = -mean_trend/0.75) %>%  # puts in comparable units to status
  left_join(status, by="rgn_id") %>%
  dplyr::mutate(trend_score = spp_trend_adj/score * 5)

5.3.2 Trend: Gapfilling missing data

Check there are data for every region. Region 232 (Bosnia) does not have a value which is not surprising because their coast is super small and results are erratic for this region. We estimate this using the mean of the 2 surrounding regions.

trend <- rgns_global %>%
  left_join(rgn_trend_score) 

summary(trend)

filter(trend, is.na(trend_score))

croatia <- filter(trend, rgn_id == 187)
mont <- filter(trend, rgn_id == 186)
bosnia <- mean(c(croatia$trend_score, mont$trend_score))

trend$trend_score[trend$rgn_id == 232] <- bosnia

5.4 Trend: Final formatting for ohi-global

trend <- trend %>%
  dplyr::select(rgn_id, score = trend_score)
dim(trend) # should be 220
summary(trend) # should be no NAs

write.csv(trend, here("globalprep/spp/v2024/output/sp_trend_global.csv"), row.names=FALSE)

5.4.1 Trend: Compare to last year

Previously we used proxy data rather than actual change in IUCN status over time. Our previous method overestimated the magnitude of the trend. It is not surprising there is poor correlation with trend estimates in previous years (before 2019), but it is reassuring the the values mainly fall in the same quadrant, and trends are largely the same now (2019 - present)

trend <- read.csv(here("globalprep/spp/v2024/output/sp_trend_global.csv"))  

old_spp <- read.csv(here("globalprep/spp/v2022/output/sp_trend_global.csv")) %>%
  rename(old_score = score) %>%
  left_join(trend) 

plot(old_spp$old_score, old_spp$score)# xlim=c(-0.35, 0.05))
abline(h=0, col="red")
abline(v=0, col="red")
abline(0,1, col="red")

ggplot(old_spp, aes(x = old_score, y = score)) +
  geom_point() + 
  geom_abline(intercept = 0) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0)

5.4.2 Status/trend: save a record of gapfilling

status <- read.csv(here("globalprep/spp/v2024/output/sp_status_global.csv")) %>%
  mutate(gapfilled = ifelse(rgn_id == 232, 1, 0)) %>%
  mutate(method = ifelse(rgn_id == 232, "mean of neighbors", NA)) %>%
  dplyr::select(rgn_id, year, gapfilled, method)
write.csv(status, here("globalprep/spp/v2024/output/sp_status_global_gf.csv"), row.names=FALSE)          

trend <- read.csv(here("globalprep/spp/v2024/output/sp_trend_global.csv")) %>%
  mutate(gapfilled = ifelse(rgn_id == 232, 1, 0)) %>%
  mutate(method = ifelse(rgn_id == 232, "mean of neighbors", NA)) %>%
  dplyr::select(rgn_id, gapfilled, method)
write.csv(trend, here("globalprep/spp/v2024/output/sp_trend_global_gf.csv"), row.names=FALSE)

5.5 Resilience data

We use species condition data as a resilience measure as well. We also calculate species condition at 3nm of shoreline, because for some goals, nearshore species condition is the relevant metric.

5.5.1 Resilience: Prepare rasters for 3nm extraction

We reproject the data to have higher resolution in order to more easily extract the data at the 3nm scale.
We modify the method a bit from above due to size of the rasters.

#3nm raster file
rgns <- raster(file.path(dir_M, "git-annex/globalprep/spatial/v2018/rgns_3nm_offshore_mol.tif"))
plot(rgns)

# relevant species files
mean_risk_comp <- raster(here("globalprep/spp/v2024/_output/mean_risk_raster_comp.tif"))
plot(mean_risk_comp)
n_comp <- raster(here("globalprep/spp/v2024/_output/n_spp_risk_raster_comp.tif"))
plot(n_comp)


risk_x_n <- mean_risk_comp*n_comp
plot(risk_x_n)

# project rasters to moll
# saved in Mazu:spp/v2022
raster::projectRaster(risk_x_n, rgns, method="ngb", overwrite=TRUE, 
              filename=file.path(dir_M, "git-annex/globalprep/spp/v2024/int/risk_x_n_comp_mol.tif"), progress = "text")

projectRaster(n_comp, rgns, method="ngb", overwrite=TRUE, 
              filename=file.path(dir_M, "git-annex/globalprep/spp/v2024/int/n_comp_mol.tif"), progress = "text")

5.5.2 Resilience: Extract data

Extract species risk data that corresponds to 3nm regions.

#3nm raster file
rgns <- rast(file.path(dir_M, "git-annex/globalprep/spatial/v2018/rgns_3nm_offshore_mol.tif"))
# plot(rgns)

risk_x_n_mol <- rast(file.path(dir_M, "git-annex/globalprep/spp/v2024/int/risk_x_n_comp_mol.tif"))
# plot(risk_x_n_mol)
n_mol <- rast(file.path(dir_M, "git-annex/globalprep/spp/v2024/int/n_comp_mol.tif"))
# plot(n_mol)

risk_stack <- c(rgns, risk_x_n_mol, n_mol)

risk_df <- terra::values(risk_stack)
 
risk_df_2 <- as.data.frame(risk_df) %>%
  filter(!is.na(rgns_3nm_offshore_mol)) %>%
  filter(!is.na(risk_x_n_comp_mol))


rgn_3nm_risk <- risk_df_2 %>%
  rowwise() %>%
  group_by(rgns_3nm_offshore_mol) %>%
  summarize(risk_x_n_comp_mol = sum(risk_x_n_comp_mol),
            n_comp_mol = sum(n_comp_mol)) %>%
  ungroup() %>%
  rename(zone = 1) %>%
  dplyr::mutate(rgn_wt_risk = risk_x_n_comp_mol/n_comp_mol) %>%
  dplyr::select(rgn_id = zone, rgn_wt_risk)

5.5.3 Resilience: estimate for previous years

We use the trend data to estimate risk values for previous years (vs. using the same values for all assessment years). The change in species status across years is based on a linear model.

5.5.3.1 Resilience: Prepare rasters for 3nm extraction (trend to estimate previous years’ data)

We reproject the data to have higher resolution in order to more easily extract the data at the 3nm scale.

rgns <- raster(file.path(dir_M, "git-annex/globalprep/spatial/v2018/rgns_3nm_offshore_mol.tif"))

trend_comp <- raster::raster(here("globalprep/spp/v2024/_output/trend_raster_comp.tif"))

n_trend_comp <- raster::raster(here("globalprep/spp/v2024/_output/n_trend_raster_comp.tif"))

trend_x_n <- trend_comp*n_trend_comp

# project rasters to moll
# saved in Mazu:spp/v2022
projectRaster(trend_x_n, rgns, method="ngb", overwrite=TRUE, 
              filename=file.path(dir_M, "git-annex/globalprep/spp/v2024/int/trend_x_n_comp_mol.tif"),
              progress="text")
projectRaster(n_trend_comp, rgns, method="ngb", over=TRUE, 
              filename=file.path(dir_M, "git-annex/globalprep/spp/v2024/int/n_trend_comp_mol.tif"),
              progress="text")

5.5.4 Resilience: Extract trend data (used to estimate previous years’ data)

Extract species risk data that corresponds to 3nm regions.

#3nm raster file
rgns_3nm <- rast(file.path(dir_M, "git-annex/globalprep/spatial/v2018/rgns_3nm_offshore_mol.tif"))
plot(rgns_3nm)

trend_3nm_mol <- rast(file.path(dir_M, "git-annex/globalprep/spp/v2024/int/trend_x_n_comp_mol.tif"))
n_trend_3nm_mol <- rast(file.path(dir_M, "git-annex/globalprep/spp/v2024/int/n_trend_comp_mol.tif"))


trend_stack <- c(rgns_3nm, trend_3nm_mol, n_trend_3nm_mol)

# trend_df <- terra::zonal(trend_stack, rgns_3nm, fun='sum')

trend_df <- terra::values(trend_stack)
 
trend_df_2 <- as.data.frame(trend_df) %>%
  filter(!is.na(rgns_3nm_offshore_mol)) %>%
  filter(!is.na(trend_x_n_comp_mol))


rgn_3nm_trend <- trend_df_2 %>%
  data.frame() %>%
  rowwise() %>%
  group_by(rgns_3nm_offshore_mol) %>%
  summarize(trend_x_n_comp_mol = sum(trend_x_n_comp_mol),
            n_trend_comp_mol = sum(n_trend_comp_mol)) %>%
  ungroup() %>%
  rename(zone = 1) %>%
  dplyr::mutate(rgn_wt_trend = trend_x_n_comp_mol/n_trend_comp_mol) %>%
  dplyr::select(rgn_id = zone, rgn_wt_trend)

5.5.5 Resilience: Calculate yearly risk scores based on trend

We estimate previous risk for each region, using the trend data. We assume change in risk is linear.

assess_years <- 2012:2023
years <- expand.grid(rgn_id = unique(rgn_3nm_trend$rgn_id), year=assess_years)

# this is what the trend will be multiplied by to get a risk estimate for each year:
year_multiplier <- data.frame(year=assess_years, multiplier = rev(0:(length(assess_years)-1))) 

rgn_risk_3nm <- rgn_3nm_risk %>%
  left_join(rgn_3nm_trend, by = "rgn_id") %>%
  left_join(years, by = "rgn_id") %>%
  left_join(year_multiplier, by="year") %>%
  rowwise() %>%
  mutate(mean_risk_per_year = rgn_wt_risk - rgn_wt_trend*multiplier) %>%
  dplyr::select(rgn_id, year, mean_risk = mean_risk_per_year)

5.5.6 Resilience: Converting regional mean risk to status

We rescale the data so a risk factor of 0.75 is equal to zero.

rgn_3nm_res <- rgn_risk_3nm %>%
  mutate(spp_status = (0.75 - mean_risk)/0.75)

# quick check
hist(rgn_3nm_res$spp_status)

5.5.7 Resilience: Gapfill missing regions

Region 19 (Tuvalu) does not have a value. This is an island. We gapfill with the value from the entire eez.

res_gf <- rgns_global %>%
  left_join(rgn_3nm_res) %>%
    dplyr::select(-mean_risk)

summary(res_gf)
filter(res_gf, is.na(spp_status))

# get eez value:
eez_status <- read.csv(here("globalprep/spp/v2024/output/sp_status_global.csv")) %>%
  filter(rgn_id == 19) %>%
  rename(spp_status = score)

res_gf <- res_gf %>%
  filter(!is.na(spp_status)) %>%
  bind_rows(eez_status)


summary(res_gf)

5.5.8 Resilience: Final formatting for ohi-global

resilience <- res_gf %>%
  dplyr::select(rgn_id, year, score = spp_status)

summary(resilience) # should be no NA values

write.csv(resilience, here("globalprep/spp/v2024/output/sp_status_3nm.csv"), row.names=FALSE)

5.5.9 Resilience: save a record of gapfilling

res <- read.csv(here("globalprep/spp/v2024/output/sp_status_3nm.csv")) %>%
  mutate(gapfilled = ifelse(rgn_id == 19, 1, 0)) %>%
  mutate(method = ifelse(rgn_id == 19, "eez scale data used", NA)) %>%
  dplyr::select(rgn_id, year, gapfilled, method)
write.csv(res, here("globalprep/spp/v2024/output/sp_status_3nm_gf.csv"), row.names=FALSE)

5.5.10 Resilience: Compare

Compared to the entire EEZ, most (but not all) coastal areas have higher species condition scores. The correlation with last year is better than it was last year.

# compare to eez values
eez_status <- read.csv(here("globalprep/spp/v2024/output/sp_status_global.csv")) %>%
  filter(year == max(year)) %>%
  dplyr::select(rgn_id, eez_score = score)

status <- read.csv(here("globalprep/spp/v2024/output/sp_status_3nm.csv")) %>%
  dplyr::select(rgn_id, year, nm3_score = score) %>%
  filter(year == max(year)) %>%
  left_join(eez_status, by = "rgn_id")

plot(status$nm3_score, status$eez_score)
abline(0,1)

ggplot(status, aes(x = nm3_score, y = eez_score)) + 
  geom_point() + 
  geom_abline(intercept = 0) +
  labs(x = "3nm score", y = "eez score")

# compare to last year's values
status <- read.csv(here("globalprep/spp/v2024/output/sp_status_3nm.csv")) %>%
  filter(year == max(year)) %>%
  dplyr::select(-year)

old_spp <- read.csv(here("globalprep/spp/v2022/output/sp_status_3nm.csv")) %>%
  filter(year == max(year)) %>%
  rename(old_score = score) %>%
  # mutate(year = 2020) %>%
  dplyr::select(-year) %>%
  left_join(status) %>%
  rename(new_score = score) %>%
  mutate(difference = new_score - old_score) %>%
  left_join(rgns_eez)

plot(old_spp$old_score, old_spp$new_score)
abline(0,1, col="red")

ggplot(old_spp, aes(x = old_score, y = new_score)) + 
  geom_point() + 
  geom_abline(intercept = 0) +
  labs(x = "v2022 resilience", y = "v2024 resilience")

old_spp_gather <- old_spp %>%
  dplyr::select(rgn_id, old_score, new_score) %>%
  tidyr::gather("assessment", "score", -1) %>%
  filter(rgn_id <= 250)

ggplot(old_spp_gather, aes(y=assessment, x=score)) + 
  geom_density_ridges()

check_3nm_old <- old_spp %>%
  filter(new_score > 0.85, old_score > 0.79, old_score <0.87) %>%
  mutate(diff = old_score - new_score)

OHI 2022 - Species subgoal (Biodiversity goal and resilience)

Compiled on Tue Aug 27 17:19:37 2024 by dduncan

1 Source

2 Summary

3 Updates from previous assessment

4 Data Sources

5 Methods

5.0.1 Code

5.0.1.1 Run this first! Setup directory: `spp/v20XX/_setup`

5.0.1.2 Then run this! Root directory: `v2024`

5.0.2 Data and output files

5.1 Compare all vs. comprehensively assessed species

5.2 SPP: Status

5.2.1 Status: Extract average species risk for each region

5.2.2 Status: estimate for previous years

5.2.3 Status: Get yearly risk scores based on trend

5.2.4 Status: Converting regional mean risk to status

5.2.5 Status: Gapfill missing regions

5.2.6 Status: Final formatting for ohi-global

5.2.7 Status: Compare to last year

5.3 SPP: Trend

5.3.1 Trend: calculating

5.3.2 Trend: Gapfilling missing data

5.4 Trend: Final formatting for ohi-global

5.4.1 Trend: Compare to last year

5.4.2 Status/trend: save a record of gapfilling

5.5 Resilience data

5.5.1 Resilience: Prepare rasters for 3nm extraction

5.5.2 Resilience: Extract data

5.5.3 Resilience: estimate for previous years

5.5.3.1 Resilience: Prepare rasters for 3nm extraction (trend to estimate previous years’ data)

5.5.4 Resilience: Extract trend data (used to estimate previous years’ data)

5.5.5 Resilience: Calculate yearly risk scores based on trend

5.5.6 Resilience: Converting regional mean risk to status

5.5.7 Resilience: Gapfill missing regions

5.5.8 Resilience: Final formatting for ohi-global

5.5.9 Resilience: save a record of gapfilling

5.5.10 Resilience: Compare

OHI 2022 - Species subgoal (Biodiversity goal and resilience)

Compiled on Tue Aug 27 17:19:37 2024 by dduncan

1 Source

2 Summary

3 Updates from previous assessment

4 Data Sources

5 Methods

5.0.1 Code

5.0.1.1 Run this first! Setup directory: spp/v20XX/_setup

5.0.1.2 Then run this! Root directory: v2024

5.0.2 Data and output files

5.1 Compare all vs. comprehensively assessed species

5.2 SPP: Status

5.2.1 Status: Extract average species risk for each region

5.2.2 Status: estimate for previous years

5.2.3 Status: Get yearly risk scores based on trend

5.2.4 Status: Converting regional mean risk to status

5.2.5 Status: Gapfill missing regions

5.2.6 Status: Final formatting for ohi-global

5.2.7 Status: Compare to last year

5.3 SPP: Trend

5.3.1 Trend: calculating

5.3.2 Trend: Gapfilling missing data

5.4 Trend: Final formatting for ohi-global

5.4.1 Trend: Compare to last year

5.4.2 Status/trend: save a record of gapfilling

5.5 Resilience data

5.5.1 Resilience: Prepare rasters for 3nm extraction

5.5.2 Resilience: Extract data

5.5.3 Resilience: estimate for previous years

5.5.3.1 Resilience: Prepare rasters for 3nm extraction (trend to estimate previous years’ data)

5.5.4 Resilience: Extract trend data (used to estimate previous years’ data)

5.5.5 Resilience: Calculate yearly risk scores based on trend

5.5.6 Resilience: Converting regional mean risk to status

5.5.7 Resilience: Gapfill missing regions

5.5.8 Resilience: Final formatting for ohi-global

5.5.9 Resilience: save a record of gapfilling

5.5.10 Resilience: Compare

5.0.1.1 Run this first! Setup directory: `spp/v20XX/_setup`

5.0.1.2 Then run this! Root directory: `v2024`