HWI appendix.Rmd

---
title: "Appendix: Human-wildlife Interactions"
author: "Dayna K Weststrate, Aimee Chhen, Stefano Mezzini, Michael J Noonan"
output: html_document
date: "This document was created on May 8, 2023. It was last modified on `r format(Sys.Date(), '%B %d, %Y')`."
---

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

# Overview

This appendix details the steps we used to analyze our data, build our model, and produce our results.

```{r}
#data, visualization
library(readr)
library(lubridate)       #date formats
library(zoo)             #date format 'year-month'
library(dplyr)           #data wrangling
library(tidyr)           #data wrangling
library(ggplot2)
library(khroma)          #colour blind friendly palette
library(viridis)
library(mgcv)            #for gams
library(MuMIn)           #for model selection, AICc()
library(gridExtra)
library(ggh4x)           #to fill in facet wrap title
library(canadianmaps)    #to download a shapefile of BC
library(elevatr)         #to download digital elevation models
library(purrr)           #for functional programming (map_***(), etc.)
library(sp)              #spatial data, SpatialPoints()
library(sf)              #for spatial data
library(terra)           #for raster data
library(progress)        #for elevation, get_elev_raster()
library(stringi)
library(climatenaR)
```


```{r include = FALSE}
#import coordinates obtained from Google maps of the national parks
#refer to 'elevation and climate.R' for script
nationalparks_bc_coordinates <- read_csv("Data/HWI/nationalparks_bc_coordinates.csv")
#import elevation data
#nationalparks_bc_dem <- read.csv('data/ClimateNA_v731/nationalparks_bc_dem.csv')
# Import cleaned historical climate data
#refer to "Historical climate data" section in 'elevation and climate.R' for script
historical_climate_data <- readRDS("Data/HWI/ClimateNA_v731/climate data/historical_climate_data.rds")
# Import model
model <- readRDS("Data/HWI/model.rds")
```


# Data Preparation

Human-wildlife coexistence incidents (HWCI) data was obtained from the Government of Canada’s Open Government database <https://open.canada.ca/data/en/dataset/cc5ea139-c628-46dc-ac55-a5b3351b7fdf/resource/b2a9f7e4-7c49-471d-8337-0192c15dd52a?inner_span=True>. 
The Open Data Record dataset is complied by Parks Canada Agency (PCA) which contains human-wildlife coexistence incident recordings from 2010 to 2021 of thirty-five national parks and historic sites. 
Our project is interested in the National Parks located within British Columbia (BC).
The HWCI data contains 14 variables and 4 variables (Incident Date, Protected Heritage Area, Incident Type, and Species Common Name) were selected for this analysis.
There is a total of 9 Incident Types assigned in the dataset. Our project is interested in the 'Human Wildlife Interaction' aspect of Human-wildlife coexistence (HWC).


```{r}
# Import Parks Canada dataset
#data obtained on Feb. 3, 2023
PCA <- read_csv("Data/HWI/pca-human-wildlife-coexistence-animals-involved-detailed-records-2010-2021.csv")

#data carpentry
BC <- PCA[which(PCA$`Protected Heritage Area` %in% c("Glacier National Park of Canada",
                                                     "Kootenay National Park of Canada",
                                                     "Mount Revelstoke National Park of Canada",
                                                     "Pacific Rim National Park Reserve of Canada",
                                                     "Yoho National Park of Canada")),]
BC <- BC[which(BC$`Incident Type` == "Human Wildlife Interaction"),]

#rename column names to make it easier for coding
names(BC)[2] <- "incident_date"
names(BC)[4] <- "park"
names(BC)[5] <- "HWI"
names(BC)[6] <- "species"

#clean BC National Parks data, removing unknown species
BC <- BC[BC$species != "None",]
BC <- BC[BC$species != "Unknown",]
BC <- BC[BC$species != "Unknown bat",]
BC <- BC[BC$species != "Unknown bear",]
BC <- BC[BC$species != "Unknown bird",]
BC <- BC[BC$species != "Unknown deer",]
BC <- BC[BC$species != "Unknown sea lion",]

# Data preparation for climate data
#Note: climate data is in monthly intervals
BC$year <- lubridate::year(BC$incident_date)
#extract month from date then create a column for month, required for climate data
BC$month <- lubridate::month(BC$incident_date)
#combine year and month column into a single column
BC$year_month <- as.yearmon(paste(BC$year, BC$month), "%Y %m") 

# Create a new dataframe for analysis
#subset the dataset for HWIs grouped monthly to correspond with the climate data 
data <- aggregate(HWI ~ year_month + park, data = BC, FUN = "length")
data$year_month <- as.Date(data$year_month, format = "%Y-%m")
data <- relocate(data, HWI, .before = year_month)
```

