---
title: "Sea level rise"
author: "Willem Stolte"
date: "5/5/2022"
output: html_document
---

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

datadir <- "p:/11202493--systeemrap-grevelingen/1_data/Noordzee"
hoofdstations <- c(
  "DELFZL",
  "DENHDR",
  "HARLGN",
  "HOEKVHLD",
  "IJMDBTHVN",
  "VLISSGN"
)

filename = paste0("maandwaterhoogte2022-03-11.csv")
# monthlyStat <- read_delim(file.path(datadir, "ddl", "standard", filename), delim = ";") %>%
  # filter(monthlyaverage < 100)
# save locally for offline use
# write_delim(monthlyStat, filename, delim = ";")
monthlyStat <- read_delim(file.path(filename), delim = ";") 

```

```{r}
p <- monthlyStat %>% 
  filter(locatie.code %in% hoofdstations) %>%
  mutate(date = lubridate::ymd(paste(year, month, "01"))) %>%
  ggplot() +
  geom_point(aes(x = date, y = n, color = locatie.code))
ggplotly(p)
```




```{r}

filteryears <- c(1970:1993)

monthlyStat %>%
  filter(locatie.code %in% hoofdstations) %>%
  filter(n > 3000) %>%
  filter(year %in% filteryears) %>%
  ggplot(aes(year, monthlyaverage)) +
  # geom_point(aes(color = month)) +
  geom_smooth(aes(color = as.factor(month)), method = "lm", alpha = 0.1) +
  facet_wrap(~ locatie.code) 

monthlyStat %>%
  filter(locatie.code %in% hoofdstations) %>%
  filter(n > 3000) %>%
  filter(year %in% filteryears) %>%
  ggplot(aes(month, monthlyaverage)) +
  geom_boxplot(aes(group = month)) +
  facet_wrap(~ locatie.code) +
  scale_x_continuous(breaks = scales::pretty_breaks()) +
  ylab("monthly average sea level")

fitted_models = monthlyStat %>% 
  filter(n > 3000) %>%
  filter(year %in% filteryears) %>%
  group_by(locatie.code, month) %>% 
  do(model = lm(monthlyaverage ~ year, data = .))

metadata <- fitted_models %>%
  select(locatie.code, month) %>%
  as_tibble()

coefficients <- map_df(fitted_models$model, broom::tidy) %>%
  filter(term == "year") %>%
  bind_cols(metadata)

coefficients %>%
  filter(locatie.code %in% hoofdstations) %>%
  mutate(average = mean(estimate)) %>%
  group_by(locatie.code) %>%
  mutate(stationaverage = mean(estimate)) %>% ungroup() %>%
  group_by(month) %>%
  mutate(monthlyaverage = mean(estimate)) %>% ungroup() %>%
  ggplot(aes(month, estimate)) +
  geom_point() +
  geom_smooth(method = "lm") +
  # geom_smooth(method = "lm") +
  geom_line(aes(month, average)) +
  annotate(geom = "text", x = 6, y = 0.17, label = "average sea level change") +
  facet_wrap(~ locatie.code) +
  scale_x_continuous(breaks = scales::pretty_breaks()) +
  ylab("change of sealevel per year for each months (cm/year)")

monthlyStat %>%
  filter(locatie.code %in% hoofdstations) %>%
  group_by(month, locatie.code) %>%
  summarize(`mean height per month` = mean(monthlyaverage)) %>% 
  left_join(coefficients) %>%
  ggplot(aes(estimate, `mean height per month`)) + 
  # geom_point(aes(color = as.factor(month)), size = 3) +
  geom_label(aes(label = month), size = 3) +
  geom_smooth(method = "lm", alpha = 0.2) +
  facet_wrap(~ locatie.code) +
  xlab("gemiddelde zeespiegelstijging, voor elke maand")


```

