03_volume_calcs.Rmd

---
title: "03_volume_calcs"
output: html_document
editor_options: 
  chunk_output_type: console
---

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

# Pull in the munged lake height and lake area data

```{r}
areas <- read_csv('data/out/areas_munged_CC30.csv') %>%
  select(area = Reconstructed, sat, name, date = timestamp) %>%
  mutate(date = as_date(date),
         areaID = row_number()) %>%
  group_by(name, date) %>%
  summarise(area = mean(area)) %>%
  ungroup()

heights <- read_csv('data/out/heights_munged.csv') %>%
  select(height, name = name.std, date)

lake.summaries <- read_csv('data/in/lake_properties.csv') %>%
  distinct(name,.keep_all = T)

```

## Find matchups (+/- 1 day pairs)

```{r}
sameday <- areas %>% inner_join(heights) %>%
  group_by(name, date) %>%
  summarise(area = mean(area),
         height = mean(height)) %>%
  mutate(obID = paste0(date,"_",name)) %>%
  ungroup()

plusone <- areas %>%
  mutate(area_date = date,
         date = date + 1) %>%
  inner_join(heights) %>%
  group_by(name, date) %>%
  summarise(area = mean(area),
         height = mean(height)) %>%
  mutate(obID = paste0(date,"_",name)) %>%
  ungroup()

minusone <- areas %>%
  mutate(area_date = date,
         date = date - 1) %>%
  inner_join(heights) %>%
  group_by(name, date) %>%
  summarise(area = mean(area),
         height = mean(height)) %>%
  mutate(obID = paste0(date,"_",name)) %>%
  ungroup()

## We prefer same day,
plusone <- plusone %>% filter(!obID %in% sameday$obID)
minusone <- minusone %>% filter(!obID %in% sameday$obID)

## Check for obs both in both plus and minus and average them
overlap <- plusone %>% filter(obID %in% minusone$obID) %>%
  group_by(name, date, obID) %>%
  summarise(height = mean(height),
            area = mean(area))

matchups <- sameday %>% mutate(matchup = 'same') %>%
  bind_rows(plusone %>% filter(!obID %in% overlap$obID) %>% mutate(matchup = 'plusone')) %>%
  bind_rows(minusone %>% filter(!obID %in% overlap$obID) %>% mutate(matchup = 'minusone')) %>%
  bind_rows(overlap %>% mutate(matchup = 'average')) %>%
  mutate(area_m = area*1000000)


## This leaves us with only 73 lakes where there is at least 1 matchup
dplyr::count(matchups,name) %>% filter(n>2)
## and 68 lakes with at least 3
```

## Calculate volume changes

```{r}
## Simple volumetric change based on trapezoidal volume
volChange <- function(lake){

  date.first <- lake[lake$date == min(lake$date),]
  
  lake %>% arrange(date) %>%
    mutate(area_diff = area_m + date.first$area_m[1],
           height_diff = height - date.first$height[1],
           vol_change = height_diff*(area_diff/2))
  }
  
volumes <- matchups %>% group_by(name) %>%
  nest() %>%
  mutate(volCalcs = map(data, volChange)) %>%
  select(-data) %>%
  unnest(cols = c(volCalcs))


## Rating curve calculations to expand volume change to all the height measurements
filter <- dplyr::count(volumes, name) %>% filter(n<2)

fits <- volumes %>%
  filter(!name %in% filter$name) %>%
  group_by(name) %>%
  nest() %>%
  mutate(fit = map(data, function(df) lm(vol_change ~ height, data = df))) %>%
  select(-data)

volumes_fitted <- heights %>%
  nest(-name) %>%
  inner_join(fits) %>%
  mutate(vol_fit = map2(fit, data, ~predict(.x, newdata = .y) %>% unname())) %>%
  select(-fit) %>%
  unnest(cols = c(data, vol_fit))


write_csv(volumes_fitted, 'data/out/volumes_fitted_CC30.csv')
```