# Historical climate data

Climate data (temperature, precipitation) was obtained from the open source ```ClimatenaR``` R package. 
Elevation data is required for the ```ClimatenaR``` package to work. 
The coordinates of each national parks were obtained from Google maps and the Digital Elevation Model (DEM) was obtained using ```get_elev_raster()``` function from the ```elevatr``` package.

```{r eval = FALSE}
#import a shapefile of British Columbia
bc_shape <- st_as_sf(PROV) %>%  # convert to spatial features (sf) object
  filter(PRENAME == 'British Columbia') %>% # filter to BC only
  st_geometry() # extract boundaries only

#import coordinates obtained from Google maps of the national parks
nationalparks_bc_coordinates <- read_csv("Data/HWI/ClimateNA_v731/nationalparks_bc_coordinates.csv")

#convert all telemetry dataset to spatial data points
nationalparks_location <- SpatialPoints(select(nationalparks_bc_coordinates, longitude, latitude))

ctmm::projection(nationalparks_location) <- '+proj=longlat' 
#do not load the the ctmm package for this or you will run into errors

#check locations
plot(bc_shape)
sp::plot(nationalparks_location, add = TRUE, col = 'red', pch = 19, cex = 0.5) 
#need sp:: in front of plot because function will try to use plot() from another package
```


## Elevation data

```{r eval = FALSE}
#obtain elevation data required for climate data
#import a Digital Elevation Model (DEM) for the region(s) of interest
dem <- get_elev_raster(locations = nationalparks_location,
                       z = 3,
                       clip = 'bbox',
                       expand = 0.1)

plot(dem)
plot(nationalparks_location, add = TRUE, col = 'red', pch = 17, cex = 1.75)

# write the csv
#(circumventing climatenaR's functions because we only need a specific number of locations)
nationalparks_bc_coordinates <- mutate(nationalparks_bc_coordinates,
                                       el = terra::extract(dem, nationalparks_location),
                                       ID1 = 1:n(),
                                       ID2 = ID1) %>%
  relocate(ID1:ID2)
#save elevation file into 'ClimateNA_v731' folder for downloading climate data (next step)
write.csv(nationalparks_bc_coordinates, file = 'Data/HWI/ClimateNA_v731/nationalparks_bc_dem.csv', row.names = FALSE) 

#check the csv
nationalparks_bc_dem <- read.csv('Data/HWI/ClimateNA_v731/nationalparks_bc_dem.csv') %>%
  head()
```

```{r}
#plot coordinates for visual
ggplot(nationalparks_bc_dem, aes(latitude, longitude, el)) +
  geom_tile()
```



## Download Historical climate data

```{r eval = FALSE}
# download historical climate data
for(y in 2010:2021) {
  cat('Downloading ', y, '...\n', sep = '') # to track progress
  histClimateNA(
    file = '/Data/HWI/ClimateNA_v731/nationalparks_bc_dem.csv', #elevation file, error, remove "park" column in the csv to fix
    dateR = as.character(y),
    tFrame = 'M', # monthly averages
    exe = '/Data/HWI/ClimateNA_v731/ClimateNA_v7.31.exe',
    outdir = '/Data/HWI/ClimateNA_v731/climate data/historical climate data') #create a output file where the data will be saved to within the 'ClimateNA_v731' folder
}
```

## Clean historical climate data

