3a_Linking_HRA_to_REC_habitat_2016_scenario_2025.Rmd

---
title: "Converting the outputs from HRA to inputs in the Recreation model using the habitat maps from 2016 - Scenario 2025"
author: "Jade Delevaux"
created: June 2022
output:
  html_document:
    df_print: paged
  pdf_document: default
---

To run this script, install the packages below

```{r Install packages}
library(sf)
library(foreign)
library(raster)
library(rgdal)
library(dplyr)
library(exactextractr) # quick zonal stat
library(ggplot2)
library(gridExtra)
library(cowplot)
library(rgeos)
library(sp)
library(ggpubr)
library(tigris)
library(maptools)
```

Note: The working directory (WD) needs to reflect your local machine. Use the function "Find and replace" in R studio to replace the root of the working directory by your local path (e.g., "E:/GreenFin") to apply this script.


# A/ Create shapefiles of habitat risk maps (Low and Medium risk maps) based on the 2025 Human stressors:

HRA outputs a risk raster map for each habitat (corals, seagrass, and mangroves), which classifies the level of risk as: no risk (0), low risk (1), medium risk (2), and high risk (3). We created rasters representing each habitat under each risk level. We need to convert those rasters to  shapefiles for the low and medium risk levels for each habitat type.

To import the HRA risk raster maps by habitat, point towards your working directory towards the folder with the HRA 2025 outputs. 

## A.1/ Coral habitat 

These section of codes will create two polygon shapefiles representing the coral habitat portions under no/low risk and medium risk.

First we convert the output coral risk map from HRA into a dissolved polygon and clip it by the footprint of the coral habitat shapefile
```{r Convert HRA raster to polygon}
raster_risk <-raster("E:/GreenFin/01_hra_model/5_hra_outputs/2016_habitats/2025/outputs/RECLASS_RISK_corals.tif")
r.to.poly<-rasterToPolygons(raster_risk, dissolve = T)
habitat_shp <- st_read("E:/GreenFin/01_hra_model/1_habitat_maps/2016/corals.shp")
int.r <-raster::intersect(r.to.poly, habitat_shp) 
writeOGR(obj=int.r, dsn="E:/GreenFin/03_recreation/hra_data/2016_habitats", layer="hra_output_corals_habitat_risk_2025", driver="ESRI Shapefile", overwrite_layer =  T)
```

Then we need to split the newly created polygon by no/low and medium risks maps and export those two maps into our Gdrive in the REC model folder. First we import the habitat risk polygon created above as a shapefile and set "risk" as a factor. THen, we create a map which only represents the no and low risk areas (i.e, the values = 0 & 1). We do the same process for the medium risk coral habitat (i.e., the value = 2):
```{r Split the polygon by each risk level}
habitat_risk <- st_read("E:/GreenFin/03_recreation/hra_data/2016_habitats/hra_output_corals_habitat_risk_2025.shp")
habitat_risk$RECLASS <- as.factor(habitat_risk$RECLASS)
target <- c('0','1')
habitat_risk_low <- habitat_risk[habitat_risk$RECLASS %in% target,]
habitat_risk_low_sp <- as(habitat_risk_low, "Spatial")
class(habitat_risk_low_sp)
writeOGR(obj=habitat_risk_low_sp, dsn="E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level", layer="hra_output_corals_habitat_risk_2025_low", driver="ESRI Shapefile", overwrite_layer =  T)
target <- c('2')
habitat_risk_med <- habitat_risk[habitat_risk$RECLASS %in% target,]
habitat_risk_med_sp <- as(habitat_risk_med, "Spatial")
class(habitat_risk_med_sp)
writeOGR(obj=habitat_risk_med_sp, dsn="E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level", layer="hra_output_corals_habitat_risk_2025_med", driver="ESRI Shapefile", overwrite_layer =  T)
```

## A.2/ Seagrass habitat

These section of codes will create two polygon shapefiles representing the seagrass habitat portions under no/low risk and medium risk.

