5a_map_impacts.Rmd

---
title: 'Map impacts for all impacted species'
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'))

library(animation)

reload <- FALSE
```


# Summary

Create a species richness map for threatened marine species.

Combine taxa-level impact  maps to create maps of number of threatened species impacted per cell, per year.  This is also done by range priority (see the `get_spp_priority()` function for details)

Plot maps as animated rasters, both by absolute numbers of species impacted in a location and by proportion of species (impact counts / species richness).

This is all done at the level of cumulative impact category, not individual stressors.

# Methods

## Create species richness map for this assessment

Set up the dataframe of species to include.  By taxon, load all species range maps, smash down to taxon-level species richness; then combine these for a total species richness.

```{r Identify species in this assessment}

spp_incl <- get_incl_spp() %>%
  mutate(taxon = tolower(rank))

spp_sens <- spp_incl %>%
  filter(!is.na(stressor))

# spp_incl$iucn_sid %>% n_distinct()
# spp_incl %>% filter(!is.na(stressor)) %>% .$iucn_sid %>% n_distinct()
### 1036 if we limit to just sensitive spp, 1271 if we include all
### spp that are mapped/threatened/comp assessed
```


``` {r create species richness map}

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

spp_rich_mapfile <- here('_output/rasters/n_spp_map.tif')

if(!file.exists(spp_rich_mapfile) | reload) {
  
  spp_range_dir <- file.path(dir_bd_anx, 'spp_rasts_mol_2020')
  
  taxa_gps <- spp_incl$taxon %>% unique() %>% sort()
  
  taxa_map_list <- vector('list', length = length(taxa_gps))
  for(i in seq_along(taxa_gps)) {
    # i <- 2
    taxa_gp <- taxa_gps[i]
    spp_ids_taxa <- spp_incl %>%
      filter(taxon == taxa_gp) %>%
      .$iucn_sid %>%
      unique()
    
    map_list <- parallel::mclapply(spp_ids_taxa, mc.cores = 12,
      FUN = function(id) {
        # id <- spp_ids_taxa[1]
        spp_f <- file.path(spp_range_dir, sprintf('iucn_sid_%s.csv', id))
        spp_df <- read_csv(spp_f, col_types = 'ii') %>%
          mutate(iucn_sid = id) %>%
          mutate(present = (presence != 5) & !is.na(presence))
      })
    map_df <- bind_rows(map_list)
    map_nspp_df <- map_df %>%
      group_by(cell_id) %>%
      summarize(n_spp = sum(present))
    taxa_map_list[[i]] <- map_to_rast(map_nspp_df, cell_val = 'n_spp')
  }
  
  taxa_map_stack <- stack(taxa_map_list)
  spp_rich_rast <- raster::calc(taxa_map_stack, 
                                fun = sum, na.rm = TRUE) %>%
    mask(ocean_a_rast)
  
  writeRaster(spp_rich_rast, spp_rich_mapfile, overwrite = TRUE)
}

spp_rich_rast <- raster(spp_rich_mapfile)
plot(log(spp_rich_rast), 
     axes = F, col = hcl.colors(n = 20),
     main = 'log(Species richness)')

```


## Setup for mapping across all taxa

For each impact category (including all):

* gather all taxa-level map files into a raster stack
* collapse the stack using `calc(fun = sum, na.rm = TRUE)`
* animate over the years

```{r identify taxa map files}
spp_range_dir <- file.path(dir_bd_anx, 'spp_rasts_mol_2020')
impact_map_dir <- here('_output/rasters/impact_maps')
taxon_impact_map_dir <- file.path(dir_bd_anx, 'taxon_impact_rasts')

all_maps <- list.files(taxon_impact_map_dir,
                       pattern = '.tif',
                       full.names = TRUE)

taxa_map_df <- data.frame(map = all_maps) %>%
  mutate(year   = str_extract(basename(map), '[0-9]{4}') %>% as.integer(),
         impact = str_replace_all(basename(map), '.+[0-9]{4}_|.tif', ''))
```


```{r animate function}
make_ramp <- function(n, palette) {
  if(n > 101) {
    n_even <- round(n + 50, -2)
    breaks <- seq(.5, n_even, length.out = 100)
    labels <- seq(0, n_even, by = 100)
  } else {
    n_even <- n
    breaks <- seq(min(0.5, n / 100), n, length.out = 100)
    labels <- seq(0, n, by = (n / 10))
  }
  colors <- hcl.colors(length(breaks), palette = palette)
  
  ### For breaks and colors, add a different color for 
  ### a value of exactly zero
  return(list('colors' = c('grey40', colors),
              'breaks' = c(0, breaks),
              'labels' = labels))
}

make_gifs <- function(map_stack, filename, layer_names = NULL) {
  if(is.null(layer_names)) layer_names = names(map_stack)
  
  ramp <- make_ramp(max(maxValue(map_stack)), palette = 'viridis')
  
  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$labels),
             main = layer_names[i])
      }}, 
      interval = 0.5, movie.name = filename, 
      ani.width = 700, ani.height = 420)
  })
  
  return(invisible(NULL))
}

nspp_rast <- spp_rich_rast
values(nspp_rast)[values(nspp_rast) == 0] <- NA
```

## Create impact map for all included spp, by stressor category

```{r impact all species by stressor category}

str_cats <- c('fishing', 'climate', 'land-based', 'ocean', 'all')
yrs <- 2003:2013

outfile_stem <- file.path(impact_map_dir, 'impact_%s_%s.tif')

