scenario_4_combined.Rmd

---
title: "Scenario 4 - Combined Fishing"
output: html_document
---

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

library(raster)
library(tidyverse)
library(sf)
library(here)
library(fasterize)
library(stringr)
library(janitor)
library(prioritizr)
library(rgdal)
library(gurobi)
```

## Scenario 4 - Combined Fishing Pressure 
The objective of this scenario was to protect important areas for sharks and rays based on their distribution from IUCN and location of suitable habitats while minimizing the impact on both artisanal and industrial fisherman.  For this run, we used area as the cost, set targets of 20% for all conservation features, locked in planning units that included existing MPAs, and set a boundary penalty to zero. 

Step 1 reads in all the necessary data
Step 2 run to see results of varying targets from 10-50%
Step 3 to see varying targets using a boundary penalty
Step 4 to run sensitivity analysis of locking in vs locking out MPAs
Step 5 to run sensitivity analysis of key biodiversity areas
Step 6 to run sensitivity analysis of aggregation sites

----

## **STEP 1: Assign planning units, cost layer, and conservation feature**

#### Planning Unit
```{r}
### reading in the baseline MZ EEZ raster (with cell ids) 
pu_4 <- raster(here("set_up", "final_rasters", "mz_eez_templates", "mz_rast_id.tif"))

### plot to make sure it looks okay
plot(pu_4)
```

#### Cost Layer
```{r}
### pull in the combined fishing pressure tif file from costs folder in repo
cost_combined <- raster(here("set_up", "final_rasters", "costs", "combined_stack_sum.tif"))

### plot to make sure it looks okay
plot(cost_combined)
```

#### Conservation Features
This includes all species distributions and habitat rasters
```{r}
# List of conservation features that will have been loaded into the global environment in other
### pull in species stack that we created in set-up step 3
features_stack_4 <- brick(here("set_up", "final_rasters", "conservation_features", "all_features_stack.tif"))

### plot to make sure it looks okay
plot(features_stack_4)
```

### Additional Files (Locked-in & Locked-out)

Marine Protected Areas:
```{r}
### read in existing mpas raster
exist_mpas <- raster(here("set_up", "final_rasters", "planning_unit", "exist_mpa.tif"))

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

Key Biodiversity Areas:
```{r}
### read in existing kbas raster
exist_kbas <- raster(here("set_up", "final_rasters", "planning_unit", "kba.tif"))

### plot to make sure it looks okay
plot(exist_kbas)
```

Aggregation Sites:
```{r}
### read in existing aggregation sites raster
exist_aggs <- raster(here("set_up", "final_rasters", "planning_unit", "aggregation_sites.tif"))

### plot to make sure it looks okay
plot(exist_aggs)
```


### Set the some of the global model parameters like gap, time limit, and efficient boundary penalty to use in soem runs:
```{r}
# if you decide to change the gap, time-limit or boundary penalty value, change it here and rerun the subsequent model runs
gap <- 0.01
time_limit <- 21600
bp <- 0.000001
```


## **STEP 2: Run Models with varying targets**

#### RUN 1: Combined - 10% target, Boundary Penalty = 0 
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_10 <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.1) %>%
  add_locked_in_constraints(exist_mpas) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit)

## Solve problem
sprob_combined_10 <- solve(prob_combined_10)

## Plot the solution to see what it looks like
plot(sprob_combined_10,  main = c("combined- 10% Targets"))

## Save as tif
writeRaster(sprob_combined_10, here("final_results", "scenario_4_combined", "sprob_combined_10_6hr_blm0.tif"), overwrite = TRUE)
```


#### RUN 2: Combined - 20% target, Boundary Penalty = 0
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_20 <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.2) %>%
  add_locked_in_constraints(exist_mpas) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit) 

## Solve problem
sprob_combined_20 <- solve(prob_combined_20)

## Plot the solution to see what it looks like
plot(sprob_combined_20,  main = c("combined- 20% Targets"))

## Save as tif
writeRaster(sprob_combined_20, here("final_results", "scenario_4_combined", "sprob_combined_20_6hr_blm0.tif"), overwrite = TRUE)

```

#### RUN 3: Combined - 30% target, Boundary Penalty = 0,  
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_30 <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.3) %>%
  add_locked_in_constraints(exist_mpas) %>%
  add_gurobi_solver(gap = 0.01, time_limit = 120) 

## Solve problem
sprob_combined_30 <- solve(prob_combined_30)

