README.Rmd

---
title: "Habits"
author: Gilles Dauby
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: md_document
---


This document document R scripts and functions used for the analysis presented in the manuscript entitled __A third of the tropical African flora is potentially threatened with extinction__ published in Sciences Advances (Stévart _et al._ 2018)



### Load packages (to be installed if needed)

```{r, include=TRUE, echo=TRUE, message=F, warning=F}
library(tidyverse)
library(sf)
library(raster)
library(ConR)
library(broom)
library(rredlist)
library(rgbif)
library(doParallel)
library(maps)
library(tmap)
```


### load dataset

```{r, include=TRUE, echo=TRUE, message=F, warning=F}
dataset <- read_csv("data/dataset_rainbio_used.csv")
```


```{r, include=TRUE, echo=TRUE, message=F, warning=F}
dataset
```


### extracting data from gbif for all species 
```{r, include=TRUE, echo=TRUE, message=F, warning=F, eval=F}


species_list <-
  dataset %>%
  distinct(tax_sp_level)

list_data <- as.list(pull(species_list))


#### Searching for gbif occurrences for all species, takes around 5 hours in parallel on 4 cores.
doParallel::registerDoParallel(4) ## define here the number of core on which you want to run the search

system.time(results <-
  plyr::llply(list_data, .fun=function(x) {
    # source("functions_criteria_A.R")
    source("gbif_search_filtered_function.R")
  
    gbif_search_filtered(taxa = x)
  }
    , .progress = "text", .parallel=T, .paropts=list(.packages=c('rgbif','raster','tidyverse', 'ConR', 'sf'))))


doParallel::stopImplicitCluster()

gbif_data <- bind_rows(results)

# write_csv(gbif_data, "gbif.searched.all.taxa.csv")

```

```{r, include=F, echo=F, message=F, warning=F}
gbif_data <- read_csv("data/gbif.searched.all.taxa.csv")
```


This table provide for each species the total number of occurence in gbif, the total number of occupied cells in 10 km cell size grid and the list of Continent where the species is recorded on gbif.

```{r, include=T, echo=F, message=F, warning=F}
gbif_data
```

### load SIG dataset
```{r, include=T, echo=T, message=F, warning=F}
## protected areas network
protected_areas_network <- sf::read_sf("data/NationalParks_filtered_corrected.shp")
## mineral deposits plus a buffer of 10 km around each point
mineral_deposits <- sf::read_sf("data/poly_mineral_deposit.shp")
raster_mayaux <- raster("data/sum.rasters.human.impacted.tif")

### getting a nice and updated world map
world <- maps::map("world", fill=TRUE, plot=FALSE)
IDs <- sapply(strsplit(world$names, ":"), function(x) x[1])
world <- maptools::map2SpatialPolygons(world, IDs=IDs, proj4string=sp::CRS("+proj=longlat +datum=WGS84"))
suppressWarnings(world_map <- broom::tidy(world))
world_id <- tibble(id_nbe=1:nrow(distinct(world_map, id)), id=pull(distinct(world_map, id)))

world_map <- as(world, "sf")
world_map <- world_map %>% 
  add_column(id_nbe=world_id$id_nbe, id=world_id$id)
```



### Preliminary Automated Conservation Assessment following Criterion B

```{r CriterionB-Analysis, eval=F, echo=T, warning=F, message=F}

protected_areas_network_sp <- 
  as(protected_areas_network, "Spatial")

criterion_b <-
  IUCN.eval(dplyr::select(dataset, ddlat, ddlon, tax_sp_level)[1:1000,], 
            protec.areas = protected_areas_network_sp, DrawMap = F, SubPop = F, file_name = "dataset_criterion_B", write_file_option = "csv", parallel = T, NbeCores = 4)


```

```{r, include=F, echo=F, message=F, warning=F}
criterion_b <- read_csv("data/dataset_Criterion_B.csv")
```


```{r, include=T, echo=F, message=F, warning=F}
criterion_b
```



### Preliminary Automated Conservation Assessment following Criterion A


