Final_SDM_project_workflow.qmd

---
title: "Team_SDM_OHW24"
format: html
editor: visual
---

```{r message = FALSE, warning = FALSE}
knitr::opts_chunk$set(message = FALSE, warning = FALSE)

```

# Maine Icons: A Species Distribution Model and Educational Tool to Highlight Gulf of Maine Creatures

**Authors: Hallie Arno, Mugdha Chiplunkar, Jessie Maguire, Camille Ross, Ben Tupper**

------------------------------------------------------------------------

## Project Overview

This project was developed during [OceanHackWeek 2024](https://oceanhackweek.org/ohw24/) with the goal of creating an interactive tool using `R` to show the predicted distribution of some iconic Gulf of Maine species.

------------------------------------------------------------------------

## Project goals

To create a species distribution model for various marine species native to the Gulf of Maine that can be modified by switching out the year, RCP climate scenario, and species of interest. We wanted to turn our results into an interactive app that can be used for educational purposes to visualize the maps of forecasted distribution shiftsunder several scenarios.

------------------------------------------------------------------------

## What is a species distribution model?

A species distribution model (SDM) is a tool used to predict the geographic distribution of a species based on environmental conditions and species occurrence data. By analyzing the relationship between known locations where a species is found and the environmental factors at those locations (like temperature, precipitation, or habitat type), SDMs can estimate the suitability of different areas for that species. These models are commonly used in ecology and conservation to understand species ranges, assess the impact of climate change, and guide habitat management efforts.

For our project, we decided to focus on the Gulf of Maine as our geographic area, since it is warming faster than 99% of our world's oceans.

------------------------------------------------------------------------

## Data

**Biological**

We obtained species occurrence data from the [Ocean Biodiversity Information System (OBIS)](https://obis.org/) via the `robis` package. Since absences were not recorded in this data, we generated 1000 random background points (pseudo-absences).

**Environmental**

We obtained environmental data from [Bio-ORACLE](https://www.bio-oracle.org/) and used the `stars` package to wrangle the raster files.

We used predictions of future environmental conditions from the Atmosphere/Ocean General Circulation Model (AOGCM).

------------------------------------------------------------------------

## References

-   OBIS (2023) Ocean Biodiversity Information System. Intergovernmental Oceanographic Commission of UNESCO. www.obis.org. Accessed: 2024-08-27

-   Assis, J., Fernández Bejarano, S.J., Salazar, V.W., Schepers, L., Gouvêa, L., Fragkopoulou, E., Leclercq, F., Vanhoorne, B., Tyberghein, L., Serrão, E.A., Verbruggen, H., De Clerck, O. (2024) Bio-ORACLE v3.0. Pushing marine data layers to the CMIP6 Earth system models of climate change research. Global Ecology and Biogeography. DOI:

------------------------------------------------------------------------

## Project Workflow

Entire workflow of SDM using Maxnet package in R

## Installing libraries

```{r}
  #required libraries for data retrieval, cleaning, and modeling
library(tidyverse)
library(yaml)
library(raster)
library(robis)
library(stars)
library(dismo)
library(terra)
library(maxnet)
library(ggplot2)
library(leaflet)
library(leafem)
```

Create YAML file to define all variables used in functions (spatial and temporal bounds, species names, and env layers)

```{r}
vars <- read_yaml("config2.yaml")
```

## Get species data from Obis using function

```{r}
  #Input: Latin name of a marine species, Output: data frame with occurrence points
get_species_data <- function(spec) {
  
    species_data <- robis::occurrence(spec)
    withDates <- species_data %>%
      separate(eventDate, into = c("Year", "Month"), sep = "-") %>%
      filter(!is.na(Year)) %>%
      filter(grepl("^\\d{4}$", as.character(Year)))
    
    filtered_data <- subset(withDates, date_year >= vars$start_year & date_year <= vars$end_year)
    
    bounded_data <- filtered_data %>%
      filter(decimalLatitude >= vars$latmin & decimalLatitude <= vars$latmax &
               decimalLongitude >= vars$lonmin & decimalLongitude <= vars$lonmax) %>% 
      dplyr::select(datasetName, decimalLatitude, decimalLongitude, Year, Month, individualCount, vernacularName)
    
    obs_sf <- bounded_data %>% 
    sf::st_as_sf(
      coords = c("decimalLongitude", "decimalLatitude"),
      crs = st_crs(4326))
    return(obs_sf)
  
}

```

## Obtain and extract environmental data from BioOracle using functions

```{r}
  #Input: variables, output: raster stack
get_enviro_data <- function(envvars) {
  #layercodes <- var
  #dir = "ohw24_proj_sdm_us"
  env <- sdmpredictors::load_layers(envvars, equalarea = FALSE, rasterstack = TRUE)
  #Crop
  env <- st_as_stars(env)
  extent <- st_bbox(c(xmin = vars$lonmin, xmax = vars$lonmax, ymin = vars$latmin, ymax = vars$latmax), crs = st_crs(env))
  rc <- st_crop(x = env, y = extent)
  return(rc)
}

  #extract specific environmental covariates 
  #ph <- get_enviro_data("BO_ph")
extractEnvData <- function(rasterStack, points) {
  env.stars <- terra::split(rasterStack)
  spec.env <- stars::st_extract(env.stars, sf::st_coordinates(points))
  na.omit(spec.env)
  return(spec.env)
  
}
```

## Creating pseudo-absence points and cropping it to environmental layer

```{r}
  #generate pseudo-absence points
getNegativePoints <- function(croppedRaster, nsamp = 1000) {
  bbox_sf <- st_as_sfc(st_bbox(c(xmin = vars$lonmin, xmax = vars$lonmax, ymin = vars$latmin, ymax = vars$latmax), crs = st_crs(croppedRaster)))
  set.seed(42)  # For reproducibility
  random_points <- st_sample(bbox_sf, size = nsamp)
  
  # Convert points to a data frame and then to an sf object
  random_points_sf <- st_as_sf(as.data.frame(st_coordinates(random_points)), coords = c("X", "Y"), crs = st_crs(croppedRaster))
  
  # Crop the points to the extent of the environmental layer
  cropped_points <- st_intersection(random_points_sf, st_as_sf(croppedRaster, as_points = FALSE, merge = TRUE))
  return(cropped_points)
}
```

## Call functions to prep the data for model

```{r}
layers <- c(vars$envVars)
envRasterStack <- get_enviro_data(layers[1:5])
speciesPoints <- get_species_data(vars$species[1])
absPoints <- getNegativePoints(envRasterStack) 

pres <- extractEnvData(envRasterStack, speciesPoints) |> mutate(pa=1)

abs <- extractEnvData(envRasterStack, absPoints) |> mutate(pa=0)

allData <- rbind(pres, abs)

head(allData)
```

## Fit the model

```{r}
  #sdm.model \<- maxnet::maxnet(presence_absence_df, environmental_df)}
  #responses \<- plot(sdm.model, type = "cloglog")
presence_absence_df <- allData %>%
  dplyr::select(pa)

environmental_df <- allData %>%
  dplyr::select(-c(pa))

  # Ensure that 'Presence' column is extracted as a numeric vector
presence_absence_vector <- presence_absence_df$pa

  # Fit the MaxEnt model
sdm.model <- maxnet::maxnet(p = presence_absence_vector, data = environmental_df)

  # Plot the response curves
responses <- plot(sdm.model, type = "logistic")

```

## Process future environmental scenario data to predict species distribution

```{r}
future_layers <- sdmpredictors::list_layers_future(marine = TRUE) %>% 
  filter(current_layer_code %in% c(vars$envVars)) %>% 
  filter(year == "2050") %>% 
  filter(scenario == "RCP26") %>% 
  filter(model == "AOGCM")
future_layers_list <- future_layers$layer_code

  #concatenate and stack all environmental raster layers 
BO21_tempmean_bdmax <- get_enviro_data(future_layers_list[1])
BO21_tempmean_ss <- get_enviro_data(future_layers_list[2])
BO22_chlomean_ss <- get_enviro_data(future_layers_list[3])
BO22_salinitymean_bdmax <- get_enviro_data(future_layers_list[4])
BO22_salinitymean_ss <- get_enviro_data(future_layers_list[5])
concat <- c(BO21_tempmean_bdmax, BO21_tempmean_ss, BO22_chlomean_ss, BO22_salinitymean_bdmax, BO22_salinitymean_ss)

names(concat) <- future_layers$current_layer_code

  #logistic model type
clamp <- TRUE      
type <- "logistic"

  # Predict species distribution within the cropped area
predicted <- predict(sdm.model, concat)

plot(predicted)
```

## Plot the SDM using leaflet

```{r}
m <- leaflet() %>% 
  addTiles() %>% 
  leafem::addStarsImage(predicted, 
                        colors = viridis::viridis(256), 
                        opacity = 0.8)

m
```

## Plotting SDM for alternate year and climate scenario

```{r}
  #plotting SDM for the year 2100 in climate scenario RCP85
future_layers_RCP85 <- sdmpredictors::list_layers_future(marine = TRUE) %>% 
  filter(current_layer_code %in% c(vars$envVars)) %>% 
  filter(year == "2100") %>% 
  filter(scenario == "RCP85") %>% 
  filter(model == "AOGCM")
future_layers_list_RCP85 <- future_layers_RCP85$layer_code

  #concatenate and stack all environmental raster layers 
BO21_tempmean_bdmax <- get_enviro_data(future_layers_list_RCP85[1])
BO21_tempmean_ss <- get_enviro_data(future_layers_list_RCP85[2])
BO22_chlomean_ss <- get_enviro_data(future_layers_list_RCP85[3])
BO22_salinitymean_bdmax <- get_enviro_data(future_layers_list_RCP85[4])
BO22_salinitymean_ss <- get_enviro_data(future_layers_list_RCP85[5])
concat <- c(BO21_tempmean_bdmax, BO21_tempmean_ss, BO22_chlomean_ss, BO22_salinitymean_bdmax, BO22_salinitymean_ss)

names(concat) <- future_layers_RCP85$current_layer_code

  #logistic model type
clamp <- TRUE      
type <- "logistic"

  # Predict species distribution within the cropped area
predicted_RCP85 <- predict(sdm.model, concat)

plot(predicted_RCP85)
```

## Map SDM for RCP85 in year 2100 in leaflet and visualize results

```{r}
library(leaflet)
m <- leaflet() %>% 
  addTiles() %>% 
  leafem::addStarsImage(predicted_RCP85, 
                        colors = viridis::viridis(256), 
                        opacity = 0.8)

m
```

## Interactive R Shiny website

We created an interactive website to select desired species, year, and climate scenario using drop-down menus to visualize the corresponding species distribution model!