04_lake_correlations.Rmd

---
title: "04_lake_correlations"
output: html_document
editor_options: 
  chunk_output_type: console
---

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

```

## Read in our volume data and calculate correlations between lake pairs

```{r}
lake.summaries <- read_csv('data/in/lake_properties.csv') %>%
  distinct(name, region)

volumes <- read_csv('data/out/volumes_fitted_CC30.csv') %>%
  left_join(lake.summaries)

corrs <- function(lakeOne, lakeTwo){
  
  lake1 <- volumes %>% filter(name == lakeOne) %>% 
    select(date, vol1 = vol_fit)
  
  lake2 <- volumes %>% filter(name == lakeTwo) %>% 
    select(date, vol2 = vol_fit)
  
  sameday <- lake1 %>% inner_join(lake2)
    
  plusone <- lake1 %>% 
    inner_join(lake2 %>%
                 mutate(date = date + 1)) %>%
    filter(!date %in% sameday$date)
  
  minusone <- lake1 %>% 
    inner_join(lake2 %>%
                 mutate(date = date - 1)) %>%
    filter(!date %in% sameday$date)

  
  ## Check for obs both in both plus and minus and average them
  overlap <- plusone %>% filter(date %in% minusone$date) %>%
    group_by(date) %>%
    summarise(vol1 = mean(vol1),
              vol2 = mean(vol2))
  
  matchups <- sameday %>%
    bind_rows(plusone %>% filter(!date %in% overlap$date)) %>%
    bind_rows(minusone %>% filter(!date %in% overlap$date)) %>%
    bind_rows(overlap) %>%
    arrange(date)
  
  if(nrow(matchups) < 3){
    tibble(lake1 = lakeOne, lake2 = lakeTwo, rho = NA, 
           p.value = NA, matchups = nrow(matchups))
  }else{
    s.corr <- cor.test(matchups$vol1, matchups$vol2, method = 'spearman')
    
    tibble(lake1 = lakeOne, lake2 = lakeTwo, rho = s.corr$estimate, 
         p.value = s.corr$p.value, matchups = nrow(matchups))
  }
}

lakes.filtered <- lake.summaries %>% filter(name %in% volumes$name)

lakeCombs <- combn(lakes.filtered$name, 2) %>%
  t() %>% as_tibble() %>% left_join(lakes.filtered %>% rename(V1 = name)) %>%
  left_join(lakes.filtered %>% rename(V2 = name, region2 = region)) %>%
  filter(region == region2)

lakeCorrs <- map2_dfr(lakeCombs$V1,lakeCombs$V2, corrs) 

lakeCorrs <- lakeCorrs %>% left_join(lakes.filtered %>% rename(lake1 = name)) %>%
  filter(lake1 !=  'Lake Washington', lake2 != 'Lake Washington') ## Turns out Lake Washington is hightly managed
## So remove it

lakeCorrs <- lakeCorrs %>%
  mutate(region = factor(region, levels = c('IL','NC','WA','WI')))

mean(lakeCorrs$rho[lakeCorrs$matchups > 9])

ggplot(lakeCorrs %>% filter(matchups > 9), aes(x = region, y = rho, fill = region)) + 
  geom_boxplot() +
  scale_fill_brewer() +
  labs(y = bquote(rho ~ 'of volumetric change in lake pairs'), x = 'Region', fill = 'Regions') +
  theme_classic()
  

ggsave('figures/RegionalCorrs.png', width = 5, height = 3.5, units = 'in', dpi = 600)

dplyr::count(lakeCorrs %>% na.omit(), region)
```


## Calculate lake distances

```{r}
library(sf)
lakes.sf <- st_read('data/in/Locss_Lakes.shp') %>%
  select(name = GNIS_NAME, area = AREASQKM) %>%
  st_transform(crs = 5070) %>%
  filter(name != 'Lake Washington')

## Fix the names
lakes.sf$name[lakes.sf$name == 'Lake Mattamuskeet E'] <- 'Lake Mattamuskeet East'
lakes.sf$name[lakes.sf$name == 'Lake Mattamuskeet W'] <- 'Lake Mattamuskeet West'
lakes.sf$name[lakes.sf$name == 'Beaver Lakes'] <- 'Beaver Lake'
lakes.sf$name[lakes.sf$name == 'Defiance Lake'] <- 'Lake Defiance'
lakes.sf$name[lakes.sf$name == 'Deep Quarry'] <- 'Deep Quarry Lake'
lakes.sf$name[lakes.sf$name == 'Huntley Lake'] <- 'Timber Lake'
lakes.sf$name[lakes.sf$name == 'Deep Lake' & lakes.sf$area == 0.135] <- 'Deep Lake Wisc.'
lakes.sf$name[lakes.sf$name == 'Deep Lake' & lakes.sf$area == 0.157] <- 'Deep Lake Wash.'
lakes.sf$name[lakes.sf$name == 'Phantom Lake' & lakes.sf$area == 0.261] <- 'Phantom Lake Wash.'
lakes.sf$name[lakes.sf$name == 'Phantom Lake' & lakes.sf$area == 0.184] <- 'Phantom Lake Wisc.'
lakes.sf$name[lakes.sf$name == 'Loon Lake' & lakes.sf$area == 0.672] <- 'West Loon Lake'