```{r CriterionA-Analysis, eval=F, echo=TRUE, warning=F, message=F}

source("criteria_A_function.R") 


list_data <- split(dplyr::select(dataset, ddlat, ddlon, tax_sp_level), 
                   f =dataset$tax_sp_level)

### Take about 1 hour

doParallel::registerDoParallel(4)

system.time(Results <-
  plyr::llply(list_data, .fun=function(x) {

    source("criteria_A_function.R")
    sum.rasters.human.impacted <- raster("data/sum.rasters.human.impacted.tif")
    mineral_deposits <- sf::read_sf("data/poly_mineral_deposit.shp")
    protected_areas_network <- sf::read_sf("data/NationalParks_filtered_corrected.shp")
    
    IUCN_eval_CA(data = x, rasters = sum.rasters.human.impacted, mineral = mineral_deposits, protected_areas = protected_areas_network)
    
  }
    , .progress = "text", .parallel=T, .paropts=list(.packages=c('raster','tidyverse', 'sf'))))

criterion_a<- bind_rows(Results)

criterion_a <-
  criterion_a %>%
  mutate(AOO_without_impact=AOO_all-AOO_left) %>%
  mutate(AOO_decline=(AOO_all-AOO_without_impact)/AOO_all*100)

criterion_a <- criterion_a %>%
  add_column(Category_code_CA=plyr::aaply(criterion_a$AOO_decline, 1, Cat_Crition_A))

write_csv(criterion_a, "dataset_Criterion_A.csv")


```



```{r, include=F, echo=F, message=F, warning=F}
criterion_a <- read_csv("data/dataset_Criterion_A.csv")
```


```{r, include=T, echo=F, message=F, warning=F}
criterion_a
```

### Obtain published full IUCN assessment from the Red List