```{r eval = FALSE}
historical_climate_data <-
  # list all files, and import each of the CSVs
  map_dfr(
    list.files('Data/HWI/ClimateNA_v731/climate data/historical climate data', full.names = TRUE), #folder where the downloaded CSV of historical data is located
    \(.fname) {
      readr::read_csv(.fname, col_types = '?') %>%
        # add a column of the file name
        mutate(file = .fname)
    }) %>%
  mutate(year = substr(file,
                       start = stri_locate_first(file, regex = 'dem_')[1] + 4,
                       stop = nchar(file) - nchar('.csv'))) %>%
  # only keep relevant columns 
  select(year, Latitude, Longitude, Elevation, Tave01, Tave02, Tave03,
         Tave04, Tave05, Tave06, Tave07, Tave08, Tave09, Tave10, Tave11, Tave12,
         PPT01, PPT02, PPT03, PPT04, PPT05, PPT06, PPT07, PPT08, PPT09, PPT10,
         PPT11, PPT12) %>%
  # pivot from wide to long format (only one column of precip and temp)
  pivot_longer(-c(year, Latitude, Longitude, Elevation),
               names_to = 'parameter', values_to = 'value') %>%
  # extract month and year out of parameters
  mutate(month = map_chr(parameter,
                         \(.chr) substr(.chr, nchar(.chr) - 1, nchar(.chr))),
         dec_date = decimal_date(date(paste(year, month, '15', sep = '-'))),
         month = as.numeric(month),
         year = as.numeric(year),
         parameter = map_chr(parameter,
                             \(.chr) substr(.chr, 1, nchar(.chr) - 2))) %>%
  # pivot wider to make separate columns of temperature and precipitation
  pivot_wider(names_from = parameter, values_from = value) %>%
  # convert monthly total precipitation to average daily precipitation
  mutate(first_day = as.Date(paste(year, month, '01', sep = '-')),
         next_month = if_else(month != '12', as.numeric(month + 1), 1),
         next_year = if_else(month != '12', year, year + 1),
         last_day = as.Date(paste(next_year, next_month, '01', sep = '-')),
         samples = as.numeric((last_day - first_day)),
         avgprecip = PPT / samples) %>% # convert to millimeters per day
  # drop temporary columns
  select(-c(first_day, next_month, next_year, last_day, samples, PPT)) %>%
  # change to names used in the models
  rename(avgtemp = Tave,
         latitude = Latitude,
         longitude = Longitude,
         elevation = Elevation) %>%
  relocate(c(month, dec_date), .after = year)
saveRDS(historical_climate_data, file = "Data/HWI/ClimateNA_v731/climate data/historical_climate_data.rds")
```


```{r}
#adding park latitude and longitude coordinates to dataframe
data <- left_join(data, nationalparks_bc_coordinates, by = 'park')
data$year <- lubridate::year(data$year_month)
#extract month from date then create a column for month, required for climate data
data$month <- lubridate::month(data$year_month)
# add historical climate data to the dataframe
data <- left_join(data, historical_climate_data, by = c("latitude","longitude", "year", "month"))
#data <- relocate(data, dec_date, .after = month)
```


```{r}
#Plot historical number of interaction recordings over time
plot(data$HWI ~ data$year_month, xlab = "Time", ylab = "Human-wildlife interaction")

#Plot historical seasonal trends of interaction recordings over time
#Glacier: light blue; Kootenay: orange; Revelstoke: green; Pacific Rim: dark blue; Yoho: purple
colour_park <- c("#66CCEE", "#EE7733", "#228833", "#004488", "#AA4499")

ggplot() +
  geom_jitter(data = data, 
              aes(y = HWI, x = month, col = park),
              alpha = 0.05, size = 1, width = 0.25, shape = 17) + #point for every observation
  geom_smooth(data = data, 
              aes(y = HWI, x = month, col = park),
              linewidth = 0.6, se = F) + #line for each park
  scale_x_continuous(breaks = seq_along(month.abb), labels = month.abb) +
  scale_y_continuous(limits = c(0, 10)) +
  scale_color_manual(name = "National Parks", values = colour_park,
                     labels = c("Glacier",
                                "Kootenay",
                                "Revelstoke",
                                "Pacific Rim",
                                "Yoho")) +
  xlab("Month") +
  ylab("Human-wildlife interaction") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=10, family = "sans", face = "bold"),
        axis.title.x = element_text(size=10, family = "sans", face = "bold"),
        axis.text.y = element_text(size=7, family = "sans"),
        axis.text.x  = element_text(size=7, family = "sans"),
        axis.ticks.x = element_blank(),
        plot.title = element_text(size = 25, family = "sans", face = "bold",
                                  vjust = -5.5, hjust = 0.02),
        #legend.position = "none",
        legend.box.background = element_rect(color = "black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0,0.2,0,0.2), "cm"))
```

# Modelling

```{r include = FALSE}
data$park <- as.factor(data$park)
#data$species <- as.factor(data$species)
data$year_month <- as.Date(data$year_month, format = "%Y-%m-%d")
```