First we convert the output seagrass risk map from HRA into a dissolved polygon and clip it by the footprint of the seagrass habitat shapefile
```{r Convert HRA raster to polygon}
raster_risk <-raster("E:/GreenFin/01_hra_model/5_hra_outputs/2016_habitats/2025/outputs/RECLASS_RISK_seagrass.tif")
r.to.poly<-rasterToPolygons(raster_risk, dissolve = T)
habitat_shp <- st_read("E:/GreenFin/01_hra_model/1_habitat_maps/2016/seagrass.shp")
int.r <-raster::intersect(r.to.poly, habitat_shp) 
writeOGR(obj=int.r, dsn="E:/GreenFin/03_recreation/hra_data/2016_habitats", layer="hra_output_seagrass_habitat_risk_2025", driver="ESRI Shapefile", overwrite_layer =  T)
```

Then we need to split the newly created polygon by no/low and medium risks maps and export those two maps into our Gdrive in the REC model folder. First we import the habitat risk polygon created above as a shapefile and set "risk" as a factor. THen, we create a map which only represents the no and low risk areas (i.e, the values = 0 & 1). We do the same process for the medium risk coral habitat (i.e., the value = 2):
```{r Split the polygon by each risk level}
habitat_risk <- st_read("E:/GreenFin/03_recreation/hra_data/2016_habitats/hra_output_seagrass_habitat_risk_2025.shp")
habitat_risk$RECLASS <- as.factor(habitat_risk$RECLASS)
target <- c('0','1')
habitat_risk_low <- habitat_risk[habitat_risk$RECLASS %in% target,]
habitat_risk_low_sp <- as(habitat_risk_low, "Spatial")
class(habitat_risk_low_sp)
writeOGR(obj=habitat_risk_low_sp, dsn="E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level", layer="hra_output_seagrass_habitat_risk_2025_low", driver="ESRI Shapefile", overwrite_layer =  T)
target <- c('2')
habitat_risk_med <- habitat_risk[habitat_risk$RECLASS %in% target,]
habitat_risk_med_sp <- as(habitat_risk_med, "Spatial")
class(habitat_risk_med_sp)
writeOGR(obj=habitat_risk_med_sp, dsn="E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level", layer="hra_output_seagrass_habitat_risk_2025_med", driver="ESRI Shapefile", overwrite_layer =  T)
```

## A.3/ Mangroves habitat

These section of codes will create two polygon shapefiles representing the mangroves habitat portions under no/low risk and medium risk.

First we convert the output mangroves risk map from HRA into a dissolved polygon and clip it by the footprint of the mangroves habitat shapefile
```{r Convert HRA raster to polygon}
raster_risk <-raster("E:/GreenFin/01_hra_model/5_hra_outputs/2016_habitats/2025/outputs/RECLASS_RISK_mangroves.tif")
r.to.poly<-rasterToPolygons(raster_risk, dissolve = T)
habitat_shp <- st_read("E:/GreenFin/01_hra_model/1_habitat_maps/2016/mangroves.shp")
int.r <-raster::intersect(r.to.poly, habitat_shp) 
writeOGR(obj=int.r, dsn="E:/GreenFin/03_recreation/hra_data/2016_habitats", layer="hra_output_mangroves_habitat_risk_2025", driver="ESRI Shapefile", overwrite_layer =  T)
```

Then we need to split the newly created polygon by no/low and medium risks maps and export those two maps into our Gdrive in the REC model folder. First we import the habitat risk polygon created above as a shapefile and set "risk" as a factor. THen, we create a map which only represents the no and low risk areas (i.e, the values = 0 & 1). We do the same process for the medium risk coral habitat (i.e., the value = 2):
```{r Split the polygon by each risk level}
habitat_risk <- st_read("E:/GreenFin/03_recreation/hra_data/2016_habitats/hra_output_mangroves_habitat_risk_2025.shp")
habitat_risk$RECLASS <- as.factor(habitat_risk$RECLASS)
target <- c('0','1')
habitat_risk_low <- habitat_risk[habitat_risk$RECLASS %in% target,]
habitat_risk_low_sp <- as(habitat_risk_low, "Spatial")
class(habitat_risk_low_sp)
writeOGR(obj=habitat_risk_low_sp, dsn="E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level", layer="hra_output_mangroves_habitat_risk_2025_low", driver="ESRI Shapefile", overwrite_layer =  T)
target <- c('2')
habitat_risk_med <- habitat_risk[habitat_risk$RECLASS %in% target,]
habitat_risk_med_sp <- as(habitat_risk_med, "Spatial")
class(habitat_risk_med_sp)
writeOGR(obj=habitat_risk_med_sp, dsn="E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level", layer="hra_output_mangroves_habitat_risk_2025_med", driver="ESRI Shapefile", overwrite_layer =  T)
```