You will need a token from the Red List to use the API. Request a token [HERE](http://apiv3.iucnredlist.org/api/v3/token)


```{r, eval=F, echo=TRUE, warning=F, message=F}

species_list <- distinct(dataset, tax_sp_level)

# i <- 1
all.cons.status <- list()

for (i in 1:nrow(species_list)) { # 
  cat(i, " ")
res_iucn <-
  rredlist::rl_search(name = species_list$tax_sp_level[i], parse =F, 
            key = "YourTokenID") ### replace your own Token here

  if(length(res_iucn$result)>0) {
    
    full <- lapply(res_iucn$result[[1]], function(x) ifelse(is.null(x), NA, x))
    full$tax_sp_level <- species_list$tax_sp_level[i]
      
    all.cons.status[[length(all.cons.status)+1]] <- do.call("cbind", full) %>%
      as_tibble()
    
    # print(all.cons.status[[length(all.cons.status)]])
  }
}

IUCN.cons.status.tb <- bind_rows(all.cons.status)
```

```{r, include=F, echo=F, message=F, warning=F}
IUCN.cons.status.tb <- read_csv("data/IUCN.cons.status.tb.final.csv")
```

```{r, include=T, echo=F, message=F, warning=F}
IUCN.cons.status.tb
```




```{r eval=F, echo=F}
### Intersecting potentially/likely threatened species with phytographical regions
dataset <- read_csv("data/dataset_rainbio_used.csv")
dataset_ca <- read_csv("data/dataset_Criterion_A.csv")
dataset_cb <- read_csv("data/dataset_Criterion_B.csv")

dataset_ca <- dataset_ca %>%
  mutate(code_ca = dplyr::case_when(
    Category_code_CA == "LC or NT" ~ 4,
    Category_code_CA == "CR" ~ 1,
    Category_code_CA == "VU" ~ 3,
    Category_code_CA == "EN" ~ 2
  ))

dataset_cb <- dataset_cb %>%
  mutate(code_cb = dplyr::case_when(
    Category_CriteriaB == "LC or NT" ~ 4,
    Category_CriteriaB == "CR" ~ 1,
    Category_CriteriaB == "VU" ~ 3,
    Category_CriteriaB == "EN" ~ 2
  ))


phyto_regions <- st_read("data/RAINBIO_phytogeography.shp")

dataset <- dataset %>%
  left_join(dataset_ca %>%
              dplyr::select(tax_sp_level, code_ca)) %>%
  left_join(dataset_cb %>%
              dplyr::select(tax_sp_level, code_cb))

dataset_sf <- st_as_sf(dataset, coords = c("ddlon", "ddlat"))
st_crs(dataset_sf) <- 4326

dataset_sf_inter <- st_intersection(dataset_sf, phyto_regions)

dataset_phytoreg <- 
  dataset_sf_inter %>%
  st_set_geometry(NULL) %>% 
  mutate(code_c_min = pmin(code_ca, code_cb))

nbr_species_phytoreg <- 
  dataset_phytoreg %>%
  distinct(Bioreg_ID, tax_sp_level, code_c_min) %>%
  pivot_wider(
    names_from = c("code_c_min"),
    values_from = "tax_sp_level",
    values_fn = ~ length(unique(.x)),
    names_prefix = "status_"
  )

nbr_species_phytoreg <- 
  nbr_species_phytoreg %>% 
  mutate(nbr_sp = status_1 + status_2 + status_3 + status_4) %>% 
  mutate(prop_threat = (status_1 + status_2 + status_3)/nbr_sp*100,
         prop_lik_threat = (status_1 + status_2)/nbr_sp*100)

dataset_phytoreg <- 
  dataset_sf_inter %>%
  st_set_geometry(NULL) %>% 
  mutate(code_c_min = pmin(code_ca, code_cb))

nbr_species_phytoreg_cb <- 
  dataset_phytoreg %>%
  distinct(Bioreg_ID, tax_sp_level, code_cb) %>%
  pivot_wider(
    names_from = c("code_cb"),
    values_from = "tax_sp_level",
    values_fn = ~ length(unique(.x)),
    names_prefix = "status_"
  )

nbr_species_phytoreg_cb <- 
  nbr_species_phytoreg_cb %>% 
  mutate(nbr_sp = status_1 + status_2 + status_3 + status_4) %>% 
  mutate(prop_threat_cb = (status_1 + status_2 + status_3)/nbr_sp*100,
         prop_lik_threat_cb = (status_1 + status_2)/nbr_sp*100)

nbr_species_phytoreg_ca <- 
  dataset_phytoreg %>%
  distinct(Bioreg_ID, tax_sp_level, code_ca) %>%
  pivot_wider(
    names_from = c("code_ca"),
    values_from = "tax_sp_level",
    values_fn = ~ length(unique(.x)),
    names_prefix = "status_"
  )

nbr_species_phytoreg_ca <- 
  nbr_species_phytoreg_ca %>% 
  mutate(nbr_sp = status_1 + status_2 + status_3 + status_4) %>% 
  mutate(prop_threat_ca = (status_1 + status_2 + status_3)/nbr_sp*100,
         prop_lik_threat_ca = (status_1 + status_2)/nbr_sp*100)

phyto_regions <- phyto_regions %>% 
  left_join(nbr_species_phytoreg %>% dplyr::select(Bioreg_ID, prop_threat, prop_lik_threat)) %>% 
  left_join(nbr_species_phytoreg_cb %>% dplyr::select(Bioreg_ID, prop_threat_cb, prop_lik_threat_cb)) %>% 
  left_join(nbr_species_phytoreg_ca %>% dplyr::select(Bioreg_ID, prop_threat_ca, prop_lik_threat_ca))

st_tr_map <- tm_shape(phyto_regions) + 
  tm_rgb() +
  tm_shape(all_waypoints_tr_sf) +
  tm_dots(size = 1, alpha = 0.6, shape = 19, col = "blue") + 
  tm_compass(type = "8star", position = c("left", "top")) +
  tm_scale_bar(breaks = c(0, 5, 10, 15), text.size = 1, position = "RIGHT") + 
  tm_layout(frame = T)

data("World")
africa <- 
  World %>% 
  filter(continent == "Africa", name != "Madagascar")

cb_prop <- tm_shape(africa) +
  tm_polygons("continent",
              legend.show = F,
              palette = "grey") +
  tm_shape(phyto_regions %>% filter(Bioreg_ID != "T")) +
  tm_polygons("prop_threat_cb",
              title = "Proportion \nof threatened \nspecies following \ncriterion B (%)") +
  tm_layout(frame = T) +
  tm_legend(legend.title.size = 2.,
            legend.text.size = 1.5)

tmap_save(cb_prop, "cb_prop_phytoreg.png")

ca_prop <- tm_shape(africa) +
  tm_polygons("continent",
              legend.show = F,
              palette = "grey") +
  tm_shape(phyto_regions %>% filter(Bioreg_ID != "T", !is.na(prop_threat_ca))) +
  tm_polygons("prop_threat_ca",
              title = "Proportion \nof threatened \nspecies following \ncriterion A (%)") +
  tm_layout(frame = T) +
  tm_legend(legend.title.size = 2.,
            legend.text.size = 1.5)

tmap_save(ca_prop, "ca_prop_phytoreg.png")

cb_prop <- tm_shape(africa) +
  tm_polygons("continent",
              legend.show = F,
              palette = "grey") +
  tm_shape(phyto_regions %>% filter(Bioreg_ID != "T")) +
  tm_polygons("prop_threat",
              title = "Proportion \nof threatened \nspecies (%)") +
  tm_layout(frame = T) +
  tm_legend(legend.title.size = 2.,
            legend.text.size = 1.5)

tmap_save(cb_prop, "prop_phytoreg.png")


```