```{r eval = FALSE}
data$dec_year_month <- decimal_date(data$year_month)
#Build the model
model <- gam(HWI ~
               #s(dec_year_month, k = 3) +                #global effect of time
               s(avgtemp) +                               
               s(log(avgprecip + 1e-10)) +                
               avgtemp:log(avgprecip + 1e-10) +
               month:avgtemp +
               month:log(avgprecip + 1e-10) +
               s(month, park, bs = 'fs', xt = list(bs = 'cc'), k = 5),     
             family = poisson(link = "log"),              
             data = data,
             method = "REML",
             control = gam.control(nthreads = 8, trace = TRUE),
             knots = list(month = c(0.5, 12.5)))
summary(model)
```

# Model validation

```{r}
#Create new dataframe for each park
group <- data %>% 
  group_by(park) %>%
  groups %>%
  unlist %>% 
  as.character

#Glacier National Park of Canada
park_glacier <- data[which(data$park == "Glacier National Park of Canada"),] %>% 
  group_by(park) %>% 
  summarise() %>% 
  slice_sample(n = 2) %>% 
  mutate(unique_id=1:NROW(.))
#create new dataframe
newdata_glacier <- data[which(data$park == "Glacier National Park of Canada"),] %>% 
  group_by(park)  %>% 
  right_join(park_glacier, by=group) %>%
  group_by_(group)

#Kootenay National Park of Canada
park_kootenay <- data[which(data$park == "Kootenay National Park of Canada"),] %>% 
  group_by(park) %>% 
  summarise() %>% 
  slice_sample(n = 2) %>% 
  mutate(unique_id=1:NROW(.))
newdata_kootenay <- data[which(data$park == "Kootenay National Park of Canada"),] %>% 
  group_by(park)  %>% 
  right_join(park_kootenay, by=group) %>%
  group_by_(group)

#Mount Revelstoke National Park of Canada
park_revelstoke <- data[which(data$park == "Mount Revelstoke National Park of Canada"),] %>% 
  group_by(park) %>% 
  summarise() %>% 
  slice_sample(n = 2) %>% 
  mutate(unique_id=1:NROW(.))
newdata_revelstoke <- data[which(data$park == "Mount Revelstoke National Park of Canada"),] %>% 
  group_by(park)  %>% 
  right_join(park_revelstoke, by=group) %>%
  group_by_(group)

#Pacific Rim National Park Reserve of Canada
park_pacific_rim <- data[which(data$park == "Pacific Rim National Park Reserve of Canada"),] %>% 
  group_by(park) %>% 
  summarise() %>% 
  slice_sample(n = 2) %>% 
  mutate(unique_id=1:NROW(.))
newdata_pacific_rim <- data[which(data$park == "Pacific Rim National Park Reserve of Canada"),] %>% 
  group_by(park)  %>% 
  right_join(park_pacific_rim, by=group) %>%
  group_by_(group)

#Yoho National Park of Canada
park_yoho <- data[which(data$park == "Yoho National Park of Canada"),] %>% 
  group_by(park) %>% 
  summarise() %>% 
  slice_sample(n = 2) %>% 
  mutate(unique_id=1:NROW(.))
newdata_yoho <- data[which(data$park == "Yoho National Park of Canada"),] %>% 
  group_by(park)  %>% 
  right_join(park_yoho, by=group) %>%
  group_by_(group)

#Generate predictions for each park
newdata_glacier$predict <- predict(model, newdata = newdata_glacier, type = "response")
newdata_kootenay$predict <- predict(model, newdata = newdata_kootenay, type = "response")
newdata_revelstoke$predict <- predict(model, newdata = newdata_revelstoke, type = "response")
newdata_pacific_rim$predict <- predict(model, newdata = newdata_pacific_rim, type = "response")
newdata_yoho$predict <- predict(model, newdata = newdata_yoho, type = "response")

```