Clear the working environment 
```{r Clear the workspace}
rm(list = ls(all.names=TRUE))
```


# B/ Calculate the functional habitat per habitat for the human stressors from 2025:
To link HRA tothe REC model, we need to translate habitat risk into functional habitat per hexagon using the following assumptions: 
High risk habitat -> 0% functional = no services
Medium risk habitat ->50% functional = some services 
Low risk habitat-> 100% functional = all delivery

First import the Recreation model AOI and reorder the attribute names and create predictor dataframe with the ID key (which i created to avoid the multiple PID with the same number from the REC gridded AOI shapefile which i intersected with the coastal planning units). We will add the habitat area from each habitat, under each risk level to eventually calculate the functional habitat.

```{r Create the REC database using the gridded AOI}
rec_aoi <- st_read("E:/GreenFin/03_recreation/_base/rec_grid_hex_w_cpu.shp")
rec_aoi <- rec_aoi %>%
  dplyr::select(id, tot_m2 = area_m2)
predictors <- rec_aoi
```

## B.1/ Calculate the functional habitat for the Coral habitat from 2016 and human stressors from 2025
First, we calculate the area (m2) of coral habitat under "no/low" risk within each hexagon from the REC model (i.e., PID).
```{r Calculate the area of habitat under low/no risk by PID}
poly_layer <- st_read("E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level/hra_output_corals_habitat_risk_2025_low.shp")
i_layer <- st_intersection(rec_aoi, poly_layer) 
i_layer$m2 <- unclass(st_area(i_layer)) # m2
pred_areas <- i_layer %>%
  st_set_geometry(NULL) %>%
  group_by(id) %>%
  summarise(pred_area = sum(m2)) ## 
```

We assume that under low risk, the habitat remains fully (100%) functional
```{r Apply the low risk 100% functional assumption}
preds_rec <- predictors %>%
  left_join(pred_areas, by = "id") %>%
  mutate(pct_pred = if_else(is.na(pred_area), 0, pred_area/tot_m2)) %>%
  dplyr::select(-pred_area)
predictors$pctcol <- preds_rec$pct_pred*100 
summary(predictors)
```

Second, we calculate the area (m2) of coral habitat under "medium" risk within each hexagon from the REC model (i.e., PID). 
```{r Calculate the area of habitat under medium risk by PID}
poly_layer <- st_read("E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level/hra_output_corals_habitat_risk_2025_med.shp")
i_layer <- st_intersection(rec_aoi, poly_layer) 
i_layer$m2 <- unclass(st_area(i_layer)) # m2
pred_areas <- i_layer %>%
  st_set_geometry(NULL) %>%
  group_by(id) %>%
  summarise(pred_area = sum(m2)) ## 
```

We assume that under medium risk, the habitat is 50% functional, so we divide those areas by 2.
```{r Apply the medium risk 50% functional assumption}
preds_rec <- predictors %>%
  left_join(pred_areas, by = "id") %>%
  mutate(pct_pred = if_else(is.na(pred_area), 0, pred_area/tot_m2)) %>%
  dplyr::select(-pred_area)
predictors$pctcom <- (preds_rec$pct_pred*100)/2 
summary(predictors)
```

Last, we combine the functional habitat areas from the "no/low" and "medium" risk to obtain the total functional habitat by hexagon: 
```{r Calculate the total functional habitat}
predictors$pctcor <- predictors$pctcom + predictors$pctcol
summary(predictors$pctcor)
```

