4a_map_impacts_per_taxon.Rmd

---
title: 'Map impacts per taxon'
author: "*Compiled on `r date()` by `r Sys.info()['user']`*"
output: 
  html_document:
    code_folding: hide
    toc: true
    toc_depth: 3
    toc_float: yes
    number_sections: true
    theme: cerulean
    highlight: haddock
    includes: 
      in_header: '~/github/src/templates/ohara_hdr.html'
  pdf_document:
    toc: true
---

``` {r setup, echo = TRUE, message = FALSE, warning = FALSE}

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

library(raster)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')

source(here('common_fxns.R'))

```


# Summary

Create taxa-level maps of number of species impacted per cell, per year.

# Methods

* Gather species maps by taxonomic group (IUCN comprehensively assessed groups)
* read in all impact maps for species for all years, flatten to impacted/not for each year
* add maps together for a taxa-level map for each year of number of impacted species per cell
* repeat similar process for refugia.

## Match threatened species to taxonomic groups per IUCN

Identify the species included in the impacted list (using the IUCN species IDs of impact map files, from prior scripts!).  Join this with the list of mapped species, to identify which taxonomic group the species is in (based on the downloaded IUCN shapefiles).  Taxa groups not represented in the impacted species list will be dropped.

```{r identify taxa groups}

impact_map_dir <- file.path(dir_bd_anx, 'spp_impact_rasts')

### identify spp with impact maps
spp_impacted <- list.files(impact_map_dir, recursive = TRUE) %>%
  str_extract('[0-9]+') %>%
  as.integer() %>%
  unique()

### Get inclusion list; drop species with no sensitivities to stressors.
### Also append the class level taxa instead of assessment group taxa
spp_sens <- get_incl_spp() %>%
  mutate(taxon = str_replace_all(desc, ' ', ''))

### look for dropped taxonomic groups...
comp_file <- here('_data', sprintf('iucn_comp_assessed_%s.csv', api_version))
spp_comp <- read_csv(comp_file)

taxa_dropped <- spp_comp %>%
  filter(!assess_gp %in% spp_sens$assess_gp) %>%
  .$assess_gp %>% unique()

```

Taxonomic groups dropped from this analysis - i.e., group does not have any species that fits the "impacted" criteria:

`r taxa_dropped`

The one vulnerable tarpon species (i.e. not LC, EX, or DD) is not sensitive to any of the stressors and is not mapped, so the whole taxon is dropped from this analysis.

``` {r setup directories and parameters}
### directories for inputs
spp_impact_map_dir   <- file.path(dir_bd_anx, 'spp_impact_rasts')
spp_range_dir        <- file.path(dir_bd_anx, 'spp_rasts_mol_2020')

### directories for outputs
taxon_impact_map_dir <- file.path(dir_bd_anx, 'taxon_impact_rasts')
### unlink(file.path(taxon_impact_map_dir, 'taxon_impact_birds_*.*'))
### unlink(file.path(taxon_impact_map_dir, '*.*'))

cell_id_rast <- raster(here('_spatial/cell_id_mol.tif'))
ocean_a_rast <- raster(here('_spatial/ocean_area_mol.tif'))

```

## Loop over taxonomic groups to create taxa-level impact maps

These maps ignore the number of impacts occurring on each species, and count up the number of species impacted by the aggregated stressor group (land-based, ocean, climate, and fishing).

We will calculate unweighted impacts and range-weighted impacts based on protected area inclusion targets from Butchart et al 2015, see function `get_spp_priority()` for details.

For a given taxon:

* identify full list of impacted species
* select a subgroup (~10 spp) to process at a time
    * read in all spp in group
    * flatten all spp impact maps (for impacts > 0, set "impacted" = TRUE)
        * no need to worry about refugia cells at this stage.
    * add all maps for this group: by cell/year per stressor
    * also add priority values for each cell/year per stressor
* for all subgroups in a taxa, add them up.
* Repeat this but focusing only on cells where species are impacted by multiple stressors.
    
```{r set up functions}

get_spp_vec <- function(spp_ids, i, subgp_size) {
  gp_first <- (i - 1) * subgp_size + 1
  gp_last  <- min(i * subgp_size, length(spp_ids))
  spp_id_vec <- spp_ids[gp_first:gp_last]
}

get_spp_impact <- function(spp_id, imp_file, n_imp = 1) {
  ### n_imp is the threshold for minimum number of impacts to include
  ### e.g. n_imp == 2 counts all cells with at least two stressors
  impact_map <- read_csv(imp_file, col_types = 'iii') %>%
    filter(!is.na(year) & !is.na(cell_id)) %>%
    filter(n_impacts >= n_imp & !is.na(n_impacts)) %>%
    select(cell_id, year) %>%
    mutate(iucn_sid = spp_id)
  return(impact_map)
}

taxon_gps <- spp_sens$taxon %>% unique() %>% sort()

subgp_size <- 24 ### number of spp to process at a time
year_span <- 2003:2013
reload <- FALSE
```