```{r}
#Plot to visually compare historical recordings vs. model-predicted interactions to make sure model is behaving
#Glacier
ggplot() +
  geom_boxplot(data = newdata_glacier, aes(x = month, y = HWI, group = cut_width(month, 1))) +
  geom_point(data = newdata_glacier, aes(x = month, y = predict), size = 2, alpha = 0.5, col = "#66CCEE") +
  scale_color_viridis(discrete = T) +
  scale_x_continuous(breaks = seq_along(month.abb), labels = month.abb) +
  xlab("Month") +
  ylab("Number of Interaction (Monthly)") +
  ggtitle("Glacier National Park of Canada: Historical vs. model-predicted interaction recordings") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=9, family = "sans"),
        axis.text.x  = element_text(size=9, angle = 45, family = "sans"),
        plot.title = element_text(hjust = -0.05, size = 12, family = "sans", face = "bold"),
        legend.position = "right",
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

#Kootenay
ggplot() +
  geom_boxplot(data = newdata_kootenay, aes(x = month, y = HWI, group = cut_width(month, 1))) +
  geom_point(data = newdata_kootenay, aes(x = month, y = predict), size = 2, alpha = 0.5, col = "#EE7733") +
  scale_color_viridis(discrete = T) +
  scale_x_continuous(breaks = seq_along(month.abb), labels = month.abb) +
  xlab("Month") +
  ylab("Number of Interaction (Monthly)") +
  ggtitle("Kootenay National Park of Canada: Historical vs. model-predicted interaction recordings") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=9, family = "sans"),
        axis.text.x  = element_text(size=9, angle = 45, family = "sans"),
        plot.title = element_text(hjust = -0.05, size = 12, family = "sans", face = "bold"),
        legend.position = "right",
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

#Revelstoke
ggplot() +
  geom_boxplot(data = newdata_revelstoke, aes(x = month, y = HWI, group = cut_width(month, 1))) +
  geom_point(data = newdata_revelstoke, aes(x = month, y = predict), size = 2, alpha = 0.5, col = "#228833") +
  scale_color_viridis(discrete = T) +
  scale_x_continuous(breaks = seq_along(month.abb), labels = month.abb) +
  xlab("Month") +
  ylab("Number of Interaction (Monthly)") +
  ggtitle("Mount Revelstoke National Park of Canada: Historical vs. model-predicted interaction recordings") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=9, family = "sans"),
        axis.text.x  = element_text(size=9, angle = 45, family = "sans"),
        plot.title = element_text(hjust = -0.05, size = 12, family = "sans", face = "bold"),
        legend.position = "right",
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

#Pacific Rim
ggplot() +
  geom_boxplot(data = newdata_pacific_rim, aes(x = month, y = HWI, group = cut_width(month, 1))) +
  geom_point(data = newdata_pacific_rim, aes(x = month, y = predict), size = 2, alpha = 0.5, col = "#004488") +
  scale_color_viridis(discrete = T) +
  scale_x_continuous(breaks = seq_along(month.abb), labels = month.abb) +
  xlab("Month") +
  ylab("Number of Interaction (Monthly)") +
  ggtitle("Mount Revelstoke National Park of Canada: Historical vs. model-predicted interaction recordings") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=9, family = "sans"),
        axis.text.x  = element_text(size=9, angle = 45, family = "sans"),
        plot.title = element_text(hjust = -0.05, size = 12, family = "sans", face = "bold"),
        legend.position = "right",
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

#Yoho
ggplot() +
  geom_boxplot(data = newdata_yoho, aes(x = month, y = HWI, group = cut_width(month, 1))) +
  geom_point(data = newdata_yoho, aes(x = month, y = predict), size = 2, alpha = 0.5, col = "#AA4499") +
  scale_color_viridis(discrete = T) +
  scale_x_continuous(breaks = seq_along(month.abb), labels = month.abb) +
  xlab("Month") +
  ylab("Number of Interaction (Monthly)") +
  ggtitle("Yoho National Park of Canada: Historical vs. model-predicted interaction recordings") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=9, family = "sans"),
        axis.text.x  = element_text(size=9, angle = 45, family = "sans"),
        plot.title = element_text(hjust = -0.05, size = 12, family = "sans", face = "bold"),
        legend.position = "right",
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

```