## B.2/ Calculate the functional habitat for the Seagrass habitat from 2016 and human stressors from 2025

First, we calculate the area (m2) of seagrass habitat under "no/low" risk within each hexagon from the REC model (i.e., PID).
```{r Calculate the area of habitat under low/no risk by PID}
poly_layer <- st_read("E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level/hra_output_seagrass_habitat_risk_2025_low.shp")
i_layer <- st_intersection(rec_aoi, poly_layer) 
i_layer$m2 <- unclass(st_area(i_layer)) # m2
pred_areas <- i_layer %>%
  st_set_geometry(NULL) %>%
  group_by(id) %>%
  summarise(pred_area = sum(m2)) ## 
```

We assume that under low risk, the habitat remains fully (100%) functional
```{r Apply the low risk 100% functional assumption}
preds_rec <- predictors %>%
  left_join(pred_areas, by = "id") %>%
  mutate(pct_pred = if_else(is.na(pred_area), 0, pred_area/tot_m2)) %>%
  dplyr::select(-pred_area)
predictors$pctsgl <- preds_rec$pct_pred*100 
summary(predictors)
```

Second, we calculate the area (m2) of seagrass habitat under "medium" risk within each hexagon from the REC model (i.e., PID). 
```{r Calculate the area of habitat under medium risk by PID}
poly_layer <- st_read("E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level/hra_output_seagrass_habitat_risk_2025_med.shp")
i_layer <- st_intersection(rec_aoi, poly_layer) 
i_layer$m2 <- unclass(st_area(i_layer)) # m2
pred_areas <- i_layer %>%
  st_set_geometry(NULL) %>%
  group_by(id) %>%
  summarise(pred_area = sum(m2)) ## 
```

We assume that under medium risk, the habitat is 50% functional, so we divide those areas by 2.
```{r Apply the medium risk 50% functional assumption}
preds_rec <- predictors %>%
  left_join(pred_areas, by = "id") %>%
  mutate(pct_pred = if_else(is.na(pred_area), 0, pred_area/tot_m2)) %>%
  dplyr::select(-pred_area)
predictors$pctsgm <- (preds_rec$pct_pred*100)/2 
summary(predictors)
```

Last, we combine the functional habitat areas from the "no/low" and "medium" risk to obtain the total functional habitat by hexagon: 
```{r Calculate the total functional habitat}
predictors$pctsgr <- predictors$pctsgm + predictors$pctsgl
summary(predictors$pctsgr)
```

## B.3/ Calculate the functional habitat for the Mangroves habitat from 2016 and human stressors from 2025

First, we calculate the area (m2) of mangroves habitat under "no/low" risk within each hexagon from the REC model (i.e., PID).
```{r Calculate the area of habitat under low/no risk by PID}
poly_layer <- st_read("E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level/hra_output_mangroves_habitat_risk_2025_low.shp")
i_layer <- st_intersection(rec_aoi, poly_layer) 
i_layer$m2 <- unclass(st_area(i_layer)) # m2
pred_areas <- i_layer %>%
  st_set_geometry(NULL) %>%
  group_by(id) %>%
  summarise(pred_area = sum(m2)) ## 
```

We assume that under low risk, the habitat remains fully (100%) functional
```{r Apply the low risk 100% functional assumption}
preds_rec <- predictors %>%
  left_join(pred_areas, by = "id") %>%
  mutate(pct_pred = if_else(is.na(pred_area), 0, pred_area/tot_m2)) %>%
  dplyr::select(-pred_area)
predictors$pctmgl <- preds_rec$pct_pred*100 
summary(predictors)
```