```{r create taxa maps}
str_cats <- c('land-based', 'ocean', 'fishing', 'climate', 'all')

for(str_cat in str_cats) {
  # str_cat <- str_cats[5]
  
  for(taxon_gp in taxon_gps) {
    ### taxon_gp <- taxon_gps[6]
    
    taxon_map_files <- file.path(taxon_impact_map_dir, 
                                sprintf('taxon_impact_%s_%s_%s.tif',
                                        taxon_gp, year_span, str_cat))

    if(any(!file.exists(taxon_map_files)) | reload) {
      spp_ids_taxa <- spp_sens %>%
        filter(taxon == taxon_gp) %>%
        filter(str_cat == 'all' | category == str_cat) %>%
        .$iucn_sid %>% 
        unique() %>% sort()
      
      n_subgps <- ceiling(length(spp_ids_taxa) / subgp_size)
      message('Processing ', str_cat, ' impacts for ', taxon_gp, 
              ': ', length(spp_ids_taxa), 
              ' species broken into ', n_subgps, ' groups.')
      
      if(length(spp_ids_taxa) == 0) {
        ### uh oh, no species in this taxon affected by this stressor.
        taxon_impacts_df <- data.frame()
      } else {
        ### loop over all spp in this taxon affected by this stressor
        ### initialize a list to store 
        taxon_impacts_list <- vector('list', length = n_subgps)
        for(i in 1:n_subgps) {
          # i <- 1
          message('  processing group ', i)
          spp_id_vec <- get_spp_vec(spp_ids_taxa, i, subgp_size)
          
          impact_maps <- parallel::mclapply(spp_id_vec,
              FUN = function(spp_id) {
                ### spp_id <- spp_id_vec[1]
                imp_file <- file.path(spp_impact_map_dir, 
                                   sprintf('spp_impact_map_%s_%s.csv', spp_id, str_cat))
                imp_map <- get_spp_impact(spp_id, imp_file)
              }, mc.cores = subgp_size) %>%
            bind_rows()
          
          message('  summarizing group ', i)
          # system.time({
          impacts_gp <- impact_maps %>%
            group_by(cell_id, year) %>%
            summarize(n_spp = n()) %>%
            ungroup()
          # })
          
          taxon_impacts_list[[i]] <- impacts_gp
        } ### end of for loop
      
      taxon_impacts_df <- bind_rows(taxon_impacts_list) %>%
        group_by(cell_id, year) %>%
        summarize(n_spp = sum(n_spp, na.rm = TRUE)) %>%
        ungroup() %>%
        mutate(taxon = taxon_gp)
      
      } ### end of if statement for zero-length species vector
      
      if(nrow(taxon_impacts_df) == 0) {
        ### either no species in the list, or no impacts for any of the spp
        taxon_impacts_df <- data.frame(cell_id = NA,
                                       year    = NA,
                                       n_spp   = NA,
                                       taxon   = taxon_gp)
      }
      
      ### write out as count by year
      for(i in seq_along(year_span)) {
        # i <- 1
        taxon_nspp_map_file <- taxon_map_files[i]
        
        taxon_impact_yr_df <- taxon_impacts_df %>%
          filter(year == year_span[i])
        
        nspp_yr_rast <- map_to_rast(taxon_impact_yr_df, 
                                    cell_val = 'n_spp')
        writeRaster(nspp_yr_rast, taxon_nspp_map_file, overwrite = TRUE)
      }
      
    } else {
      # message('Files exist... skipping!')
    }
  }
}
```

```{r create taxa maps with multiple stressors - all 2013 only}

for(taxon_gp in taxon_gps) {
  ### taxon_gp <- taxon_gps[6]
  
  taxon_2plus_map_file <- file.path(taxon_impact_map_dir, 
                              sprintf('taxon_impact_%s_2013_all_2plus.tif',
                                      taxon_gp))

  if(any(!file.exists(taxon_2plus_map_file)) | reload) {
    spp_ids_taxa <- spp_sens %>%
      filter(taxon == taxon_gp) %>%
      filter(str_cat == 'all' | category == str_cat) %>%
      .$iucn_sid %>% 
      unique() %>% sort()
    
    n_subgps <- ceiling(length(spp_ids_taxa) / subgp_size)
    message('Processing ', str_cat, ' impacts for ', taxon_gp, 
            ': ', length(spp_ids_taxa), 
            ' species broken into ', n_subgps, ' groups.')
    
    if(length(spp_ids_taxa) == 0) {
      ### uh oh, no species in this taxon affected by this stressor.
      taxon_impacts_df <- data.frame()
    } else {
      ### loop over all spp in this taxon affected by this stressor
      ### initialize a list to store 
      taxon_impacts_list <- vector('list', length = n_subgps)
      for(i in 1:n_subgps) {
        # i <- 1
        message('  processing group ', i)
        spp_id_vec <- get_spp_vec(spp_ids_taxa, i, subgp_size)
        
        impact_maps <- parallel::mclapply(spp_id_vec,
            FUN = function(spp_id) {
              ### spp_id <- spp_id_vec[1]
              imp_file <- file.path(spp_impact_map_dir, 
                                 sprintf('spp_impact_map_%s_all.csv', spp_id))
              imp_map <- get_spp_impact(spp_id, imp_file, n_imp = 2)
            }, mc.cores = subgp_size) %>%
          bind_rows()
        
        message('  summarizing group ', i)
        # system.time({
        impacts_gp <- impact_maps %>%
          filter(year == 2013) %>%
          group_by(cell_id) %>%
          summarize(n_spp = n()) %>%
          ungroup()
        # })
        
        taxon_impacts_list[[i]] <- impacts_gp
      } ### end of for loop
    
    taxon_impacts_df <- bind_rows(taxon_impacts_list) %>%
      group_by(cell_id) %>% ### no need to group by year
      summarize(n_spp = sum(n_spp, na.rm = TRUE)) %>%
      ungroup() %>%
      mutate(taxon = taxon_gp)
    
    } ### end of if statement for zero-length species vector
    
    if(nrow(taxon_impacts_df) == 0) {
      ### either no species in the list, or no impacts for any of the spp
      taxon_impacts_df <- data.frame(cell_id = NA,
                                     n_spp   = NA,
                                     taxon   = taxon_gp)
    }
    
    nspp_imp2plus_rast <- map_to_rast(taxon_impacts_df, 
                                 cell_val = 'n_spp')
    writeRaster(nspp_imp2plus_rast, taxon_2plus_map_file, overwrite = TRUE)

  } else {
    # message('Files exist... skipping!')
  }
}
```

