prioritizr_demo.Rmd

---
title: "EcoDataScience Workshop: Conservation Planning with Prioritizr"
author: "Anna Abelman, Vanessa Rathbone, & Rachel Rhodes"
date: "5/7/2021"
output: 
  html_document:
    code_folding: show
    toc: true
    toc_depth: 3
    toc_float: yes
    number_sections: false
    theme: cerulean
    highlight: haddock
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

### Load packages
library(raster)
library(tidyverse)
library(sf)
library(here)
library(fasterize)
library(stringr)
library(janitor)
library(prioritizr)
library(rgdal)
library(slam)
library(gurobi)

```

# Session summary
For this EcoDataScience session, we will cover some basics of working with  `prioritizr` R package to build and solve conservation planning problems using marine spatial data presented in gridded format (rasters) and vector format (line, point, polygons). We will go through the basic functionality and necessary inputs to create and solve a conservation planning problem, as well as some guidance on how you can custom-tailor this to the specific needs of different conservation planning exercises. For brevity, we will not go through all the pre-processing steps of cleaning up the raw spatial data we used for this analysis, but will include a separate set-up RMD that outlines what we did which you are welcome to check out on your own. 

The analysis presented in this tutorial comes from a USCB Bren School of Environmental Science & Management master's project by Anna Abelman, Courtney Krone, Vanessa Rathbone, Rachel Rhodes & Erin Ristig in collaboration with Wildlife Conservation Society. The project aimed to identify priority areas to protect sharks and rays in Mozambique. For more information on the project you can visit: https://bren.ucsb.edu/projects/framework-designing-marine-protected-areas-sharks-and-rays-mozambique

# Background
## What is Prioritizr?
The *prioritizr* R package is a systematic conservation planning tool that is used to design an optimal system of protected areas that meet conservation objectives using ecological, social, and economic criteria. Specifically, this tool is used to conduct spatial prioritization that finds the efficient spatial allocation of priority areas by formulating a mathematical optimization problem and then solving it to find an efficient solution. This type of tool is useful for conservation planning because it can help inform decisions on how and where to focus conservation efforts given limited resources and can make for a more transparent and reproducible decision making process. 

For those of you familiar with other conservation planning software like Marxan, *prioritizr* is similar, however, instead of using simulated annealing, *prioritizr* uses integer linear programming (ILP) and an algorithm solver to find the exact optimal solution and can help find cheaper solutions in a shorter period of time. It requires similar inputs which are described below. 

## Necessary Inputs
The following parameters are required to conduct and solve a conservation planning exercise with *prioritizr*: 

1. Study Area: define area of interest

2. Planning Units: divide the study area into set of discrete areas called planning units

3. Cost: specifies the cost of including/managing each planning unit in the reserve system.  This value can be a measure of the actual fiscal cost required to purchase a piece of land or it can be a surrogate for cost that represents management costs or opportunity costs of foregone activities. 

4. Conservation Features: biodiversity elements that are of conservation interest - often include species, populations, or habitats 

5. Problem objective: defines the overall goal of the conservation plan & whether a specific property of the solution should be maximized or minimized (there are a number of different types of objectives that can be used which are described more later)

## Example data

**Species Data:** 
IUCN 2021. The IUCN Red List of Threatened Species. Version 2021-1. <https://www.iucnredlist.org>

**Habitat:** 
Seamount & Knolls: Yesson C, Clark MR, Taylor M, Rogers AD (2011). The global 
distribution of seamounts based on 30-second bathymetry data. Deep Sea Research Part I: Oceanographic Research Papers 58: 442-453. DOI: http://dx.doi.org/10.1016/j.dsr.2011.02.004. Data DOI: http://doi.pangaea.de/10.1594/PANGAEA.757564

**Fishing Pressure:** 
Global Fishing Watch
D.A. Kroodsma, J. Mayorga, T. Hochberg, N.A. Miller, K. Boerder, F. Ferretti, A. Wilson, B. Bergman, T.D. White, B.A. Block, P. Woods, B. Sullivan, C. Costello, and B. Worm. "Tracking the global footprint of fisheries." Science 361.6378 (2018).

# Basic Analysis

## Step 1. Project Area & Planning  Units
The first step is to define the project area - for this analysis we will use the Mozambique Exclusive Economic Zone, and divide the area up into planning units. 

The size of the planning unit can impact the outcomes of the reserve design, so it is
important to carefully consider the spatial scale of the planning unit - for this exercise we will rasterize the project area to produce a grid of cells at 10km^2 resolution. 

```{r}
# Load raster of Mozambique Exclusive Economic Zone (EEZ) 
mozambique_eez <- raster(here("cleaned_data", "mz_eez_rast.tif"))
```


## Step 2. Conservation Features 
For this exercise we will use spatial data of 5 shark and ray species distribution from IUCN and two critical habitats that are important for sharks and rays. For brevity, we will load in the species and critical habitat spatial data that has already been pre-processed (rasterized & projected in the same coordinate system) - if you would like to see how we did this you can check out the `set-up.RMD`. 

```{r}
# Pull in species & habitat rasters from 'cleaned_data' folder
## Ray species
devil_ray_rast <- raster(here("cleaned_data", "devil_ray.tif"))
sawfish_rast <- raster(here("cleaned_data", "sawfish.tif"))
reef_manta_rast <- raster(here("cleaned_data", "reef_manta.tif"))