## Join the the areas and look at some timeseries

```{r}
areas <- read_csv('data/out/areas_munged_CC30.csv') %>%
  select(area = Reconstructed, sat, name, date = timestamp) %>%
  mutate(date = as_date(date))

volsLong <- volumes_fitted %>%
  full_join(areas) %>%
  #filter(area > 5) %>%
  select(name, area, date, height, vol_fit) %>%
  pivot_longer(c(area, height, vol_fit), names_to = 'Metric') %>%
  mutate(Metric = factor(Metric, levels = c('height', 'area', 'vol_fit')))

lake <- 'Salters Lake'

volsLong %>% filter(name == lake,  date > '2017-04-01') %>%
  ggplot(aes(x = date, y = value)) + geom_path() + geom_point()   +
  facet_wrap(~Metric, scales = 'free_y', ncol = 1) +
  ggtitle(lake)


## Make them individually for correct lines

p1 <- ggplot(heights %>% filter(name == lake), aes(x = date, y = height)) + geom_line(color = 'blue') +
  labs(y = 'Level (m)', subtitle = 'a) Lake Level') +
  theme_classic() +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_text(angle = 45, hjust = .5),
        plot.subtitle = element_text(hjust = 0.5))
  
p2 <- ggplot(areas %>% filter(name == lake, date > '2017-04-18'), aes(x = date, y = area)) + geom_line(color = 'red') +
  labs(y = expression(Surface~Area~(km^2)), subtitle = 'b) Lake Surface Area') +
  theme_classic() +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_text(angle = 45, hjust = .5),
        plot.subtitle = element_text(hjust = 0.5))

p3 <- ggplot(volumes_fitted %>% filter(name == lake), aes(x = date, y = vol_fit)) + geom_line(color = 'black') +
  labs(y = expression(Delta~Volume~(m^2)), subtitle = 'c) Change in Lake Volume') +
  theme_classic() +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_text(angle = 45, hjust = .5),
        plot.subtitle = element_text(hjust = 0.5))

g <- gridExtra::grid.arrange(p1, p2,p3, ncol = 1)

ggsave('figures/Salters3panel.png', plot = g, width = 5, height = 4, units = 'in', dpi = 600)

```

## Rating curve

```{r}
ggplot(volumes %>% filter(name == "Bay Tree Lake"), aes(x = height_diff, y = vol_change)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  ggpmisc::stat_poly_eq(formula = y~x, rr.digits = 3,coef.digits = 5, label.x = 'right', label.y = 'bottom',
               aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")), 
               parse = TRUE, color = 'red') +
  labs(y = expression(Delta~Volume~(m^2)), x = expression(Delta~Height~(m))) +
  theme_bw() +
  theme(axis.text.y = element_text(angle = 45, hjust = .5))

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

```

## Make the summary table

```{r}

satCounts <- areas %>%
  filter(timestamp > '2015-01-01',timestamp < '2020-01-01') %>%
  group_by(name) %>%
  summarise(countAreas = n(),
            mean = mean(Reconstructed))

heightsCounts <- heights %>% group_by(name) %>%
  summarise(countLevels = n(),
            dateFirst = min(date),
            dateLast = max(date))

tabOut <- lake.summaries %>% select(name, region) %>%
  left_join(heightsCounts) %>%
  left_join(satCounts)

## West Loon and East Loon are the same
tabOut$countLevels[tabOut$name == 'West Loon Lake'] <- tabOut$countLevels[tabOut$name == 'East Loon Lake']
tabOut$dateFirst[tabOut$name == 'West Loon Lake'] <- tabOut$dateFirst[tabOut$name == 'East Loon Lake']
tabOut$dateLast[tabOut$name == 'West Loon Lake'] <- tabOut$dateLast[tabOut$name == 'East Loon Lake']

write_csv(tabOut, 'data/out/SummaryTable.csv')

```