4b_map_trends_per_taxon.Rmd

---
title: 'Map trends 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 with intensifying or deintensifying impacts per cell, per year.

# Methods

* Gather species maps by taxonomic group (IUCN comprehensively assessed groups)
* read in all trend maps for all stressors, flatten to increasing/decreasing/indeterminate for each year (with different thresholds for incr/decr)
* add maps together for a taxa-level map of number of increasingly/decreasingly impacted species per cell

Individual stressor trend layers are clipped to their respective footprints based on the 95% contour volume applied to the 2013 stressor map.  Therefore, cells with a significant trend but no footprint are dropped.  In physical terms, this is a forward-looking approach, and means one of:

* cells with a stressor with decreasing trend, and below the threshold, have little room to continue decreasing - effectively the stressor is eliminated and can't get "better"
* cells with a stressor with an increasing trend, but still below the threshold, can't be increasing all that rapidly and do not seem to be a priority.

## 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}

### Get inclusion list; drop species with no sensitivities to stressors
spp_sens <- get_incl_spp() %>%
  filter(!is.na(stressor)) %>%
  mutate(taxon = str_replace_all(desc, ' ', ''))

### directories for inputs
spp_range_dir <- file.path(dir_bd_anx, 'spp_rasts_mol_2020')

### directories for outputs
taxon_trend_map_dir <- file.path(dir_bd_anx, 'taxon_trend_rasts')
### unlink(file.path(taxon_trend_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 trend maps

These maps count up the number of species with increasing impacts by the aggregated stressor group (land-based, ocean, climate, fishing, and all); similarly for species with decreasing impacts.

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 trend maps to increasing or decreasing or neutral
    * clip trend map to footprint for 2013
    * add all maps for each category group: by cell/year per stressor
* for all subgroups in a taxa, add them up.
    
```{r set up functions and dataframes}

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_rangemap <- function(spp_id) {
  range_map_filestem <- file.path(spp_range_dir, 'iucn_sid_%s.csv')
  rangemap <- read_csv(sprintf(range_map_filestem, spp_id),
                       col_types = 'ii') %>%
    mutate(iucn_sid = spp_id)
  return(rangemap)
}

get_str_footprint <- function() {
  ### This retrieves the stressor footprints for the 2013 year, and
  ### converts the rasters to a footprint dataframe of cell IDs
  message('Retrieving stressor footprints...')
  str_vec <- get_str_cats() %>%
    .$stressor
  ### set up file stem for stressor footprint rasters
  str_footprint_dir <- file.path(dir_bd_anx, 'layers/stressors_10km_95vc')
  str_ft_stem <- file.path(str_footprint_dir, '%s_2013_mol10km_95vc.tif')
    ### stressor - use 2013 footprint for clipping
  
  str_ft <- raster::stack(sprintf(str_ft_stem, str_vec))
  str_ft_df <- data.frame(values(str_ft)) %>%
    setNames(str_vec) %>%
    mutate(cell_id = 1:ncell(str_ft)) %>%
    gather(stressor, value, -cell_id) %>%
    filter(!is.na(value)) %>%
    select(-value)
  
  return(str_ft_df)
}

### function to read in all the pace files
get_pace_df <- function() {
  pace_fs <- list.files(file.path(dir_bd_anx, 'layers/stressor_pace'), 
                        full.names = TRUE)
  pace_df <- parallel::mclapply(pace_fs, FUN = function(f) {
        df <- read_csv(f, col_types = c(cell_id = 'i', pace = 'i', stressor = 'c'))
      }, mc.cores = 14) %>%
    bind_rows()
}
```

