OHI Science | Citation policy
This script prepares scores (status and trend) for Iconic Species in each global coastal region. For each iconic marine species, the countries of occurrence are pulled from the IUCN API. Extinction risk categories for each species are pulled based on current and past assessments; by tracking the assessed extinction risk over time, we can understand the trends of extinction risk for iconic species directly rather than using the “population trend” method from prior OHI assessments.
The Iconic Species sub-goal model calculates a region’s status based upon an unweighted average of species health for all ‘iconic’ species found within each reporting region.
From Halpern et al (2012):
Iconic species are those that are relevant to local cultural identity through a species’ relationship to one or more of the following: 1) traditional activities such as fishing, hunting or commerce; 2) local ethnic or religious practices; 3) existence value; and 4) locally-recognized aesthetic value (e.g., touristic attractions/common subjects for art such as whales). Habitat-forming species are not included in this definition of iconic species, nor are species that are harvested solely for economic or utilitarian purposes (even though they may be iconic to a sector or individual). …
Ultimately, almost any species can be iconic to someone, and so the intent with this goal was to focus on those species widely seen as iconic within a country, and iconic from a cultural or existence value (rather than for a livelihoods or extractive reason). …
The reference point is to have the risk status of all assessed species as Least Concern (i.e., a goal score = 1.0)
The Status of this sub-goal (XICO) is then the % of iconic species in each threat category (as defined by the IUCN Red List), such that:
\[X_{ICO} = \frac{\displaystyle\sum_{category}S_{cat}*w_{cat}}{\displaystyle\sum_{category}S_{cat}}\]
where for each IUCN threat category:
ICO trend is calculated in a similar manner, but weightings are assigned according to IUCN population trend: ‘Decreasing’ = -0.5, ‘Stable’ = 0.0, ‘Increasing’ = +0.5.
No changes have been made to the data prep methods for this goal, since the 2016 assessment, though there were some slight changes to improve the code. API key is read from a file, the ‘get from API’ function now includes error checking (it reports a little on what the error is, whether it’s a zero-length data frame or a try error, etc) and multiple attempts (if it gets an error it tries again a few times), and there is also a delay parameter to space out requests to the API.
Changes were made to the iconic species trend calculation. Trends are now calculated (in ohi-global, functions.R) using status values based only on species with two or more IUCN evaluations in the last 20 years. This change was made because including species which have not been re-evaluated for a long time tends to positively bias trend values, or at least dampen the trends by essentially including as zeros, values that are not truely known to be zero… To accommodate this change in the trend calculation, years 1992-present are included in the prepped data as well as in the ohi-global scenario_data_years.csv.
List of iconic species:
Species native country information:
Using the IUCN API, we accessed the full IUCN species list at http://apiv3.iucnredlist.org/api/v3/speciescount?token=
iucn_sid | kingdom | phylum | class | order | family | genus | sciname | population | category
Get all pages and bind into total species list, using the IUCN API.
api_file <- file.path(dir_anx, 'api_key.csv')
api_key <- scan(api_file, what = 'character') # fb71ae836f415f04f41176f6d30c4a9e4cea620d46b9e5021bf2fb142ea51bf5
spp_list_from_api_file <- file.path(dir_goal, 'int/spp_list_from_api.csv')
reload <- TRUE
if(!file.exists(spp_list_from_api_file) | reload) {
message('Using API to create full species list from scratch')
spp_npage_url <- sprintf('http://apiv3.iucnredlist.org/api/v3/speciescount?token=%s', api_key)
n_spp <- fromJSON(spp_npage_url) %>%
.$count %>% as.integer()
n_pages <- ceiling(n_spp/10000)
spp_page_url <- 'http://apiv3.iucnredlist.org/api/v3/species/page/%s?token=%s'
spp_df_all <- mc_get_from_api(spp_page_url, c(0:(n_pages - 1)), api_key) # mc_get_from_api function is defined in the ico_fxn.R script
spp_df_all <- spp_df_all %>%
dplyr::select(-infra_rank, -infra_name, -count, -page) %>%
rename(iucn_sid = taxonid, sciname = scientific_name) %>%
setNames(names(.) %>%
str_replace('_name', ''))
write_csv(spp_df_all, spp_list_from_api_file)
} else {
message('reading file of API species list: \n ', spp_list_from_api_file)
spp_df_all <- read_csv(spp_list_from_api_file)
}The list of Iconic Species is based upon the original ICO list generated in 2011, using species identified as globally iconic (WWF Flagship species and Priority species) or regionally iconic (based upon WWF regional/local priority species and nation-specific lists).
ico_list_raw <- get_ico_list(reload = TRUE) # get_ico_list is defined in ico_fxn.R
## head of ico_list_raw for v2018:
## comname sciname ico_gl ico_rgn_id
## <chr> <chr> <lgl> <int>
## 1 Blue Shark Prionace glauca T NA
## 2 Whale Shark Rhincodon typus T NA
## 3 Shortfin Mako Isurus oxyrinchus T NA
## 4 Olive Ridley Turtle Lepidochelys olivacea T NA
## 5 Irrawaddy Dolphin Orcaella brevirostris T NA
## 6 Humphead wrasse, Napoleon Wrasse Cheilinus undulatus T NA
## ico_rgn_id: rgn_id in which species is iconic by regional/national lists;
## if globally iconic, ico_rgn_id <- NA
write_csv(ico_list_raw, file.path(dir_goal, 'raw/ico_list_raw.csv'))DT::datatable(ico_list_raw, caption = 'Iconic Species List')Filtering the complete IUCN species list to include only the identified Iconic Species, we then use the IUCN API to access the list of countries in which each species occurs, from http://apiv3.iucnredlist.org/api/v3/species/countries/id/
spp_df_all <- read_csv(file.path(dir_goal, 'int/spp_list_from_api.csv'))
ico_list_raw <- read_csv(file.path(dir_goal, 'raw/ico_list_raw.csv')) # static list; same each year
spp_ico <- spp_df_all %>%
filter(sciname %in% ico_list_raw$sciname)
spp_missing <- ico_list_raw %>%
filter(!sciname %in% spp_ico$sciname) # for v2018, missing only Sprattus sprattus
ico_list <- ico_list_raw %>%
left_join(spp_ico %>%
select(iucn_sid, sciname, subpop = population, cat = category),
by = 'sciname') %>%
filter(!is.na(iucn_sid))
write_csv(ico_list, file.path(dir_goal, 'int/ico_list_prepped.csv')) # 144 in v2018For each of these species, use the IUCN API to gather a list of countries in which it is present.
## for each species ID, get country list
ico_country_url <- 'http://apiv3.iucnredlist.org/api/v3/species/countries/id/%s?token=%s'
ico_spp_countries <- mc_get_from_api(ico_country_url, ico_list$iucn_sid, api_key)
rgn_iucn2ohi <- read_csv(file.path(dir_goal, 'raw/rgns_iucn2ohi.csv'))
ico_spp_rgn_raw <- ico_spp_countries %>%
select(-code, -count, iucn_sid = name, iucn_rgn_name = country) %>%
mutate(iucn_sid = as.integer(iucn_sid),
iucn_rgn_name = str_trim(iucn_rgn_name)) %>%
left_join(rgn_iucn2ohi,
by = 'iucn_rgn_name')
## error check on region name matching; landlocked --> not matched
non_match <- ico_spp_rgn_raw %>%
filter(is.na(ohi_rgn_name))
if(nrow(non_match) > 0) {
cat('The following IUCN countries did not match with OHI region names:\n ')
print(paste(non_match$iucn_rgn_name %>% unique(), collapse = ', '))
}
ico_spp_rgn_raw <- ico_spp_rgn_raw %>%
rename(rgn_name = ohi_rgn_name) %>%
select(-iucn_rgn_name) %>%
filter(!is.na(rgn_id)) %>%
distinct()
write_csv(ico_spp_rgn_raw, file.path(dir_goal, 'raw/ico_spp_rgn_raw.csv'))ico_spp_rgn_raw <- read_csv(file.path(dir_goal, 'raw/ico_spp_rgn_raw.csv'))
ico_list <- read_csv(file.path(dir_goal, 'int/ico_list_prepped.csv'))
ico_spp_rgn_prepped <- ico_spp_rgn_raw %>%
left_join(ico_list, by = 'iucn_sid')
## filter this for species who are global (so all instances are iconic)
## OR ico_rgn_id matches rgn_id (so locally iconic matches with location)
ico_spp_rgn_prepped <- ico_spp_rgn_prepped %>%
filter(ico_gl == TRUE | ico_rgn_id == rgn_id)
write_csv(ico_spp_rgn_prepped, file.path(dir_goal, 'int/ico_spp_rgn_prepped.csv'))
ico_spp_rgn_prepped <- ico_spp_rgn_prepped %>%
select(iucn_sid, sciname, comname, presence, origin, distribution_code, cat, ico_gl, rgn_id, rgn_name, subpop, ico_rgn_id, api_error, param_id)`ico_spp_rgn_prepped` saved at `~/github/ohiprep/globalprep/ico/v2016/int/ico_spp_rgn_prepped.csv`DT::datatable(ico_spp_rgn_prepped, width = 400)We accessed the IUCN API to determine past IUCN assessments for each of the identified iconic species: http://apiv3.iucnredlist.org/api/v3/species/history/id/
Each assessment includes a year and an extinction risk, along with additional information on the assessment.
ico_past_assess_url <- 'http://apiv3.iucnredlist.org/api/v3/species/history/id/%s?token=%s'
ico_list <- read_csv(file.path(dir_goal, 'int/ico_list_prepped.csv'))
## for each species ID, get past assessments
ico_assess_raw <- mc_get_from_api(ico_past_assess_url, ico_list$iucn_sid, api_key)
ico_assess_raw <- ico_assess_raw %>%
rename(iucn_sid = name) %>%
mutate(iucn_sid = as.integer(iucn_sid),
year = as.integer(year)) %>%
left_join(ico_list %>%
select(iucn_sid, sciname) %>%
distinct(),
by = 'iucn_sid')
write_csv(ico_assess_raw, file.path(dir_goal, 'raw/ico_assessments_raw.csv'))DT::datatable(ico_assess_raw, caption = 'ICO species and past IUCN assessments')These raw assessments are cleaned up to standardize codes and categories, and category scores are assigned based on a scale from Least Concern = 1 to Extinct = 0, as per the Species OHI subgoal. Note that past assessments used different coding, especially prior to 1994; we reclassified older codes according to this chart:
iucn_sid | year | code | category | sciname
ico_assess_raw <- read_csv(file.path(dir_goal, 'raw/ico_assessments_raw.csv'))
ico_assess <- ico_assess_raw %>%
rename(cat = code, cat_txt = category) %>%
mutate(cat = toupper(cat),
cat = str_replace(cat, 'LR/', ''),
cat = ifelse(cat %in% c('K', 'I'), 'DD', cat),
cat = ifelse(cat == 'NR', 'NE', cat),
cat = ifelse(str_detect(toupper(cat_txt), 'VERY RARE'), 'CR', cat),
cat = ifelse(str_detect(toupper(cat_txt), 'LESS RARE'), 'T', cat),
cat = ifelse(str_detect(toupper(cat_txt), 'STATUS INADEQUATELY KNOWN'), 'DD', cat),
cat = ifelse(cat == 'V', 'VU', cat),
cat = ifelse(cat == 'E', 'EN', cat))
pop_cat <- data.frame(cat = c("LC", "NT", "VU", "EN", "CR", "EX", "T", "CD", "NE", "DD"),
cat_score = c( 0, 0.2, 0.4, 0.6, 0.8, 1.0, 0.6, 0.3, NA, NA),
stringsAsFactors = FALSE)
ico_assess <- ico_assess %>%
left_join(pop_cat, by = 'cat') %>%
filter(!is.na(cat_score)) %>%
distinct() %>%
arrange(iucn_sid, year)
write_csv(ico_assess, file.path(dir_goal, 'int/ico_assess_clean.csv'))Using tidyr::complete() and tidyr::fill(), we create a full time series for all species from the earliest assessment to the most recent year.
ico_assess <- read_csv(file.path(dir_goal, 'int/ico_assess_clean.csv'))
ico_list <- read_csv(file.path(dir_goal, 'int/ico_list_prepped.csv'))
ico_assess_full <- ico_assess %>%
mutate(eval_yr = year) %>%
select(-sciname) %>%
arrange(iucn_sid, year) %>%
complete(year = full_seq(year, 1), nesting(iucn_sid)) %>%
group_by(iucn_sid) %>%
fill(cat, cat_txt, cat_score, eval_yr) %>% ## fills all the way to latest year
ungroup()
ico_spp_cat <- ico_list %>%
rename(cat_2016 = cat) %>%
left_join(ico_assess_full, by = c('iucn_sid'))
## if no time series available, time series years will be NA. Assign a list to
## those NAs, then unnest it to create observations for those years.
ico_spp_cat <- ico_spp_cat %>%
mutate(year = ifelse(is.na(year),
list(c(min(year, na.rm = TRUE):max(year, na.rm = TRUE))),
year)) %>%
unnest(year)
## NAs will be filled backward in time by starting from the most recent non-NA.
## To do this, we'll swap any current-year NAs with the cat_score (meaning no
## time series fill), and fill upwards instead of downwards.
ico_spp_cat <- ico_spp_cat %>%
left_join(pop_cat %>%
rename(cat_2016 = cat, cat_2016_score = cat_score),
by = 'cat_2016') %>%
mutate(cat_score = ifelse(year == max(year, na.rm = TRUE) & is.na(cat),
cat_2016_score,
cat_score)) %>%
arrange(iucn_sid, year) %>%
group_by(iucn_sid) %>%
fill(cat, cat_score, cat_txt, .direction = 'up') %>%
ungroup() %>%
distinct()
write_csv(ico_spp_cat, file.path(dir_goal, 'int/ico_spp_cat.csv'))Using dplyr::full_join() we combine the ico_spp_rgn dataframe (iconic species by OHI region) with the ico_spp_cat dataframe (iconic species by category and year, with species info, year, and category info).
ico_cat_ts_abbr <- read_csv(file.path(dir_goal, 'int/ico_spp_cat.csv')) %>%
select(iucn_sid, sciname, year, cat, cat_score, eval_yr) %>%
filter(year >= 1992)
ico_spp_rgn <- read_csv(file.path(dir_goal, 'int/ico_spp_rgn_prepped.csv')) %>%
select(rgn_id, rgn_name, iucn_sid, comname, sciname, ico_gl, ico_rgn_id, presence)
ico_spp_rgn_cat <- ico_cat_ts_abbr %>%
full_join(ico_spp_rgn, by = c('iucn_sid', 'sciname'))
### No information on when species go extinct locally, so just set all years to extinct for that region
ico_spp_rgn_cat <- ico_spp_rgn_cat %>%
mutate(cat = ifelse(str_detect(presence, '^Extinct'), 'EX', cat), # ^ indicates start of string
cat_score = ifelse(cat == 'EX', 1, cat_score)) %>%
filter(ico_gl | ico_rgn_id == rgn_id) %>% # keep (all globally iconic) and (regionally iconic in region only)
distinct()
write_csv(ico_spp_rgn_cat, file.path(dir_goal, 'int/ico_spp_rgn_cat.csv'))knitr::kable(head(ico_spp_rgn_cat, 10))The toolbox wants rgn_id, species sciname, and extinction risk category for the basic calculations. Since some regions contain multiple subpops (or parent/subpop) we also include iucn_sid to differentiate. This information is included for each year, filtered back to the year 2000.
While the official calculations are performed in the toolbox, we perform the same basic calcs here to get a sense of the ICO status and trend ahead of time. Report and summarize estimate of regional iconic species status:
ico_spp_rgn_cat <- read_csv(file.path(dir_goal, 'int/ico_spp_rgn_cat.csv'))
## Report out for toolbox format (rgn_id | sciname | category or popn_trend for each species within a region).
## Note: in toolbox, group_by(rgn_id, sciname) and then summarize(category = mean(category)) to
## average any parent/subpop species listings before aggregating to overall average per region.
ico_status_raw <- ico_spp_rgn_cat %>%
select(rgn_id, rgn_name, sciname, iucn_sid, cat, cat_score, year, eval_yr) %>%
arrange(rgn_id, desc(year), sciname) %>%
ungroup()
ico_status_calc <- ico_status_raw %>%
group_by(rgn_id, rgn_name, sciname, year) %>%
filter(!is.na(cat_score)) %>% # remove any DDs
summarize(cat_score = mean(cat_score)) %>%
group_by(rgn_id, rgn_name, year) %>%
summarize(mean_cat = round(mean(cat_score), 5),
ico_status = (1 - mean_cat) * 100,
n_spp = n()) %>%
ungroup()
ico_trend <- data.frame()
for (i in 1993:max(ico_status_calc$year, na.rm = TRUE)) { # i <- 2013
tmp_status <- ico_status_calc %>%
filter(year <= i & year > (i - 10)) # trend based on 10-year average since assessments are sporadic
tmp_trend <- tmp_status %>%
group_by(rgn_id) %>%
do(trend_lm = lm(ico_status ~ year, data = .)$coefficients[2]) %>%
mutate(year = i,
trend_lm = as.numeric(trend_lm)/100, # status is 0 - 100; trend should be +1 to -1
ico_trend = round(trend_lm * 5, 5)) %>% # trend prediction five years out
ungroup()
ico_trend <- ico_trend %>%
bind_rows(tmp_trend)
}
ico_sum <- ico_status_raw %>%
left_join(ico_status_calc, by = c('rgn_id', 'rgn_name', 'year')) %>%
left_join(ico_trend, by = c('rgn_id', 'year'))
write_csv(ico_sum, file.path(dir_goal, 'summary/ico_summary.csv'))
### Report out for finalized status and trend values per region
ico_status_raw1 <- ico_status_raw %>%
dplyr::select(rgn_id, sciname, iucn_sid, year, eval_yr, category = cat)
write_csv(ico_status_raw1, file.path(dir_goal, 'output/ico_spp_iucn_status.csv'))
write_csv(ico_status_calc, file.path(dir_goal, 'output/ico_status_calc.csv'))
write_csv(ico_trend, file.path(dir_goal, 'output/ico_trend.csv'))
ico_status_raw1[duplicated(ico_status_raw1 ), ]
### NOTE: if iucn_sid were removed, this would show duplicates due to subpops with same category.
table(ico_status_raw1$category)DT::datatable(ico_status_raw %>% filter(year == 2017))DT::datatable(ico_status_calc %>% filter(year == 2017), caption = 'ICO status')
DT::datatable(ico_trend %>% filter(year == 2017) %>% select(-trend_lm), caption = 'ICO trend')To examine results of the new methods (including API-based data scraping and trend based on time series), we plot the estimated status and trend over time.
library(ggplot2)
library(plotly)
status_ts_plot <- ggplot(ico_sum %>%
filter(!is.na(rgn_id)),
aes(x = year, y = ico_status, color = rgn_id, group = rgn_id)) +
geom_line(size = 1.2, alpha = .4) +
labs(x = 'year',
y = 'ICO status',
title = 'ICO status over time',
color = 'Region')
ggplotly(status_ts_plot)
trend_ts_plot <- ggplot(ico_sum %>%
filter(!is.na(rgn_id) &!is.na(ico_trend)),
aes(x = year, y = ico_trend, color = rgn_id, group = rgn_id)) +
geom_line(size = 1.2, alpha = .4) +
labs(x = 'year',
y = 'ICO trend',
title = 'ICO trend over time',
color = 'Region')
ggplotly(trend_ts_plot)Plot the estimated status and trends for 2016 from assessments v2017 and v2018.
status_trend <- ico_sum %>% # status and trend in this year's assessment
filter(year == max(year)-2) %>% # max year of ico_sum is 2018 not 2017
select(rgn_id, ico_status, ico_trend, n_spp) %>%
distinct() # because removing iucn_sid, cat, cat_score, and species name, results in duplicates
ico_compare <- read_csv(file.path(dir_previous, 'summary/ico_summary.csv')) %>%
filter(year == max(year)) %>% # max year of v2017 ico_summary is 2016; make sure to update for v2019
select(rgn_id, rgn_name, ico_status, ico_trend, n_spp) %>%
distinct() %>%
rename(previous_status = ico_status, previous_trend = ico_trend, previous_nspp = n_spp) %>%
left_join(status_trend, by = 'rgn_id') %>%
mutate(d_nspp = as.factor(previous_nspp - n_spp))
compare_status_plot <- ggplot(ico_compare) +
geom_point(aes(previous_status, ico_status, color = d_nspp, label = rgn_id), alpha = 0.6) +
geom_abline(slope = 1, intercept = 0, color = 'red') +
labs(x = sprintf('ICO status v%s', previous_scen),
y = sprintf('ICO status %s', scenario),
title = 'ICO Status Comparison')
ggplotly(compare_status_plot)
trend_compare_plot <- ggplot(ico_compare) +
geom_point(aes(previous_trend, ico_trend, color = d_nspp, label = rgn_id), alpha = .6) +
geom_abline(slope = 1, intercept = 0, color = 'red') +
labs(x = sprintf('ICO trend v%s', previous_scen),
y = sprintf('ICO trend %s', scenario),
title = 'ICO Trend Comparison')
ggplotly(trend_compare_plot)