## Plot the solution to preview
plot(sprob_combined_30,  main = c("combined- 30% Targets"))

## Save as tif
writeRaster(sprob_combined_30, here("final_results", "scenario_4_combined", "sprob_combined_30_6hr_blm0.tif"), overwrite = TRUE)

```

#### RUN 4: Combined - 50% target, Boundary Penalty = 0,  
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_50 <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.5) %>%
  add_locked_in_constraints(exist_mpas) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit) 

## Solve problem
sprob_combined_50 <- solve(prob_combined_50)

## Plot the solution to preview
plot(sprob_combined_50,  main = c("combined- 50% Targets"))

## Save as tif
writeRaster(sprob_combined_50, here("final_results", "scenario_4_combined", "sprob_combined_50_6hr_blm0.tif"), overwrite = TRUE)

```


## **STEP 3: Run Models with efficient boundary penalty (0.000001)**

#### Run 5: Combined - 10% target, Boundar Penaly = 0.000001
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_10_blm <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.1) %>%
  add_locked_in_constraints(exist_mpas) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit) %>% 
  add_boundary_penalties(bp)

## Solve problem
sprob_combined_10_blm <- solve(prob_combined_10_blm)

## Plot the solution to see what it looks like
plot(sprob_combined_10_blm,  main = c("combined - 10% Targets"))

## Save as tif
writeRaster(sprob_combined_10_blm, here("final_results", "scenario_4_combined", "sprob_combined_10_6hr_blm0.000001.tif"), overwrite = TRUE)
```


#### RUN 6: Combined - 20% target, Boundary Penalty = 0.000001 
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_20_blm <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.2) %>%
  add_locked_in_constraints(exist_mpas) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit) %>% 
  add_boundary_penalties(bp)

## Solve problem
sprob_combined_20_blm <- solve(prob_combined_20_blm)

## Plot the solution to see what it looks like
plot(sprob_combined_20_blm,  main = c("combined - 20% Targets"))

## Save as tif
writeRaster(sprob_combined_20_blm, here("final_results", "scenario_4_combined", "sprob_combined_20_6hr_blm0.000001.tif"), overwrite = TRUE)
```

#### RUN 7: Combined - 30% target, Boundary Penalty = 0.000001 
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_30_blm <- problem(cost_combined, features = features_stack_4) %>%
  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_combined_30_blm <- solve(prob_combined_30_blm)

## Plot the solution to preview
plot(sprob_combined_30_blm,  main = c("combined- 30% Targets"))

## Save as tif
writeRaster(sprob_combined_30_blm, here("final_results", "scenario_4_combined", "sprob_combined_30_6hr_blm0.000001.tif"), overwrite = TRUE)
```

#### RUN 8: Combined - 50% target, Boundary Penalty = 0.000001 
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_50_blm <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.5) %>%
  add_locked_in_constraints(exist_mpas) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit) %>% 
  add_boundary_penalties(bp)

## Solve problem
sprob_combined_50_blm <- solve(prob_combined_50_blm)

## Plot the solution to preview
plot(sprob_combined_50_blm,  main = c("combined- 50% Targets"))

## Save as tif
writeRaster(sprob_combined_50_blm, here("final_results", "scenario_4_combined", "sprob_combined_50_6hr_blm0.000001.tif"), overwrite = TRUE)
```


### **STEP 4: Run Models with MPA sensitivity analysis**

#### Run 9: Combined - 20% target, boundary = 0 and no MPAs locked in
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_20_nompa <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.2) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit)

## Solve problem
sprob_combined_20_nompa <- solve(prob_combined_20_nompa)

## Plot the solution to see what it looks like
plot(sprob_combined_20_nompa,  main = c("combined- 20% Targets without mpas"))

## Save as tif
writeRaster(sprob_combined_20_nompa, here("final_results", "sensitivity", "mpas", "sprob_combined_20_6hr_blm0_nompa.tif"), overwrite = TRUE)
```

#### Run 10: Combined - 20% target, boundary = 0.000001 and no MPAs locked in
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_20_blm_nompa <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.2) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit) %>% 
  add_boundary_penalties(bp)

## Solve problem
sprob_combined_20_blm_nompa <- solve(prob_combined_20_blm_nompa)

## Plot the solution to see what it looks like
plot(sprob_combined_20_blm_nompa,  main = c("combined- 20% Targets without MPAS and with boundary penalty"))

## Save as tif
writeRaster(sprob_combined_20_blm_nompa, here("final_results", "sensitivity", "mpas", "sprob_combined_20_6hr_blm0.000001_nompa.tif"), overwrite = TRUE)
```