```{r}
# Plot monthly interactions in response to temperature
temp_HWI <- 
  ggplot(data = data, 
         aes(y = HWI, x = avgtemp)) +
    geom_point(aes(col = park), alpha = 0.1, pch = 16) +
  geom_smooth(method = "gam",
              formula = y ~ s(x, bs = "cs",
                              k = 4),
              method.args = list(family = poisson),
              se = FALSE,
              col = "black",
              size = 1.2) + # line for overall trend
  geom_smooth(aes(col = park),
              method = "gam",
              formula = y ~ s(x, bs = "cs",
                              k = 4),
              method.args = list(family = poisson),
              se = FALSE,
              linewidth = 0.7) + 
  scale_x_continuous(breaks = seq(-15,15,5)) + 
  scale_y_continuous(limits = c(0,30), expand = c(0,0.2)) +
  scale_color_manual(name = "National Parks", values = colour_park,
                     labels = c("Glacier",
                                "Kootenay",
                                "Revelstoke",
                                "Pacific Rim",
                                "Yoho")) + 
  xlab("Average Monthly Temperature (ºC)") +
  ylab("Total Monthly HWIs") +
  ggtitle("A") +
  guides(col = guide_legend(override.aes = list(alpha=1))) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=8, family = "sans"),
        axis.text.x  = element_text(size=8, family = "sans"),
        plot.title = element_text(hjust = 0.04, vjust = -7, size = 16, 
                                  family = "sans", face = "bold"),
        legend.position=c(0.1,0.5),
        legend.box.background = element_rect(color = "black"),
        legend.title = element_text(face = "bold"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

# Plot monthly interactions in response to precipitation
precip_HWI <- 
  ggplot(data = data, 
         aes(y = HWI, x = avgprecip)) +
    geom_point(aes(col = park), alpha = 0.1, pch = 16) +
    geom_smooth(method = "gam",
                formula = y ~ s(x, bs = "cs",
                                k = 4),
                method.args = list(family = poisson),
                se = FALSE,
                col = "black",
                size = 1.2) + # line for overall trend
    geom_smooth(aes(col = park),
                method = "gam",
                formula = y ~ s(x, bs = "cs",
                                k = 4),
                method.args = list(family = poisson),
                se = FALSE,
                linewidth = 0.7) + 
    scale_y_continuous(limits = c(0,30), expand = c(0,0.2)) +
    scale_color_manual(name = "National Parks", values = colour_park,
                       labels = c("Glacier",
                                  "Kootenay",
                                  "Revelstoke",
                                  "Pacific Rim",
                                  "Yoho")) + 
    xlab("Average Monthly Precipitation (mm)") +
    ylab("Total Monthly HWIs") +
    ggtitle("B") +
    guides(col = guide_legend(override.aes = list(alpha=1))) +
    theme_bw() +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.title.y = element_text(size=12, family = "sans", face = "bold"),
          axis.title.x = element_text(size=12, family = "sans", face = "bold"),
          axis.text.y = element_text(size=8, family = "sans"),
          axis.text.x  = element_text(size=8, family = "sans"),
          plot.title = element_text(hjust = 0.04, vjust = -7, size = 16, 
                                    family = "sans", face = "bold"),
          legend.position = "none",
          panel.background = element_rect(fill = "transparent"),
          plot.background = element_rect(fill = "transparent", color = NA),
          plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

TOP <- grid.arrange(temp_HWI, precip_HWI, 
                     ncol=2)
```

# Projecting

```{r}
# Import clean projected climate data
#refer to "projected climate data" section within the 'elevation and climate.R' script
projection_climate_data <- readRDS("Data/HWI/ClimateNA_v731/climate data/projection_climate_data_monthly.rds")
projection_climate_data <- merge(projection_climate_data, nationalparks_bc_coordinates,
                                 by = c("latitude","longitude"), all.x=TRUE)
names(projection_climate_data)[7] <- "avgtemp"
names(projection_climate_data)[6] <- "dec_year_month"
```


