4_trendAnalysis.Rmd

---
title: "4_trendAnalysis"
output: html_document
date: "2023-03-13"
editor_options: 
  chunk_output_type: console
---

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

# Libraries
```{r}
library(tidyverse)
library(sf)
library(lubridate)
library(grDevices)
library(mapview)
library(extrafont)
library(ggpubr)
library(ggmap)
library(RgoogleMaps)
library(broom)
library(feather)
library(tidyhydat)
library(sp)
library(data.table)
library(ggalluvial)
library(patchwork)
library(magick)
library(units)
library(Kendall)
library(ggspatial)
library(dtplyr)
#Import libraries for Random Forest
library(caret) 
library(e1071)
library(Boruta)
library(tidymodels)
library(skimr)
library(vip)
```

# Import files / set constants
```{r}
# dates for version control
todayDate  = "20230324" # the first data join phase

# Names of files and folders for reflectance data
import.filePath = "C:/Users/whyana/OneDrive/DocumentsLaptop/001_GraduateSchool/Research/Connectivity/Mackenzie/Data/GEE Downloads"


# intermediate working directory
int.wd="C:/Users/whyana/OneDrive/DocumentsLaptop/001_GraduateSchool/Research/Connectivity/Mackenzie/Data/intermediaryDownloads"
refl.import = 'srCorrected_mackLakes_202303138.feather'

#Name of file and folder for lake shapefiles & island polygon shapefiles
shapeFiles.filePath = "C:/Users/whyana/OneDrive/DocumentsLaptop/001_GraduateSchool/Research/Connectivity/Mackenzie/Data/shapeFiles"
lakes.shapeFile = "mackenzieGoodLakes.shp"
setwd(shapeFiles.filePath)
lakes.sf = st_read(lakes.shapeFile)
import.sword = "na_sword_reaches_hb82_v14.shp"

images.wd = "C:/Users/whyana/OneDrive/DocumentsLaptop/001_GraduateSchool/Research/Connectivity/Mackenzie/images"
```

# Import river centerlines and set the projection for all future plots, import classifications
```{r}
crs.plot = "+proj=tcea +lon_0=-134.3847656 +datum=WGS84 +units=m +no_defs"
setwd(shapeFiles.filePath)

study.area.large=cbind.data.frame(lon=c(-136.80, -136.80, -133.47, -133.47), 
                 lat=c(67.25, 69.55, 69.55, 67.46)) %>% 
  st_as_sf(coords=c("lon", "lat")) %>% st_set_crs(4326) %>% 
  st_bbox() %>% st_as_sfc() %>% 
  st_transform(crs = crs.plot)

mack.basin.large = st_read(import.sword) %>% 
  st_transform(crs = crs.plot) %>% 
  st_intersection(study.area.large) %>% dplyr::filter(width>90)

# import classifications
setwd(int.wd)
all.classified.filter = read_feather(paste0("final.class_", todayDate, ".feather")) 
```

# Trend analysis for calibrated reflected
## Prep lake connectivity classifications for trend analysis
```{r}
# group lakes within two temporal groups and summarize mean annual connectivity in each month
results.summary.subgroups = all.classified.filter %>% lazy_dt() %>% 
  dplyr::select(.pred_class, OBJECTID, year, month) %>% 
  mutate(yeargroup = case_when(
    year>=1984 & year<=2002 ~ "1984-2002",
    year>=2003 & year<=2022 ~ "2003-2022"
  )) %>% filter(!is.na(yeargroup)) %>% 
  group_by(OBJECTID, month, year, yeargroup)%>% 
  summarise(class.mean = mean(as.numeric(as.character(.pred_class)), na.rm=T),
            count=n()) %>% ungroup()
# calculate average yearly connectivity for all years
results.summary.all = all.classified.filter %>% lazy_dt() %>% 
  dplyr::select(.pred_class, OBJECTID, year, month) %>% 
  mutate(yeargroup = "all") %>% 
  group_by(OBJECTID, month, year, yeargroup)%>% 
  summarise(class.mean = mean(as.numeric(as.character(.pred_class)), na.rm=T),
            count=n()) %>% ungroup()

# Combine the two dataframes together
results.summary = rbind.data.frame(results.summary.subgroups %>% as_tibble(), 
                                   results.summary.all %>% as_tibble())

# group by time period, count number of years of data each lake has in each month in each period
good.ids = results.summary %>% group_by(OBJECTID, month, yeargroup) %>%count() %>% ungroup() %>% 
  filter(n>=10) 
# select only lakes that have at least 10 obs in all periods
best.ids = good.ids %>% group_by(OBJECTID, month) %>% count() %>% ungroup() %>% filter(n==3)

# Apply the best.ids filter, and group observations by lake, month, and yeargroup
nested.data = results.summary %>%
  left_join(best.ids, by=c("OBJECTID", "month")) %>% 
  dplyr::filter(!is.na(n)) %>% 
  group_by(OBJECTID, month, yeargroup) %>% nest() %>% ungroup() %>% as_tibble()

```

