set_up.Rmd

---
title: 'EDS Workshop: Prioritizr - Setup'
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(fasterize)
library(here)
library(tmap)
```

## Planning Units 

#### **Mozambique's EEZ**
```{r}
## Grab Mozambique's EEZ shapefile
mz_eez_sf <- sf::read_sf(here("raw_data", "mz_eez"))

## Check coordinate system
st_crs(mz_eez_sf)

## Use Moz EEZ to set extent for future rasters
mz_ext <- raster::extent(mz_eez_sf)

## Create raster with cell ids and clip to MZ EEZ raster
mz_eez_rast <- raster::raster(x=mz_ext, crs=crs(mz_eez_sf), res=10000)

## Assign cell values
values(mz_eez_rast) <- 1:ncell(mz_eez_rast) 

## Create data frame to check it out
mz_eez_rast_df <- rasterToPoints(mz_eez_rast) %>%
  as.data.frame()

## Plot to make sure it looks good
plot(mz_eez_rast)

## Mask it to only include the EEZ
mz_eez_rast <- mask(mz_eez_rast, mz_eez_sf)
plot(mz_eez_rast)

## Let's save this cell id raster as tif file - hashtagged this out for now since the tif is too big
writeRaster(mz_eez_rast, here('cleaned_data/mz_eez_rast.tif'), overwrite = TRUE)
```

#### **Existing MPAs**
```{r}
## Grab protected area shapefile from computer
exist_mpa_path <- 'G:/group_project/data/existing_mpas'
mpa_shape <- list.files(exist_mpa_path, full.names = TRUE)

## Read in the shapefile as a simple feature
mpa_sf <- sf::read_sf(mpa_shape[7]) 

## Check the CRS of the simple feature to make sure it matches the Mozambique Raster 
st_crs(mpa_sf)

## CRS matches so lets make a raster using the mz_rast_id to set the extent
mpa_rast <- fasterize::fasterize(mpa_sf, mz_eez_rast)

## Double check the CRS again 
crs(mpa_rast)

## Plot to make sure it looks good, note this shapefile includes all terrestrial protected areas as well so we need to clip it to just include those in the water (the EEZ)
plot(mpa_rast)

## Mask the raster to the EEZ 
mpa_rast <- mask(mpa_rast, mz_eez_sf) 

## Double check CRS of new raster & plot to make sure it looks good
crs(mpa_rast)
plot(mpa_rast)

## Let's save this cell id raster as tif file - hashtagged this out for now since the tif is too big
writeRaster(mpa_rast, here('cleaned_data/mpa_rast.tif'), overwrite = TRUE)
```


## Conservation Features

### **Species**

#### **IUCN Species Distribution**
To download the IUCN Chronrichthyes spatial data locally, follow this [*link*](https://www.iucnredlist.org/resources/spatial-data-download).
```{r}
### grab iucn data from local computer
iucn_path <- 'G:/group_project/data/updated_iucn'
iucn_shape <- list.files(iucn_path, full.names = TRUE) 

iucn_sf <- sf::read_sf(iucn_shape[9])

### create a dataframe to look at the different variables and remove geometry to speed it up
iucn_df <- iucn_sf %>%
  as.data.frame() %>% 
  select(-geometry)
```

```{r}
# create raster function by species name
create_species_rast <- function(iucn_sf, species_name){
  
  species_df <- iucn_sf %>%
    filter(binomial == species_name) %>% #filter by species name
    st_transform(., crs = st_crs(mz_eez_sf)) #change crs to match mz_eez_sf

  species_rast <- fasterize::fasterize(species_df, mz_eez_rast) %>% #create a species presence raster
    mask(mz_eez_sf) #make everything outside the EEZ a NA

  return(species_rast)
  
}
```

```{r}
# Mobula kuhlii - shortfin devil ray
devil_ray <- create_species_rast(iucn_sf, "Mobula kuhlii")
plot(devil_ray)

# Looks good lets save as a tif
writeRaster(devil_ray, here('cleaned_data/devil_ray.tif'), overwrite = TRUE)

#Pristis pristis - largetooth sawfish
sawfish <- create_species_rast(iucn_sf, "Pristis pristis")
plot(sawfish)

# Looks good lets save as a tif
writeRaster(sawfish, here('cleaned_data/sawfish.tif'), overwrite = TRUE)

#Pristis pristis - largetooth sawfish
scalloped_hammerhead <- create_species_rast(iucn_sf, "Sphyrna lewini")
plot(scalloped_hammerhead)

# Looks good lets save as a tif
writeRaster(scalloped_hammerhead, here('cleaned_data/scalloped_hammerhead.tif'), overwrite = TRUE)
```

### **Habitats**

#### **SEAMOUNTS**
```{r}
## Grab seamount shapefile from the local computer
seamount_path <- 'G:/group_project/Data/Habitats/seamounts'
seamount_shape <- list.files(seamount_path, full.names = TRUE)

## Load shapefile as simple feature
seamount_sf <- sf::read_sf(seamount_shape[6]) 