```{r}
#Project the number of interactions under 4 SSP climate scenarios
projection_climate_data$predicted_interactions <- predict(model, newdata = projection_climate_data, type = "response")

#Calculate sum of interactions in each year then scale relative to 2020
agg_proj <- aggregate(predicted_interactions ~ year + park + scenario , data = projection_climate_data, FUN = "sum")

data2 <- agg_proj[which(agg_proj$year == 2022),]
data2 <- aggregate(predicted_interactions ~ park, data = data2, FUN = "median")

#Kootenay
agg_proj[which(agg_proj$park == "Kootenay National Park of Canada"),"predicted_interactions"] <- 
  agg_proj[which(agg_proj$park == "Kootenay National Park of Canada"),"predicted_interactions"]/data2[which(data2$park == "Kootenay National Park of Canada"),"predicted_interactions"]

#Pacific Rim
agg_proj[which(agg_proj$park == "Pacific Rim National Park Reserve of Canada"),"predicted_interactions"] <- 
  agg_proj[which(agg_proj$park == "Pacific Rim National Park Reserve of Canada"),"predicted_interactions"]/data2[which(data2$park == "Pacific Rim National Park Reserve of Canada"),"predicted_interactions"]

#Glacier
agg_proj[which(agg_proj$park == "Glacier National Park of Canada"),"predicted_interactions"] <- 
  agg_proj[which(agg_proj$park == "Glacier National Park of Canada"),"predicted_interactions"]/data2[which(data2$park == "Glacier National Park of Canada"),"predicted_interactions"]

#Yoho
agg_proj[which(agg_proj$park == "Yoho National Park of Canada"),"predicted_interactions"] <- 
  agg_proj[which(agg_proj$park == "Yoho National Park of Canada"),"predicted_interactions"]/data2[which(data2$park == "Yoho National Park of Canada"),"predicted_interactions"]

#Revelstoke
agg_proj[which(agg_proj$park == "Mount Revelstoke National Park of Canada"),"predicted_interactions"] <- 
  agg_proj[which(agg_proj$park == "Mount Revelstoke National Park of Canada"),"predicted_interactions"]/data2[which(data2$park == "Mount Revelstoke National Park of Canada"),"predicted_interactions"]
```

# Visualizing Results

```{r}
#Plot results
#Make colour strips in x-direction for panel title boxes (they will correspond to SSPs)
strip <- strip_themed(background_x = 
                        elem_list_rect(fill = c("#0571b0", "#92c5de", "#f4a582", "#ca0020")))

#Assign a letter for each SSP so we can label each panel
data_text <- data.frame(label = c("A", "B", "C", "D"),
                        scenario = names(table(projection_climate_data$scenario)),
                        x = c(2023, 2023, 2023, 2023), 
                        y = c(1.9, 1.9, 1.9, 1.9))

BOT <- 
ggplot(agg_proj, aes(x = year , y = predicted_interactions)) + 
  geom_hline(yintercept = 1, size = 0.5, color = "grey70") + # line at 1 for reference
  geom_line(linewidth=0.5, aes(group = park, col = park), alpha = 0.8) + # line for each park
  facet_wrap2(~ scenario, strip = strip, labeller = as_labeller(c(
    `ssp126_low` = "SSP 1-2.6",
    `ssp245_intermediate` = "SSP 2-4.5",
    `ssp370_high` = "SSP 3-7.0",
    `ssp585_veryhigh` = "SSP 5-8.5")),
    nrow = 1) + # create a panel for each climate change scenario
  xlab("Year") +
  ylab("Relative change in annual HWIs") +
  #scale_y_continuous(limits = c(0.75,2)) +
  guides(col = guide_legend(override.aes = list(alpha=1))) +
  scale_color_manual(name = "National Parks", values = colour_park,
                     labels = c("Glacier",
                                "Kootenay",
                                "Revelstoke",
                                "Pacific Rim",
                                "Yoho")) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text = element_text(face = "bold"),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=8, family = "sans"),
        axis.text.x  = element_text(size=8, family = "sans"),
        plot.title = element_text(size = 12, family = "sans", face = "bold"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.position = "none",
        legend.key.size = unit(0.4, 'cm'),
        legend.title = element_text(size = 8, face = "bold"),
        legend.text=element_text(size=7),
        legend.box.background = element_rect(color = "black"),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

FIG <- grid.arrange(TOP, BOT, 
                    ncol=1, heights = c(3,2.7))
```

```{r }
ggsave(temp, filename = "Figures/Supplementary Figures/HWI/HWI_temp.png", device = NULL,
       path = NULL, scale = 1, width = 6.86, height = 6, units = "in", dpi = 600)

ggsave(precip, filename = "Figures/Supplementary Figures/HWI/HWI_precip.png", device = NULL,
       path = NULL, scale = 1, width = 6.86, height = 6, units = "in", dpi = 600)

ggsave(TOP, filename = "Figures/Supplementary Figures/HWI/HWI_temp_precip.png", device = NULL,
       path = NULL, scale = 1, width = 14.91, height = 6.47, units = "in", dpi = 600)

ggsave(BOT, filename = "Figures/Supplementary Figures/HWI/projections.png", device = NULL,
       path = NULL, scale = 1, width = 14.91, height = 6.47, units = "in", dpi = 600)

ggsave(FIG, filename = "Figures/figure5.png", device = NULL,
       path = NULL, scale = 1, width = 14.91, height = 10, units = "in", dpi = 600)
```