## Shark species
scalloped_hammerhead_rast <- raster(here("cleaned_data", "scalloped_hammerhead.tif"))
dusky_shark_rast <- raster(here("cleaned_data", "dusky_shark.tif"))

## Habitat
knolls_habitat <- raster(here("cleaned_data", "knolls.tif"))
seamounts_habitat <- raster(here("cleaned_data", "seamounts.tif"))

# Create a raster stack of all the conservation features
## note to create a raster stack, all rasters must have the same extent and resolution
feature_stack <- stack(devil_ray_rast, sawfish_rast, reef_manta_rast, scalloped_hammerhead_rast, dusky_shark_rast, knolls_habitat, seamounts_habitat)

# Plot to see what it looks like
plot(feature_stack)
```


## Step 3. Cost
We will use fishing pressure as a surrogate for cost in our model which represents the cost of lost fishing opportunity. This will prioritize areas with the least fishing activity that overlap with areas that are critical for sharks and rays - this will show us the reserve designs that have the least impact on fishers (alternatively you could create an inverse of this cost layer to prioritize areas with the most fishing activity that overlap with critical areas for sharks and rays which would have the most benefit for sharks and rays, but greater impacts on fishers).

```{r}
## Pull in the fishing pressure data from the 'cleaned_data' folder
fishing_cost <- raster(here("cleaned_data", "fishing_stack.tif"))

## Plot to see what it looks like
plot(fishing_cost)

## Create data frame to see what the data looks like - cost ranges from 0.001 - 25
fishing_cost_df <- as.data.frame(fishing_cost)

```


## Step 4. Create conservation problem
Now that we have these three basic inputs, it is time to start constructing our conservation problem. 

There are a few key considerations that go into this including: 
(1) what type of problem objective we want to use 
(2) what are the conservation targets (minimum proportion  of each conservation feature that need to be included int he reserve system)?

The problem objective defines the overall goal of the conservation plan. There are different objectives for different conservation problems. For our analysis, we'll use the minimum set objective, which minimizes the cost of our reserve system, while meeting our conservation targets that we specify: `add_min_set_objective()`

The conservation targets are the minimum proportion of our conservation features’ distribution that needs to be protected and included in the prioritization (or reserve design). These targets can be thought of as a quantitative interpretation of our conservation goal, and thus should be set at the amount required to ensure the long-term persistence of our focal species, habitat, or population. For this example, we will use a uniform conservation target of 10%  across all features, so we are saying the reserve system must include at least 10% of each species distribution and at least 10% of each critical habitat: `add_relative_targets(0.1)`

Now that we have decided those two things, we can construct the conservation problem:

```{r}
# Define the problem using the cost, conservation features, objective, and targets
prob_1 <- problem(fishing_cost, feature_stack) %>% 
  add_min_set_objective() %>% 
  add_relative_targets(0.1) %>% 
  add_gurobi_solver(gap = 0.01, time_limit = 2100)

# Print problem
print(prob_1)
```

## Step 5. Solve the problem 
To solve the problem, we will run *prioritizr* using a solver or optimization software. The recommended and most common one used is Gurobi, which has special licenses available to academics at no cost. If you have not already downloaded this, please see the README.md for instructions. 
```{r}
## Solve the problem
sprob_1 <- solve(prob_1)