## Check the CRS of the seamount and compare to the mz rast
st_crs(mz_eez_rast)
st_crs(seamount_sf)

## Let's create a rsater using the mz_rast with resolution of 10000 and mask tot the MZ EEZ
seamounts <- fasterize::fasterize(seamount_sf, mz_eez_rast)%>%
  mask(mz_eez_sf)

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

# Looks good lets save as a tif
writeRaster(seamounts, here('cleaned_data/seamounts.tif'), overwrite = TRUE)
```

#### **KNOLLS**
```{r}
## Grab knolls shapefile from the local computer
knolls_path <- 'G:/group_project/Data/Habitats/knolls'
knolls_shape <- list.files(knolls_path, full.names = TRUE)

## Load shapefile as simple feature
knolls_sf <- sf::read_sf(knolls_shape[6]) 

## Check the CRS of the knolls sf
st_crs(knolls_sf)

## Let's create a rsater using the mz_rast with resolution of 10000 and mask tot the MZ EEZ
knolls <- fasterize::fasterize(knolls_sf, mz_eez_rast)%>%
  mask(mz_eez_sf)

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

# Looks good - lets save as a tif
writeRaster(knolls, here('cleaned_data/knolls.tif'), overwrite = TRUE)
```

## Cost Layer

#### **Industrial Fishing Pressure from Global Fishing Watch**

```{r}
#path to Global Fishing Watch fishing pressure data
gfw_path <- 'G:/GFW2'
gfw_files <- list.files(gfw_path, full.names = TRUE)

#create a loop to read in data as csv 
for (i in 1:length(gfw_files)) assign(gfw_files[i], read.csv(gfw_files[i]))

#create a function that takes the csv and changes into raster
fishing_rast <- function(csv, mz_eez_sf){
  
  ## Convert the data frame to a simple features object (sf) to create a geometry for each id
  fishing_sf <-  st_as_sf(csv, coords = c("cell_ll_lon", "cell_ll_lat"), crs= 4326) %>% #make sf
    st_transform(., crs = st_crs(mz_eez_sf))
  
  ## check the class of sf
  class(fishing_sf)
  
  ## change to a sp object
  fishing_sp <- as(fishing_sf, 'Spatial')
  
  ## check the class to make sure it worked
  class(fishing_sp)
  
  ## double check that this matches your RASTER file (which should be the same as the sf)
  crs(fishing_sp)
  crs(mz_eez_rast)
  
  ## make a raster 
  fishing_rast <- raster::rasterize(fishing_sp, mz_eez_rast, field = "mmsi_present")

return(fishing_rast)
  
}
```

#### **Fishing Presence from December 27-31, 2020**
```{r}
#fishing presence on 27 Dec 2020
fishing_27 <- fishing_rast(`G:/GFW2/2020-12-27.csv`, mz_eez_sf)
plot(fishing_27)

#to check what the data looks like
fishing_27_df <- fishing_27 %>% 
  as.data.frame() 

#fishing presence on 28 Dec 2020
fishing_28 <- fishing_rast(`G:/GFW2/2020-12-28.csv`, mz_eez_sf)
plot(fishing_28)

#fishing presence on 29 Dec 2020
fishing_29 <- fishing_rast(`G:/GFW2/2020-12-29.csv`, mz_eez_sf)
plot(fishing_29)

#fishing presence on 30 Dec 2020
fishing_30 <- fishing_rast(`G:/GFW2/2020-12-30.csv`, mz_eez_sf)
plot(fishing_30)

#fishing presence on 31 Dec 2020
fishing_31 <- fishing_rast(`G:/GFW2/2020-12-31.csv`, mz_eez_sf)
plot(fishing_31)

#create a stack of all fishing pressure to be used in Prioritzr
fishing_stack <- stack(fishing_27, fishing_28, fishing_29, fishing_30, fishing_31)
plot(fishing_stack)

```

```{r}
#adding all the rasters together to get a single raster
fishing_stack_sum <- calc(fishing_stack, fun = sum, na.rm = TRUE) %>% 
  mask(mz_eez_sf)#add a mask so only the values in the project area are shown

#let's plot it
plot(fishing_stack_sum)

#check it out just to double check 
fishing_sum_df <- fishing_stack_sum %>% 
  as.data.frame()

#Prioritzr does NOT like 0s, so we changed all values from 0 to 0.001
fishing_stack_sum[values(fishing_stack_sum == 0)] <- 0.001

#check it once last time!
fishing_sum_df <- fishing_stack_sum %>% 
  as.data.frame() 

## Looks good, let's save as a multilayer feature tif
writeRaster(fishing_stack_sum, filename=here("cleaned_data", "fishing_stack.tif"), overwrite=TRUE)
```



**For fun... creating a interactive map!**
```{r}
tmap_mode("view") # Set to interactive viewing
tm_shape(fishing_stack_sum) +
  tm_basemap("Stamen.Terrain")+
  tm_raster("layer",alpha=0.9, palette = "viridis", breaks = c(0,1,2,5,25), title = "Industrial Fishing")+
tm_shape(mz_eez_sf)+
  tm_borders("black")
```