## Animate select taxa impact maps 

```{r animate function, eval = FALSE}

library(animation)

make_ramp <- function(n, palette, log_out = TRUE) {
  if(log_out) {
    n_even <- ceiling(n)
    n_max <- 10^n_even
    log_steps <- c(1, 2, 5, 10, 20, 50, 100, 200, 500, 1000)
    log_steps <- log_steps[log_steps <= n_max]
    
    breaks <- seq(0, n_even, length.out = 100)
    lblpos <- log10(log_steps)
    labels <- 10^lblpos
    colors <- hcl.colors(length(breaks), palette = palette)
  } else {
    n_even <- round(n + 50, -2)
    breaks <- c(0, seq(.5, n_even, length.out = 100))
    lblpos <- seq(0, n_even, by = 50)
    labels <- lblpos
    colors <- c('grey40', hcl.colors(100, palette = palette))
  }
  
  ### For breaks and colors, add a different color for 
  ### a value of exactly zero
  return(list('colors' = colors,
              'breaks' = breaks,
              'labels' = labels,
              'lblpos' = lblpos))
}

make_gifs <- function(map_stack, filename, layer_names = NULL, log_out = TRUE) {
  if(is.null(layer_names)) layer_names = names(map_stack)

  if(log_out) map_stack <- log10(map_stack)
  
  n <- max(c(maxValue(map_stack), 1), na.rm = TRUE)

  ramp <- make_ramp(n, palette = 'viridis', log_out)
  
  capture.output({
    saveGIF({
      for(i in 1:nlayers(map_stack)){
        plot(map_stack[[i]], 
             col = ramp$colors,
             breaks = ramp$breaks,
             axes = FALSE,
             axis.args = list(at = ramp$lblpos, labels = ramp$labels),
             main = layer_names[i])
      }}, 
      interval = 0.5, movie.name = filename, 
      ani.width = 700, ani.height = 420)
  })
  
  return(invisible(NULL))
}

```

```{r impact maps by taxa, eval = FALSE}

mapfiles <- list.files(taxon_impact_map_dir, full.names = TRUE)

all_taxa_maps <- data.frame(f = mapfiles) %>%
  mutate(impact = str_replace_all(basename(f), 'taxon.+_[0-9]{4}_|.tif', ''),
         type = ifelse(str_detect(basename(f), 'priority'), 'priority', 'sum'),
         year = str_extract(basename(f), '[0-9]{4}') %>% as.integer(),
         taxon = str_replace_all(basename(f), 'taxon.+impact_|_[0-9].+', ''))

str_cats <- c('fishing', 'land-based', 'ocean', 'climate', 'all')
taxa <- all_taxa_maps$taxon %>% 
  unique()
taxa <- taxa[c(1, 2, 3, 5, 9)]

```

``` {r animate impacts, results = 'asis', eval = FALSE}
### Animate the results

for(str_cat in str_cats) {
  # str_cat <- str_cats[3]
  for(t in taxa) {
    # t <- taxa[2]
    gif_file <- here(sprintf('figs/impact_%s_%s.gif', t, str_cat))
    
    if(!file.exists(gif_file) | reload) {
      taxon_map_df <- all_taxa_maps %>%
        filter(taxon == t) %>%
        filter(impact == str_cat) %>%
        filter(type == 'sum')
      rast_files <- taxon_map_df$f
      map_stack <- stack(rast_files)
      
      make_gifs(map_stack, 
                filename = gif_file,
                layer_names = paste(t, str_cat, 'impacts', taxon_map_df$year))
    }
    
    cat(sprintf('![](%s)', gif_file))
  }
}
```