## plot the solution to see what it looks like
plot(sprob_1,  main = c("10% Targets"))
```

#----------------------------#
# Variations to the Analysis 
#----------------------------#

Now that you've defined your baseline problem, you can adjust your model with different variations. Here, we'll walk through a few of those variables you can add. 

## Changing Targets 
```{r}
# Let's change our prob_1 target to include 30% of the conservation features, instead of 10%

## Baseline -30% target
prob_30 <- problem(fishing_cost, features = feature_stack) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.3)

## Solve problem
sprob_30 <- solve(prob_30)

## Plot the solution to see what it looks like
plot(sprob_30,  main = c("Area- 30% Targets"))

```

## Locking-in/Locking-out
```{r}
# Here, we lock-in existing Marine Protected Areas

## Read in existing mpas raster
exist_mpas <- raster(here("cleaned_data", "mpa_rast.tif"))

## Plot to make sure it looks okay
plot(exist_mpas)

## Run model with locked-in constraints: 
prob_30_mpa <- problem(fishing_cost, features = feature_stack) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.3) %>%
  add_locked_in_constraints(exist_mpas)

## Solve problem
sprob_30_mpa <- solve(prob_30_mpa)

## Plot the solution to see what it looks like
plot(sprob_30_mpa,  main = c("Area- 30% Targets with MPAs locked-in"))
plot(sprob_30)

```

### Setting some global model parameters

Here, we'll be setting:
-*the gap* (the optimality gap tolerance you are willing to allow the solver to declare as an optimal solution, even though the model may have other, better solutions - we want the gap to be as small as possible)
-*time limit* (how long your model will run)
-*efficient boundary penalty* (optimal penalty to reduce fragmentation within your reserves)

Boundary penalty: Sometimes you might want to avoid fragmentation in your design for management feasibility or designing for mobile and migratory species where larger reserves will be more effective for conservation than smaller fragmented reserves.To avoid highly fragmented solution, you use a boundary penalty. But know that using a boundary penalty is trade off -  by setting a high boundary penalty the model will prioritize solutions that are less fragmented even if they cost more. 

```{r}

## If you decide to change the gap, time-limit or boundary penalty value, change it here and rerun the subsequent model runs
# Here, setting both a relative gap and time limit the solver will terminate after 24 hours OR when the relative gap is less than 1%

gap <- 0.01 
time_limit <- 120
bp <- 0.000001

# Run the model with all of the new model parameters:
prob_30 <- problem(fishing_cost, features = feature_stack) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.3) %>%
  add_locked_in_constraints(exist_mpas) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit) %>% 
  add_boundary_penalties(bp)

## Solve problem
sprob_30 <- solve(prob_30)

## Plot the solution to see what it looks like
plot(sprob_30,  main = c("Area- 30% Targets with MPAs locked-in with BP"))
plot(sprob_30_mpa)
```


## Other Considerations (NEW TO prioritizr)

#Summary stats use 'eval_target_coverage_summary(prob_30, sprob_30)':
To evaluate the performance of the solution, we can calculate summary statistics.

#For Contiguious Reserve Designs use `add_contiguity_constraints`
To maintain connectivity and dispersal capacity between reserves, we can add contiguity constraints.

#For areas that MUST be considered in the design use `irreplaceability`:
Some areas in your design might have be included in your design. These areas can now be calculated as "irreplaceable". The `irreplaceability()` function will calculate importance (irreplaceability) scores using a version of the replacement cost method. Under this method, planning units with higher scores are more important for meeting the objective of our conservation planning problem than those with lower scores. Note that we override the solver behavior in the code below to prevent lots of unnecessary text from being output.

## Different Conservation Problems

Can find more details here: https://cran.r-project.org/web/packages/prioritizr/vignettes/prioritizr.html#initialize-a-problem

*Minimum set objective:* Minimize the cost of the solution whilst ensuring that all targets are met.
*Maximum features objective:* Fulfill as many targets as possible while ensuring that the cost of the solution does not exceed a budget.
*Minimum shortfall objective:* Minimize the overall (weighted sum) shortfall for as many targets as possible while ensuring that the cost of the solution does not exceed a budget.
*Minimum largest shortfall objective:* Minimize the largest (maximum) shortfall while ensuring that the cost of the solution does not exceed a budget.
*Maximum phylogenetic diversity objective:* Maximize the phylogenetic diversity of the features represented in the solution subject to a budget 
*Maximum phylogenetic endemism objective:* Maximize the phylogenetic endemism of the features represented in the solution subject to a budget
*Maximum utility objective:* Secure as much of the features as possible without exceeding a budget.