Second, we calculate the area (m2) of seagrass habitat under "medium" risk within each hexagon from the REC model (i.e., PID). 
```{r Calculate the area of habitat under medium risk by PID}
poly_layer <- st_read("E:/GreenFin/03_recreation/hra_data/2016_habitats/one_risk_level/hra_output_mangroves_habitat_risk_2025_med.shp")
i_layer <- st_intersection(rec_aoi, poly_layer) 
i_layer$m2 <- unclass(st_area(i_layer)) # m2
pred_areas <- i_layer %>%
  st_set_geometry(NULL) %>%
  group_by(id) %>%
  summarise(pred_area = sum(m2)) ## 
```

We assume that under medium risk, the habitat is 50% functional, so we divide those areas by 2.
```{r Apply the medium risk 50% functional assumption}
preds_rec <- predictors %>%
  left_join(pred_areas, by = "id") %>%
  mutate(pct_pred = if_else(is.na(pred_area), 0, pred_area/tot_m2)) %>%
  dplyr::select(-pred_area)
predictors$pctmgm <- (preds_rec$pct_pred*100)/2 
summary(predictors)
```

Last, we combine the functional habitat areas from the "no/low" and "medium" risk to obtain the total functional habitat by hexagon: 
```{r Calculate the total functional habitat}
predictors$pctmgv <- predictors$pctmgm + predictors$pctmgl
summary(predictors$pctmgv)
```

Before exporting it, we remove the columns with the areas of low and medium risks because we do not need them for the REC model. 
```{r Clean the database}
names(predictors)
predictors <- predictors[-c(4:5,7:8,10:11)]
names(predictors)
```

Export it as a shapefile and a csv file 
```{r Export as shp and csv}
predictors_sp <- as(predictors, "Spatial")
class(predictors_sp)
writeOGR(obj=predictors_sp, dsn="E:/GreenFin/03_recreation/predictors/2016_habitats/shp", layer="preds_functional_habitats_2025", driver="ESRI Shapefile", overwrite_layer =  T)
write.csv(predictors %>% st_set_geometry(NULL), paste0("E:/GreenFin/03_recreation/predictors/2016_habitats/csv/preds_functional_habitats_2025", ".csv"))
```

Last, clear the working environment
```{r Clear the workspace}
rm(list = ls(all.names=TRUE))
```

# C/ Rasterize functional habitat predictors for human stressors 2025 

Import the polygon shapefile with the functional habitat for coral, seagrass nad mangroves in percent cover per hexagon.  
```{r Import the shp with the 2025 functional habitat values}
rec_predictors_poly <- st_read("E:/GreenFin/03_recreation/predictors/2016_habitats/shp/preds_functional_habitats_2025.shp")
names(rec_predictors_poly)
```

Rasterize the functional habitat percent cover by hexagon for each habitat type. 

Coral: 
```{r Rasterize 2025 coral functional habitat}
blank_raster <- raster(as(rec_predictors_poly, "Spatial"), ncol=6543, nrow=9750)
predictor_raster <- rasterize(as(rec_predictors_poly, "Spatial"), blank_raster, 'pctcor', progress = "text")
plot(predictor_raster)
writeRaster(predictor_raster, filename="E:/GreenFin/03_recreation/predictors/2016_habitats/corals_2025.tif", overwrite=TRUE)
```

Seagrass:  
```{r Rasterize 2025 seagrass functional habitat}
blank_raster <- raster(as(rec_predictors_poly, "Spatial"), ncol=6543, nrow=9750)
predictor_raster <- rasterize(as(rec_predictors_poly, "Spatial"), blank_raster, 'pctsgr', progress = "text")
plot(predictor_raster)
writeRaster(predictor_raster, filename="E:/GreenFin/03_recreation/predictors/2016_habitats/seagrass_2025.tif", overwrite=TRUE)
```

Mangroves:  
```{r Rasterize 2025 mangroves functional habitat}
blank_raster <- raster(as(rec_predictors_poly, "Spatial"), ncol=6543, nrow=9750)
predictor_raster <- rasterize(as(rec_predictors_poly, "Spatial"), blank_raster, 'pctmgv', progress = "text")
plot(predictor_raster)
writeRaster(predictor_raster, filename="E:/GreenFin/03_recreation/predictors/2016_habitats/mangroves_2025.tif", overwrite=TRUE)
```