### **STEP 5: Run Models with KBA sensitivity analysis**

#### Run 11: Combined - 20% target, boundary = 0 and KBAs & MPAs locked in
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_20_mpa_kba <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.2) %>%
  add_locked_in_constraints(exist_mpas) %>% 
  add_locked_in_constraints(exist_kbas) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit)

## Solve problem
sprob_combined_20_mpa_kba  <- solve(prob_combined_20_mpa_kba )

## Plot the solution to see what it looks like
plot(sprob_combined_20_mpa_kba ,  main = c("combined- 20% Targets (with KBA and MPAs locked in)"))

#Save plot
writeRaster(sprob_combined_20_mpa_kba , here("final_results", "sensitivity", "key_biodiversity_areas", "sprob_combined_20_6hr_blm0_mpas_kbas.tif"), overwrite = TRUE)
```


#### Run 12: Combined - 20% target, boundary = 0 and KBAs locked in and MPA's not locked in
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_20_kba <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.2) %>%
  add_locked_in_constraints(exist_kbas) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit) 

## Solve problem
sprob_combined_20_kba  <- solve(prob_combined_20_kba )

## Plot the solution to see what it looks like
plot(sprob_combined_20_kba ,  main = c("combined- 20% Targets (with kbas locked in)"))

## Save as tiff
writeRaster(sprob_combined_20_kba , here("final_results", "sensitivity", "key_biodiversity_areas", "sprob_combined_20_6hr_blm0_kbas.tif"), overwrite = TRUE)
```

### **STEP 6: Run Models with Aggregation site sensitivity analysis**

#### Run 13: Combined - 20% target, boundary = 0 and Aggregation Sites & MPAS locked in
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_20_agg_mpa <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.2) %>%
  add_locked_in_constraints(exist_aggs) %>% 
  add_locked_in_constraints(exist_mpas) %>% 
  add_gurobi_solver(gap = gap, time_limit = time_limit) 

## Solve problem
sprob_combined_20_agg_mpa <- solve(prob_combined_20_agg_mpa)

## Plot the solution to see what it looks like
plot(sprob_combined_20_agg_mpa,  main = c("combined- 20% Targets (with aggregation & MPA locked in)"))

## Save as tif
writeRaster(sprob_combined_20_agg_mpa, here("final_results", "sensitivity", "aggregation_sites", "sprob_combined_20_6hr_blm0_agg_mpa.tif"), overwrite = TRUE)
```

#### Run 14: Combined - 20% target, boundary = 0 and aggregation sites
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_20_agg <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.2) %>%
  add_locked_in_constraints(exist_aggs) %>% 
  add_gurobi_solver(gap = gap, time_limit = time_limit) 

## Solve problem
sprob_combined_20_agg  <- solve(prob_combined_20_agg )

## Plot the solution to see what it looks like
plot(sprob_combined_20_agg ,  main = c("Combined- 20% Targets (with aggregation sites)"))

## Save as tif
writeRaster(sprob_combined_20_agg , here("final_results", "sensitivity", "aggregation_sites", "sprob_combined_20_6hr_blm0_aggsites.tif"), overwrite = TRUE)
```

#### Run 15: Combined - 20% target, boundary = 0 with aggregation sites, MPAs and KBAs locked in
```{r}
## Define the problem with the cost, conservation features, targets, locked in and locked out constraints using the minimum set objective
prob_combined_20_agg <- problem(cost_combined, features = features_stack_4) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.2) %>%
  add_locked_in_constraints(exist_aggs) %>% 
  add_locked_in_constraints(exist_mpas) %>% 
  add_locked_in_constraints(exist_kbas) %>%
  add_gurobi_solver(gap = gap, time_limit = time_limit) 

## Solve problem
sprob_combined_20_agg  <- solve(prob_combined_20_agg )

## Plot the solution to see what it looks like
plot(sprob_combined_20_agg ,  main = c("Combined- 20% Targets (with aggregation sites)"))

## Save as tif
writeRaster(sprob_combined_20_agg , here("final_results", "sensitivity", "sprob_combined_20_6hr_blm0_agg_mpa_kba.tif"), overwrite = TRUE)
```