```{r create taxa maps}

# x <- list.files(taxon_trend_map_dir, full.names = TRUE)
# # unlink(x)

### Trend files are rasters of value (trend slope) and p_val.  This captures
### pace in three levels: trend > 0, trend > .001, and trend > .01 (and 
### less than neg values for decrease).
pace_df <- get_pace_df()

### get the stressor footprints for 2013 as a dataframe
str_ft_df <- get_str_footprint()

str_cats <- c('land-based', 'ocean', 'fishing', 'climate', 'all')
taxon_gps <- spp_sens$taxon %>% unique() %>% sort()

subgp_size <- 24 ### number of spp to process at a time
reload <- FALSE

for(str_cat in str_cats) {
  # str_cat <- str_cats[3]
  
  for(taxon_gp in taxon_gps) {
    ### taxon_gp <- taxon_gps[1]
    
    tx_incr_f <- file.path(taxon_trend_map_dir,
      sprintf('taxon_trend_%s_%s_incr%s.tif', taxon_gp, str_cat, 1:3))
    tx_decr_f <- file.path(taxon_trend_map_dir, 
      sprintf('taxon_trend_%s_%s_decr%s.tif', taxon_gp, str_cat, 1:3))

    if(any(!file.exists(tx_incr_f, tx_decr_f)) | 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, ' intensification 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_trend_df <- data.frame()
      } else {
        ### loop over all spp in this taxon affected by this stressor
        ### initialize a list to store 
        taxon_trend_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)
          
          spp_gp_sens <- spp_sens %>%
            filter(iucn_sid %in% spp_id_vec) %>%
            filter(str_cat == 'all' | category == str_cat) %>%
            select(iucn_sid, stressor) %>%
            distinct()
            
          spp_ranges <- parallel::mclapply(
              spp_id_vec,
              FUN = get_spp_rangemap,
              mc.cores = subgp_size) %>%
            bind_rows()
          
          trend_gp <- spp_ranges %>%
            left_join(spp_gp_sens, by = 'iucn_sid') %>%
            # left_join(pace_df, by = 'cell_id')
            dt_join(pace_df, by = c('cell_id', 'stressor'), type = 'inner') %>%
              ### inner join drops non-intensifying/deintensifying cells
            dt_join(str_ft_df, by = c('cell_id', 'stressor'), type = 'inner')
              ### inner join drops all cells where stressor is not "present"
              ### according to 95% contour volume
          
          trend_gp_summary <- trend_gp %>%
            group_by(cell_id, iucn_sid) %>%
            summarize(incr1 = any(pace %in% 1:3) & !any(pace < 0),
                      decr1 = any(pace %in% -1:-3) & !any(pace > 0),
                      incr2 = any(pace %in% 2:3) & !any(pace < 0),
                      decr2 = any(pace %in% -2:-3) & !any(pace > 0),
                      incr3 = any(pace %in% 3) & !any(pace < 0),
                      decr3 = any(pace %in% -3) & !any(pace > 0)) %>%
            group_by(cell_id) %>%
            summarize(n_incr1_spp = sum(incr1, na.rm = TRUE),
                      n_decr1_spp = sum(decr1, na.rm = TRUE),
                      n_incr2_spp = sum(incr2, na.rm = TRUE),
                      n_decr2_spp = sum(decr2, na.rm = TRUE),
                      n_incr3_spp = sum(incr3, na.rm = TRUE),
                      n_decr3_spp = sum(decr3, na.rm = TRUE)) %>%
            ungroup()
          
          taxon_trend_list[[i]] <- trend_gp_summary
        } ### end of for loop
        
      ### collate all the list element dataframes
      taxon_trend_df <- bind_rows(taxon_trend_list) %>%
        group_by(cell_id) %>%
        summarize(n_incr1_spp = sum(n_incr1_spp, na.rm = TRUE),
                  n_decr1_spp = sum(n_decr1_spp, na.rm = TRUE),
                  n_incr2_spp = sum(n_incr2_spp, na.rm = TRUE),
                  n_decr2_spp = sum(n_decr2_spp, na.rm = TRUE),
                  n_incr3_spp = sum(n_incr3_spp, na.rm = TRUE),
                  n_decr3_spp = sum(n_decr3_spp, na.rm = TRUE)) %>%
        ungroup() %>%
        mutate(taxon = taxon_gp)
      
      } ### end of if statement for zero-length species vector
      
      if(nrow(taxon_trend_df) == 0) {
        ### either no species in the list, or no trend for any of the spp
        taxon_trend_df <- data.frame(cell_id = 1,
                                     n_incr1_spp = 0,
                                     n_decr1_spp = 0,
                                     n_incr2_spp = 0,
                                     n_decr2_spp = 0,
                                     n_incr3_spp = 0,
                                     n_decr3_spp = 0,
                                     taxon   = taxon_gp)
      }
      
      ### convert to rasters and write out as +/- by count and priority

      nspp_incr1_rast <- map_to_rast(taxon_trend_df, cell_val = 'n_incr1_spp')
      nspp_decr1_rast <- map_to_rast(taxon_trend_df, cell_val = 'n_decr1_spp')
      nspp_incr2_rast <- map_to_rast(taxon_trend_df, cell_val = 'n_incr2_spp')
      nspp_decr2_rast <- map_to_rast(taxon_trend_df, cell_val = 'n_decr2_spp')
      nspp_incr3_rast <- map_to_rast(taxon_trend_df, cell_val = 'n_incr3_spp')
      nspp_decr3_rast <- map_to_rast(taxon_trend_df, cell_val = 'n_decr3_spp')
 
      writeRaster(stack(nspp_incr1_rast, nspp_incr2_rast, nspp_incr3_rast),
                  tx_incr_f, bylayer = TRUE, overwrite = TRUE)
      writeRaster(stack(nspp_decr1_rast, nspp_decr2_rast, nspp_decr3_rast),
                  tx_decr_f, bylayer = TRUE, overwrite = TRUE)

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

```

## Plot select taxa trend maps 

Plot taxa maps species experiencing intensifying/deintensifying impacts.  Here we plot only the trends for fishing stressors, climate stressors, and all cumulative stressors, for a subset of taxa.  We plot the count only (not priority-weighted).

Note, since we have three levels of intensification, here we will print all levels for a few select taxa and stressor categories...

```{r trend maps by taxa}

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

strs <- c('fishing', 'climate', 'all')

taxa <- c('mammals', 'marinereptiles', 'bonyfishes', 'sharksandrays', 'seabirds')

for(str_cat in strs) {
  # str_cat <- strs[1]
  for(t in taxa) {
    # t <- taxa[3]
    fstem <- file.path(taxon_trend_map_dir, 'taxon_trend_%s_%s_%s%s.tif')
    
    incr_map <- stack(sprintf(fstem, t, str_cat, 'incr', 1:3))
    plot(incr_map, 
         col = c('grey70', hcl.colors(n = 20, rev = TRUE)),
         main = sprintf('incr %s, %s, lvl %s', str_cat, t, 1:3))
    
    decr_map <- stack(sprintf(fstem, t, str_cat, 'decr', 1:3))
    plot(decr_map, 
         col = c('grey70', hcl.colors(n = 20, rev = TRUE)),
         main = sprintf('decr %s, %s, lvl %s', str_cat, t, 1:3))
    
    overlap <- stack(incr_map)
    values(overlap) <- NA
    values(overlap)[values(incr_map) > 0 & values(decr_map) > 0] <- 1
    cat('sum of overlap of', str_cat, 'for', t, 'across all levels: ', sum(values(overlap), na.rm = TRUE), 'cells\n')
  }
}
```