## Apply trend analysis calculations to each lake
```{r}
## for each lake, calculate the trend (tau) and pvalue
row.combo=NULL
for (i in 1:nrow(nested.data)){
  dat = nested.data$data[[i]] %>% arrange(year)
  OBJECTID = nested.data$OBJECTID[[i]]
  month = nested.data$month[[i]]
  yeargroup = nested.data$yeargroup[[i]]
  n.obs = nrow(dat)
  obs.count = dat %>% group_by(class.mean) %>% count() %>% ungroup() %>% 
    mutate(all.obs = n.obs,
           pct = n/n.obs)
  if(isTRUE(obs.count$pct[obs.count$class.mean<=0.66]>=0.95)){
    class = "always less than 0.66"
    col.combo = cbind.data.frame(OBJECTID, month,yeargroup,class, pval=NA, S=NA, tau=NA)
    row.combo=rbind.data.frame(row.combo, col.combo)
  } else if(isTRUE(obs.count$pct[obs.count$class.mean>0.66 |obs.count$class.mean<=1.33]>=0.95)){
    class = "always 0.66-1.33"
    col.combo = cbind.data.frame(OBJECTID, month,yeargroup,class, pval=NA, S=NA, tau=NA)
    row.combo=rbind.data.frame(row.combo, col.combo)
  }else if(isTRUE(obs.count$pct[obs.count$class.mean>1.33]>=0.95)){
    class = "always >1.33"
    col.combo = cbind.data.frame(OBJECTID, month,yeargroup,class, pval=NA, S=NA, tau=NA)
    row.combo=rbind.data.frame(row.combo, col.combo)
  } else {
    class = "trendtest"
    test.obj=MannKendall(dat$class.mean)
    S=test.obj$S[[1]]
    tau = test.obj$tau
    pval = test.obj$sl
    col.combo = cbind.data.frame(OBJECTID, month,yeargroup, class, pval, S, tau)
    row.combo=rbind.data.frame(row.combo, col.combo)
  }
}

setwd(int.wd)
write_feather(row.combo, "raw_trend.feather")

row.combo = read_feather("raw_trend.feather")

## Format trend data results
trend.data=row.combo %>% as_tibble()%>% 
  mutate(trend = case_when(
    tau>0 & pval < 0.05 ~ "increasing connectivity trend",
    tau<0 & pval < 0.05~ "decreasing connectivity trend",
    pval>0.05 ~ "no monotonic trend")) %>% 
  left_join(lakes.sf, by="OBJECTID") %>% st_as_sf() %>% 
   st_transform(crs = crs.plot)

write_feather(trend.data %>% as_tibble() %>% select(-geometry), 
              paste0(todayDate, "Treds_cal.feather"))

```

## Analyze/plot trends in connectivity
### plots not included in manuscript
```{r}
# Print a summary of trend results
setwd(int.wd)
row.combo = read_feather("raw_trend.feather") %>% as.data.table()
trend.data = read_feather(paste0(todayDate, "Treds_cal.feather")) %>% 
  left_join(lakes.sf, by="OBJECTID") %>% st_as_sf() %>% 
  st_transform(crs.plot)

june.combo = row.combo[month==6,]

june.spread = dcast(june.combo, OBJECTID ~ yeargroup, value.var=c("class")) 

june.spread %>% na.omit()

row.combo %>% as_tibble()%>% 
  mutate(trend = case_when(
    tau>0 & pval < 0.05 ~ "increasing sig. connectivity trend",
    tau<0 & pval < 0.05~ "decreasing sig. connectivity trend",
    pval>0.05 ~ "no monotonic trend",
    is.na(tau) & class == "always less than 0.66" ~ "always less than 0.66",
    is.na(tau) & class == "always 0.66-1.33" ~ "always 0.66-1.33",
    is.na(tau) & class == "always >1.33" ~ "always >1.33")) %>% 
  filter(month==6) %>% 
  group_by(month, yeargroup, trend) %>% count() %>% 
  spread(yeargroup, n) 


row.combo %>% filter(!is.na(tau)&month==6) %>% 
  mutate(group = case_when(tau<0 ~ "lt 0", tau> 0 ~ "gt 0", tau==0 ~ "0")) %>% 
  group_by(yeargroup, group) %>% count()

```