for(str_cat in str_cats) {
  # str_cat <- str_cats[1]
  impact_map_files <- taxa_map_df %>%
    filter(impact == str_cat)
  
  for(y in yrs) {
    # y <- 2003
    outfile <- sprintf(outfile_stem, str_cat, y)
    # message('Processing ', basename(outfile))
    
    if(!file.exists(outfile) | reload) {
      imp_yr_files <- impact_map_files %>%
        filter(year == y) %>%
        .$map
      impact_yr_stack <- stack(imp_yr_files)
      impact_yr_rast <- raster::calc(impact_yr_stack, 
                                     fun = sum, na.rm = TRUE) %>%
        mask(ocean_a_rast)
      
      writeRaster(impact_yr_rast, outfile, overwrite = TRUE)
      
    } #else {
      # message('File exists: ', outfile, '... skipping!')
    # }
  }
}

```

```{r impact all species by stressor category 2plus stressors 2013 only}

outfile <- file.path(impact_map_dir, 'impact_all_2013_2plus.tif')


if(!file.exists(outfile) | reload) {
      
  impact_map_files <- taxa_map_df %>%
    filter(impact == 'all_2plus')
  
  impact_2plus_stack <- stack(impact_map_files$map)
  impact_2plus_rast <- raster::calc(impact_2plus_stack, 
                                 fun = sum, na.rm = TRUE) %>%
    mask(ocean_a_rast)
  
  writeRaster(impact_2plus_rast, outfile, overwrite = TRUE)
    
}

```

<!-- ### maps by count of threatened species impacted -->

``` {r Animate impact count maps, results = 'asis', eval = FALSE, echo = FALSE}
for(str_cat in str_cats) {
  # str_cat <- str_cats[1]
  
  gif_file <- here(sprintf('figs/impact_all_spp_%s.gif', str_cat))
  
  if(!file.exists(gif_file) | reload) {

    rast_files <- sprintf(outfile_stem, str_cat, yrs)
    map_stack <- stack(rast_files)
    
    make_gifs(map_stack, 
              filename = gif_file,
              layer_names = sprintf('Impacted spp (#): %s (%s)', str_cat, yrs))
  }
  
  # knitr::include_graphics(gif_file)
  cat(sprintf('![](%s)', gif_file))
  
}
```

### maps by percent of threatened species impacted

``` {r Animate impact pct maps, results = 'asis'}
for(str_cat in str_cats) {
  # str_cat <- str_cats[1]
  
  gif_file <- here(sprintf('figs/impact_pct_all_spp_%s.gif', str_cat))
  
  if(!file.exists(gif_file) | reload) {
    rast_files <- sprintf(outfile_stem, str_cat, yrs)
    map_stack <- stack(rast_files)
    x <- map_stack / nspp_rast
    
    if(any(values(x) > 1, na.rm = TRUE)) {
      stop('pct raster values greater than 1, wtf')
    }
    
    values(x)[values(x) > 1] <- 1
    x <- x * 100
    
    make_gifs(x, 
              filename = gif_file,
              layer_names = sprintf('Impacted spp (percent): %s (%s)', str_cat, yrs))
  }
  
  # knitr::include_graphics(gif_file)
  cat(sprintf('![](%s)', gif_file))

}
```

### Plot 2013 layers

``` {r plot just 2013 layers}
str_cats <- c('fishing', 'climate', 'land-based', 'ocean', 'all')
outfile_stem <- file.path(impact_map_dir, 'impact_%s_%s.tif')

for(str_cat in str_cats) {
  # str_cat <- str_cats[1]
  r <- raster(sprintf(outfile_stem, str_cat, 2013))
  # r_pri <- raster(sprintf(outfile2_stem, str_cat, 2013))
  
  r_norm <- r / nspp_rast
  # r_pri_norm <- r_pri / priority_sum_rast
    
  # plot(r, col = hcl.colors(n = 20),
  #      axes = FALSE,
  #      main = paste('count spp impacted:', str_cat, 'stressors'))
  plot(r_norm, col = hcl.colors(n = 20),
       axes = FALSE,
       main = paste('pct spp impacted:', str_cat, 'stressors'))
  # plot(r_pri_norm, col = hcl.colors(n = 20),
  #      axes = FALSE,
  #      main = paste('priority pct spp impacted:', str_cat, 'stressors'))
  
}
```

```{r plot 2plus stressor map}
r <- raster(sprintf(outfile_stem, 'all', '2013_2plus'))

r_norm <- r / nspp_rast

plot(r_norm, col = hcl.colors(n = 20),
     axes = FALSE,
     main = paste('pct spp impacted: all stressors (2+)'))
```

<!-- How correlated are impact counts with species richness?  Examine 2013 year. -->

```{r, eval = FALSE}

for(str_cat in str_cats) {
  # str_cat <- str_cats[1]
  r_stem <- here('_output/rasters/impact_maps/impact_%s_2013.tif')
  imp_r <- raster(sprintf(r_stem, str_cat))
  df <- data.frame(impact = values(imp_r),
                   nspp   = values(nspp_rast)) %>%
    filter(!is.na(nspp))
  
  df1 <- sample_frac(df, .02)
  thingplot <- ggplot(df1, aes(x = nspp, y = impact)) +
    geom_point(alpha = .2) +
    geom_abline(slope = 1, intercept = 0, color = 'red') +
    geom_smooth(method = 'lm', color = 'green', alpha = .5) +
    labs(title = sprintf('%s impacts vs spp richness (2 pct sample)', str_cat))
  print(thingplot)
}

```