lakes.sf <- lakes.sf %>% left_join(lake.summaries)


LakeDistances <- function(lake1, lake2){
  dist <- st_distance(lakes.sf %>% filter(name == lake1), lakes.sf %>% filter(name == lake2))[1] %>%
  as.integer()
}

distances <- lakeCorrs %>%
  mutate(dist = map2_dbl(lake1,lake2, LakeDistances))

distances  <- distances %>%
  mutate(region = factor(region, levels = c('NC','IL','WA','WI'), labels = c('North Carolina', 'Illinois', 'Washington', 'Wisconsin')),
         dist = dist/1e3)

corrs <- distances %>% filter(matchups > 9) %>% group_by(region) %>%
  nest() %>%
  mutate(r = map_dbl(data, ~cor.test(.$rho, .$dist, method = 'spearman')$estimate),
         p = map_dbl(data, ~cor.test(.$rho, .$dist,method = 'spearman')$p.value)) %>%
  select(-data)


p1 <- distances %>% filter(region == 'North Carolina', matchups > 9) %>% na.omit() %>%
  ggplot(aes(x = dist, y = rho)) +
  geom_point(aes(size = matchups)) +
  geom_smooth(method = 'lm', se = F) +
  scale_size_area(breaks = c(25,50,100,150, 200), max_size = 2) +
  labs(title = ~underline('North Carolina'), subtitle = expression(~rho~'= -0.02 P = 0.90'), tag = 'b)') +
  theme_bw() +
  theme(plot.title = element_text(hjust = .5),
        plot.subtitle = element_text(hjust = .5),
        legend.position = 'none',
        axis.title = element_blank())

p2 <- distances %>% filter(region == 'Illinois', matchups > 9) %>% na.omit() %>%
  ggplot(aes(x = dist, y = rho)) +
  geom_point(aes(size = matchups)) +
  geom_smooth(method = 'lm', se = F) +
  scale_size_area(breaks = c(25,50,100,150, 200), max_size = 2) +
  labs(title = ~underline('Illinois'), subtitle = expression(~rho~'= -0.54 P = <0.01'), tag = 'c)') +
  theme_bw() +
  theme(plot.title = element_text(hjust = .5),
        plot.subtitle = element_text(hjust = .5),
        legend.position = 'none',
        axis.title = element_blank())


p3 <- distances %>% filter(region == 'Washington', matchups > 9) %>% na.omit() %>%
  ggplot(aes(x = dist, y = rho)) +
  geom_point(aes(size = matchups)) +
  geom_smooth(method = 'lm', se = F) +
  scale_size_area(breaks = c(25,50,100,150, 200), max_size = 2) +
  labs(title = ~underline('Washington'), subtitle = expression(~rho~'= -0.37 P = <0.01'), tag = 'd)') +
  theme_bw() +
  theme(plot.title = element_text(hjust = .5),
        plot.subtitle = element_text(hjust = .5),
        legend.position = 'none',
        axis.title = element_blank())


p4 <- distances %>% filter(region == 'Wisconsin', matchups > 9) %>% na.omit() %>%
  ggplot(aes(x = dist, y = rho)) +
  geom_point(aes(size = matchups)) +
  geom_smooth(method = 'lm', se = F) +
  scale_size_area(breaks = c(25,50,100,150, 200), max_size = 2) +
  labs(title = ~underline('Wisconsin'), subtitle = expression(~rho~'= -0.19 P = 0.16'), tag = 'e)') +
  theme_bw() +
  theme(plot.title = element_text(hjust = .5),
        plot.subtitle = element_text(hjust = .5),
        legend.position = 'none',
        axis.title = element_blank())

cor.test(distances$rho[distances$matchups > 9],distances$dist[distances$matchups > 9], method = 'spearman' )

p5 <- distances %>% filter(matchups > 9) %>% na.omit() %>%
  ggplot(aes(x = dist, y = rho)) +
  geom_point(aes(size = matchups)) +
  geom_smooth(method = 'lm', se = F) +
  scale_size_area(breaks = c(25,50,100,150, 200), max_size = 2) +
  labs(title = ~underline('All Regions'), subtitle = expression(~rho~'= -0.26 P = <0.01'), tag = 'a)',
       size = 'Number of\nMeasurements') +
  theme_bw() +
  theme(plot.title = element_text(hjust = .5),
        plot.subtitle = element_text(hjust = .5),
        legend.position = 'bottom',
        legend.direction = 'vertical',
        axis.title = element_blank())

layout <- rbind(c(1,2,3),
                c(1,4,5))

g <- gridExtra::grid.arrange(p5,p1,p2,p3,p4, layout_matrix = layout, left = 'Correlation Coefficient (ρ)', bottom = 'Distance Between Lakes (km)', widths = c(.4,.3,.3))

ggsave('figures/DistCorrs.png', width = 8, height = 6, units = 'in', plot = g, dpi = 600)

write_csv(distances, 'data/out/DistCorrs.csv')

```