Data_analysis.Rmd

---
title: "Data analysis"
output: html_document
date: "`r Sys.Date()`"
---

```{r setup, include=FALSE}
renv::restore()
library(tidyverse)
library(mgcv)
library(lubridate)
library(timetk)
library(RcppRoll)
library(dLagM)
library(dlnm)
library(zoo)
theme_set(theme_gray())
library(knitr)
library(kableExtra)
library(sf)
library(tmap)
library(gridExtra)
library(lubridate)
library(timetk)
library(RcppRoll)
library(dLagM)
library(dlnm)
library(zoo)
theme_set(theme_gray())
library(js)
library(sjPlot)
library(verification)
library(ROCR)
library(pROC)
# Read data
source("data_read.R")
```
# Geo Mapping

```{r Spatial analysis}
theme_map <- theme_bw()+
  theme(panel.grid.major = element_line(colour = "transparent"))

plot(st_geometry(hydro_bavaria))

ggplot()+
  geom_sf(data = hydro_bavaria)



# Plot with catchment points
ggplot()+
  geom_sf(data = hydro_bavaria, fill = "transparent") +
  geom_sf(data = pegel_prop_sf, aes(color = as.factor(unique(ID))), size = 3)+
  labs(color = "Catchment ID")+
  scale_color_viridis_d()+
  theme_map
ggsave("./graphics/Max/Catchments_hydro_bavaria_dots.png", width = 12, height = 9)

# Plot hydrobavaria and admin bavaria
ggplot()+
  geom_sf(data = hydro_bavaria, fill = "white")+
  geom_sf(data = admin_bavaria, fill = "transparent", color = "red")+
  geom_sf(data = pegel_prop_sf, aes(color = as.factor(unique(ID))), size = 3)+
  labs(color = "Catchment ID")+
  scale_color_viridis_d()+
  theme_map
ggsave("./graphics/Max/hydro_admin_bavaria.png", width = 12, height = 9)

# Catchments with rivers and Donau/Main for orientation
waterways_three %>% 
  ggplot()+
  geom_sf(data = hydro_bavaria, fill = "transparent")+
  geom_sf(data = waterways_three, color = "blue", size = 1.5)+
  geom_sf(data = pegel_prop_sf, aes(color = as.factor(unique(ID))), size = 3)+
  labs(color = "Catchment ID")+
  scale_color_viridis_d()+
  theme_map

ggsave("./graphics/Max/map_rivers_catchments.png", width = 12, height = 9)

```

### tmap geographical overview
```{r}
tmap_mode("view")
tm_shape(hydro_bavaria)+
  tm_polygons(id = "NameString")+
  tm_shape(waterways_three)+
  tm_lines(col = "blue")+
  tm_shape(pegel_prop_sf)+
  tm_markers(size = 0.3)
```

# Low flow analysis ----
```{r Analysis of low flow events}
# First visualizations
plot_drainage <- function(dataset, year = c(2000), member_input = "kbe") {
  # Abflussmenge für unterschiedliche Pegel über die Zeit
  dataset %>% filter(YY %in% year, member == member_input) %>%
  ggplot(mapping = aes(x = date, y = drainage, color = waterlevel)) +
    geom_line() +
    scale_y_continuous(trans='log10') +
    geom_point(data = dataset %>% filter(YY %in% year, member == member_input, lowlevel == TRUE), 
               aes(x = date, y = drainage), color = "purple") +
    labs(title = "Abflussmenge für unterschiedliche Pegel über die Zeit",
       y = "Abflussmenge (in m³/s, Log-Skala)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2")
}

plot_drainage_unmarked <- function(dataset, year = c(2000), member_input = "kbe") {
  # Abflussmenge für unterschiedliche Pegel über die Zeit
  dataset %>% filter(YY %in% year, member == member_input) %>%
  ggplot(mapping = aes(x = date, y = drainage, color = waterlevel)) +
    geom_line() +
    scale_y_continuous(trans='log10') +
    labs(title = "Abflussmenge für unterschiedliche Pegel über die Zeit",
       y = "Abflussmenge (in m³/s, Log-Skala)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2")
}

plot_single_drainage <- function(dataset, year = c(2000), member_input = "kbe") {
  # Abflussmenge für einen Pegel über die Zeit
  dataset %>% filter(YY %in% year, member == member_input) %>%
  ggplot(mapping = aes(x = date, y = drainage)) +
    geom_line() +
    scale_y_continuous(trans='log10') +
    geom_line(mapping = aes(y = nm7q), color = "red") +
    geom_point(data = dataset %>% filter(YY %in% year, member == member_input, lowlevel == TRUE), 
               aes(x = date, y = drainage), color = "purple") +
    labs(title = "Abflussmenge für einen Pegel über die Zeit",
       y = "Abflussmenge (in m³/s, Log-Skala)",
       x = "Datum") +
    theme(text = element_text(size = 20)) +
    scale_color_brewer(palette = "Dark2")
}

plot_single_drainage(hydro_total %>% filter(name_waterlevel == "10304"), year = 2000:2003)
# -> Daten sind saisonal
ggsave("Plots/Chris/analysis/single_drainage_Isar-Mittenwald_2000-2003.png", width = 12, height = 9)

plot_drainage(hydro_total, year = 2000)
ggsave("Plots/Chris/analysis/drainage_total_2000_kbe.png", width = 12, height = 9)

plot_drainage_unmarked(hydro_total, year = 2000)
ggsave("Plots/Chris/analysis/drainage_total_unmarked_2000_kbe.png", width = 12, height = 9)

# Visualisierung der 10 Member für 10304, Jahr 2000
hydro_total_10304 <- hydro_total %>% filter(name_waterlevel == "10304")
ggplot(data = subset(hydro_total_10304, YY %in% 2000), mapping = aes(x = date, y = drainage, color = member, alpha = 0.01)) +
  geom_line() +
  labs(title = "Abfluss des Pegels Isar Mittenwald im Jahr 2000 getrennt \nnach Simulationsreihen",
       y = "Abflussmenge (in m³/s, Log-Skala)",
       x = "Datum")  +
    scale_y_continuous(trans='log10')+
    theme(text = element_text(size = 20)) +
    scale_color_brewer(palette = "Paired")
ggsave("Plots/Chris/analysis/all_member_drainage_Isar-Mittenwald_2000.png", width = 12, height = 9)

# Tables for futher plots
table_days_low_flow <- hydro_total %>%
  group_by(waterlevel, YY) %>%
  count(lowlevel, name = "number_low_flows") %>%
  filter(lowlevel == TRUE) %>%
  mutate(avg_number_low_flows = number_low_flows / length(unique(hydro_total$member)))

table_days_low_flow_by_hydro_year <- hydro_total %>%
  group_by(hydro_year, waterlevel, YY) %>%
  count(lowlevel, name = "number_low_flows") %>%
  filter(lowlevel == TRUE) %>%
  mutate(avg_number_low_flows = number_low_flows / length(unique(hydro_total$member)))

# Jährliche durchschnittliche Anzahl von Niedrigwasser-Tagen für die 3 Wasserstände
ggplot(data = table_days_low_flow, mapping = aes(x = YY, y = avg_number_low_flows, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliche durchschnittliche Anzahl von Niedrigwasser-Tagen \nim Zeitverlauf (aller Member)",
       y = "Jährliche durchschnittliche Anzahl von Niedrigwasser-Tagen",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2")
ggsave("Plots/Chris/analysis/yearly_avg_number_of_low_flow_days.png", width = 12, height = 9)

# Jährliche durchschnittliche Anzahl von Niedrigwasser-Tagen für die 3 Wasserstände je Jahreszeit
ggplot(data = table_days_low_flow_by_hydro_year, mapping = aes(x = YY, y = avg_number_low_flows, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(
    title = "Jährliche durchschnittliche Anzahl von Niedrigwasser-Tagen je Jahreszeit \nim Zeitverlauf (aller Member)",
       y = "Jährliche durchschnittliche Anzahl von Niedrigwasser-Tagen",
       x = "Datum") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter"))) +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2")
ggsave("Plots/Chris/analysis/yearly_avg_number_of_low_flow_days_per_season.png", width = 12, height = 9)



```
# Analysis of drivers ----
```{r Analysis of drivers}
# Read tables
table_yearly_avg_min_groundwaterdepth <- 
  readRDS("data/tables/driver_analysis/table_yearly_avg_min_groundwaterdepth.RDS")
table_yearly_avg_min_soilwater <- readRDS("data/tables/driver_analysis/table_yearly_avg_min_soilwater.RDS")
table_yearly_avg_min_snowstorage <- readRDS("data/tables/driver_analysis/table_yearly_avg_min_snowstorage.RDS")
table_yearly_avg_min_airtmp <- readRDS("data/tables/driver_analysis/table_yearly_avg_min_airtmp.RDS")
table_yearly_avg_max_precip <- readRDS("data/tables/driver_analysis/table_yearly_avg_max_precip.RDS")
table_yearly_avg_max_glorad <- readRDS("data/tables/driver_analysis/table_yearly_avg_max_glorad.RDS")
table_yearly_avg_max_relhum <- readRDS("data/tables/driver_analysis/table_yearly_avg_max_relhum.RDS")
table_yearly_avg_max_infiltration <- readRDS("data/tables/driver_analysis/table_yearly_avg_max_infiltration.RDS")


# Tables for futher plots
table_yearly_max_prec_by_member <- hydro_total %>%
  group_by(waterlevel, YY, hydro_year, member) %>%
  summarise(max_precip = max(precip))

table_yearly_avg_max_prec <- hydro_total %>%
  group_by(waterlevel, YY, hydro_year, member) %>%
  summarise(max_precip = max(precip)) %>%
  group_by(waterlevel, YY, hydro_year) %>%
  summarise(avg_max_precip = mean(max_precip)) 

table_yearly_avg_max_airtmp <- hydro_total %>%
  group_by(waterlevel, YY, hydro_year, member) %>%
  summarise(max_airtmp = max(airtmp)) %>% 
  group_by(waterlevel, YY, hydro_year) %>%
  summarise(avg_max_airtemp = mean(max_airtmp)) 

table_yearly_avg_max_snowstorage <- hydro_total %>%
  group_by(waterlevel, YY, hydro_year, member) %>%
  summarise(max_snowstorage = max(snowstorage)) %>% 
  group_by(waterlevel, YY, hydro_year) %>%
  summarise(avg_max_snowstorage = mean(max_snowstorage)) 

table_yearly_avg_max_soilwater <- hydro_total %>%
  group_by(waterlevel, YY, hydro_year, member) %>%
  summarise(max_soilwater = max(soilwater)) %>% 
  group_by(waterlevel, YY, hydro_year) %>%
  summarise(avg_max_soilwater = mean(max_soilwater))

# Model drivers: avg and groundwaterdepth


# Plots
ggplot(data = table_yearly_avg_max_prec, mapping = aes(x = YY, y = avg_max_precip, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Maximum an täglichem Niederschlag \nje Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Maximum an täglichem Niederschlag \n(in mm/24h)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/yearly_avg_max_prec.png", width = 12, height = 9)

ggplot(data = table_yearly_avg_max_airtmp, mapping = aes(x = YY, y = avg_max_airtemp, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Maximum an täglicher Lufttemperatur \nje Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Maximum an täglicher Lufttemperatur \n(in °C)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/yearly_avg_max_airtemp.png", width = 12, height = 9)

ggplot(data = table_yearly_avg_max_snowstorage, mapping = aes(x = YY, y = avg_max_snowstorage, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Maximum an täglichem Schneespeicher \nje Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Maximum an täglichem Schneespeicher \n(in mm)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/yearly_avg_max_snowstorage.png", width = 12, height = 9)

ggplot(data = table_yearly_avg_max_soilwater, mapping = aes(x = YY, y = avg_max_soilwater, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Maximum an täglicher Oberflächennahe \nrelative Bodenfeuchte je Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Maximum an täglicher \noberflächennaher relative Bodenfeuchte (in %)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/yearly_avg_max_soilwater.png", width = 12, height = 9)

# Plots for model drivers

ggplot(data = table_yearly_avg_min_groundwaterdepth, mapping = aes(x = YY, y = avg_min_groundwaterdepth, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Minimum des Grundwasserstandes \nje Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Minimum des täglichem \nGrundwasserstandes (in m unter Oberfläche)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/driver_plots/yearly_avg_min_groundwaterdepth.png", width = 12, height = 9)

ggplot(data = table_yearly_avg_min_soilwater, mapping = aes(x = YY, y = avg_min_soilwater, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Minimum der täglichen oberflächennahen \nrelativen Bodenfeuchte je Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Minimum der täglichen oberflächennahen \nrelativen Bodenfeuchte (in %)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/driver_plots/yearly_avg_min_soilwater.png", width = 12, height = 9)

ggplot(data = table_yearly_avg_min_snowstorage, mapping = aes(x = YY, y = avg_min_snowstorage, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Minimum an täglichem Schneespeicher \nje Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Minimum an täglichem Schneespeicher \n(in mm)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/driver_plots/yearly_avg_min_snowstorage.png", width = 12, height = 9)

ggplot(data = table_yearly_avg_min_airtmp, mapping = aes(x = YY, y = avg_min_airtmp, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Minimum an täglicher Lufttemperatur \nje Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Minimum an täglicher Lufttemperatur \n(in °C)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/driver_plots/yearly_avg_min_airtmp.png", width = 12, height = 9)

ggplot(data = table_yearly_avg_max_precip, mapping = aes(x = YY, y = avg_max_precip, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Maximum an täglichem Niederschlag \nje Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Maximum an täglichem Niederschlag \n(in mm/24h)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/driver_plots/yearly_avg_max_precip.png", width = 12, height = 9)

ggplot(data = table_yearly_avg_max_glorad, mapping = aes(x = YY, y = avg_max_glorad, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Maximum an täglich einfallender \nkurzwelligen Strahlung je Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Maximum an täglich einfallender \nkurzwelligen Strahlung (in Wh/m²)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/driver_plots/yearly_avg_max_glorad.png", width = 12, height = 9)

ggplot(data = table_yearly_avg_max_relhum, mapping = aes(x = YY, y = avg_max_relhum, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Maximum an täglicher relativer \nLuftfeuchtigkeit je Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Maximum an täglicher relativer \nLuftfeuchtigkeit (in %)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/driver_plots/yearly_avg_max_relhum.png", width = 12, height = 9)

ggplot(data = table_yearly_avg_max_infiltration, mapping = aes(x = YY, y = avg_max_infiltration, color = waterlevel)) +
  geom_line() +
  geom_point() +
  labs(title = "Jährliches durchschnittliches Maximum an täglicher Versickerung \nje Jahreszeit im Zeitverlauf (aller Member)",
       y = "Jährliches durchschnittliches Maximum an täglicher Versickerung \n(in mm/24h)",
       x = "Datum") +
  guides(color=guide_legend(title="Pegel")) +
    theme(text = element_text(size = 20),
          legend.position="bottom") +
    scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(hydro_year), labeller = as_labeller(c(`summer` = "Sommer", `winter` = "Winter")))
ggsave("Plots/Chris/analysis/driver_plots/yearly_avg_max_infiltration.png", width = 12, height = 9)
```


## Chris: Extreme value analysis
```{r Chris: Global extreme value analysis}
### Quantile: -------------------------------------------------------------

# Define the quantiles for which you want to calculate the percentage of the target variable
quantiles <- c(0.05, 0.1, 0.25, 0.75, 0.9, 0.95)

# Create an empty matrix to store the percentage values
quantile_percents <- matrix(nrow=length(quantiles), ncol=ncol(hydro_total[, 15:22]))

# Loop over each predictor variable and calculate the percentage of the target variable for each quantile
for (i in 15:22) {
  # Get the quantile values for the current predictor variable
  quantile_values <- quantile(hydro_total[, i], probs=quantiles)
  for (j in 1:length(quantiles)) {
    # Calculate the percentage of the target variable for the current quantile and predictor variable
    quantile_percents[j, i - 14] <- 
      sum((hydro_total[hydro_total[, i] >= quantile_values[j], "lowlevel"]) == "TRUE")/sum(
        hydro_total[, "lowlevel"] == "TRUE")
  }
}

# Print the results
colnames(quantile_percents) <- colnames(hydro_total)[15:22]
rownames(quantile_percents) <- paste0(quantiles*100, "%")
quantile_percents <- round(quantile_percents, digits = 3)
quantile_percents

# Save the table as a PDF file
pdf("tables/quantile_percents.pdf", height=4, width=10)
grid.table(quantile_percents)
dev.off()

### Ranges: ---------------------------------------------------------------
# Define the quantiles for which you want to calculate the percentage of the target variable
quantiles_ranges <- seq(0, 1, by=0.1)

# Create an empty matrix to store the percentage values
quantile_percents_ranges <- matrix(nrow=length(quantiles_ranges)-1, ncol=ncol(hydro_total[, 15:22]))

# Loop over each predictor variable and calculate the percentage of the target variable for each range
for (i in 15:22) {
  # Get the quantile values for the current predictor variable
  quantile_values <- quantile(hydro_total[, i], probs=quantiles_ranges)
  for (j in 1:(length(quantiles_ranges)-1)) {
    # Calculate the percentage of the target variable for the current range and predictor variable
    quantile_percents_ranges[j, i - 14] <- 
      sum((hydro_total[hydro_total[, i] >= quantile_values[j] & hydro_total[, i] < quantile_values[j+1], "lowlevel"]) == "TRUE")/sum((hydro_total[, "lowlevel"]) == "TRUE")
  }
}

# Print the results
colnames(quantile_percents_ranges) <- colnames(hydro_total)[15:22]
rownames(quantile_percents_ranges) <- 
  paste0(quantiles_ranges[-length(quantiles_ranges)]*100, "% - ", quantiles_ranges[-1]*100, "%")
quantile_percents_ranges <- round(quantile_percents_ranges, digits = 3)
quantile_percents_ranges

# Save the table as a PDF file
pdf("tables/quantile_percents_ranges.pdf", height=4, width=10)
grid.table(quantile_percents_ranges)
dev.off()



```

```{r Chris local extreme value analysis}
pegel_list <- list(hydro_summer_20203, hydro_summer_11502, hydro_summer_10304, 
               hydro_winter_20203, hydro_winter_11502, hydro_winter_10304)
pegel_names <- c("hydro_summer_20203", "hydro_summer_11502", "hydro_summer_10304", 
               "hydro_winter_20203", "hydro_winter_11502", "hydro_winter_10304")
counter_names <- 1

# Define the quantiles for which you want to calculate the percentage of the target variable
quantiles_ranges <- seq(0, 1, by=0.1)
colum_names <- c("avg_precip", "avg_airtmp", "avg_glorad", "avg_relhum", "avg_soilwater", "avg_snowstorage", 
                 "groundwaterdepth", "avg_infiltration", "max_precip")#, "avg_snowstorage_drain")
colum_names_long <- c("Mittlerer\nNiederschlag", "Mittlere\nLufttemperatur", 
                      "Mittlere\neinfallende\nkurzwellige Strahlung", "Mittlere\nrelative\nLuftfeuchte", 
                      "Mittlere\noberflächennahe\nrelative\nBodenfeuchte", "Mittlerer\nSchneespeicher", 
                      "Grundwasserstand", "Mittlere\nVersickerung", "Maximaler\nNiederschlag")
for (pegel in pegel_list) {
  i = 1
  
  # Create an empty matrix to store the percentage values
  quantile_percents_ranges <- matrix(nrow=length(quantiles_ranges)-1, ncol=length(colum_names))

    # Loop over each predictor variable and calculate the percentage of the target variable for each range
  for (col_name in colum_names) {
    # Get the quantile values for the current predictor variable
    quantile_values <- quantile(pegel[, col_name], probs=quantiles_ranges, na.rm = TRUE)
    for (j in 1:(length(quantiles_ranges)-1)) {
      # Calculate the percentage of the target variable for the current range and predictor variable
      quantile_percents_ranges[j, i] <- 
        sum((pegel[pegel[, col_name] >= quantile_values[j] & pegel[, col_name] < quantile_values[j+1],
                         "lowlevel"]) == "TRUE", na.rm = TRUE)/sum((pegel[, "lowlevel"]) == "TRUE", na.rm = TRUE)
    }
    i = i + 1
  }
  
  # Print the results
  colnames(quantile_percents_ranges) <- colum_names_long
  rownames(quantile_percents_ranges) <- paste0(
    quantiles_ranges[-length(quantiles_ranges)]*100, "% - ", quantiles_ranges[-1]*100, "%")
  quantile_percents_ranges <- round(quantile_percents_ranges, digits = 3)
  
  # Save as table
  saveRDS(object = quantile_percents_ranges, 
          file = paste0("data/tables/extreme_values/quantile_percents_ranges_", pegel_names[counter_names], ".RDS"))
  
  # Save the table as a PDF file
  pdf(paste0("tables/extreme_values/quantile_percents_ranges_", pegel_names[counter_names], ".pdf"), height=4, width=16)
  grid.table(quantile_percents_ranges)
  dev.off()
  counter_names = counter_names + 1
}
```

```{r Extreme value analysis subset}
pegel <- hydro_summer_11502
pegel_names <- "hydro_summer_11502_subset"
counter_names <- 1

# Define the quantiles for which you want to calculate the percentage of the target variable
quantiles_ranges <- seq(0, 1, by=0.1)
colum_names <- c("avg_airtmp", "avg_soilwater")
colum_names_long <- c("Mittlere\nLufttemperatur", "Mittlere\noberflächennahe\nrelative\nBodenfeuchte")
  i = 1
  
  # Create an empty matrix to store the percentage values
  quantile_percents_ranges <- matrix(nrow=length(quantiles_ranges)-1, ncol=length(colum_names))

    # Loop over each predictor variable and calculate the percentage of the target variable for each range
  for (col_name in colum_names) {
    # Get the quantile values for the current predictor variable
    quantile_values <- quantile(pegel[, col_name], probs=quantiles_ranges, na.rm = TRUE)
    for (j in 1:(length(quantiles_ranges)-1)) {
      # Calculate the percentage of the target variable for the current range and predictor variable
      quantile_percents_ranges[j, i] <- 
        sum((pegel[pegel[, col_name] >= quantile_values[j] & pegel[, col_name] < quantile_values[j+1],
                         "lowlevel"]) == "TRUE", na.rm = TRUE)/sum((pegel[, "lowlevel"]) == "TRUE", na.rm = TRUE)
    }
    i = i + 1
  }
  
  # Print the results
  colnames(quantile_percents_ranges) <- colum_names_long
  rownames(quantile_percents_ranges) <- paste0(
    quantiles_ranges[-length(quantiles_ranges)]*100, "% - ", quantiles_ranges[-1]*100, "%")
  quantile_percents_ranges <- round(quantile_percents_ranges, digits = 3)
  
  # Save as table
  saveRDS(object = quantile_percents_ranges, 
          file = paste0("data/tables/extreme_values/quantile_percents_ranges_", pegel_names[counter_names], ".RDS"))
  
  # Save the table as a PDF file
  pdf(paste0("tables/extreme_values/quantile_percents_ranges_", pegel_names[counter_names], ".pdf"), height=4, width=16)
  grid.table(quantile_percents_ranges)
  dev.off()
  counter_names = counter_names + 1

```


## Stationary analysis
```{r}
library(tseries)

stationary_drainage_mean_p <- function(waterlevel, season) {
  name <- paste0("stationar_", season, "_", waterlevel, "_")
  member_names <- levels(hydro_total$member)
  for (i in 1:length(member_names)) {
    assign(paste0(name, i),
           kpss.test(get(paste0("hydro_", season, "_", waterlevel, "_", member_names[i]))$drainage)$p.value)
  }
  mean(get(paste0(name, "1")), get(paste0(name, "2")), get(paste0(name, "3")), get(paste0(name, "4")),
       get(paste0(name, "5")), get(paste0(name, "6")), get(paste0(name, "7")), get(paste0(name, "8")),
       get(paste0(name, "9")), get(paste0(name, "10")))
}

stationary_drainage_mean_p("10304", "winter")
stationary_drainage_mean_p("10304", "summer")
stationary_drainage_mean_p("11502", "winter")
stationary_drainage_mean_p("11502", "summer")
stationary_drainage_mean_p("20203", "winter")
stationary_drainage_mean_p("20203", "summer")
# nur die Zeitreihe von drainage für 20203/winter ist zum Signifikanzniveau 0.05 signifikant nicht-stationär;
# bei allen anderen Zeitreihen von drainage kann die Annahme der Stationarität nicht abgelehnt werden

stationary_lowlevel_mean_p <- function(waterlevel, season) {
  name <- paste0("stationar_", season, "_", waterlevel, "_")
  member_names <- levels(hydro_total$member)
  for (i in 1:length(member_names)) {
    assign(paste0(name, i),
           kpss.test(get(paste0("hydro_", season, "_", waterlevel, "_", member_names[i]))$lowlevel)$p.value)
  }
  mean(get(paste0(name, "1")), get(paste0(name, "2")), get(paste0(name, "3")), get(paste0(name, "4")),
       get(paste0(name, "5")), get(paste0(name, "6")), get(paste0(name, "7")), get(paste0(name, "8")),
       get(paste0(name, "9")), get(paste0(name, "10")))
}

stationary_lowlevel_mean_p("10304", "winter")
stationary_lowlevel_mean_p("10304", "summer")
stationary_lowlevel_mean_p("11502", "winter")
stationary_lowlevel_mean_p("11502", "summer")
stationary_lowlevel_mean_p("20203", "winter")
stationary_lowlevel_mean_p("20203", "summer")
# die Zeitreihen von lowlevel für 10304/winter, 20203/winter, 20203/summer sind zum Signifikanzniveau 0.05 signifikant nicht-stationär:
# für die Zeitreihen von lowlevel für 10304/summer, 11502/winter, 11502/summer kann die Annahme der Stationarität nicht abgelehnt werden

```


# Several member visualization
```{r}
hydro_summer_10304 %>% 
  filter(YY %in% seq(1990, 2020, 5)) %>% 
  ggplot(aes(x = date, y = avg_airtmp,col = member)) +
  geom_line(alpha = 0.2)+
  facet_wrap(~ YY, scales = "free_x")
```

```{r}
plot_ts_facet <- function(variable, data, period = 5, ...) {
  data %>% 
  filter(YY %in% seq(1990, 2020, period)) %>% 
  ggplot(aes_string(x = "date", y = variable, col = "member")) +
  geom_path(alpha = 0.2)+
  facet_wrap(~ YY, scales = "free_x")+
    labs(...)
}
```
 
# Summer
```{r}
plot_ts_facet("avg_airtmp", hydro_summer_10304, period = 5, title = 10304)
plot_ts_facet("avg_airtmp", hydro_summer_11502, period = 5, title = 11502)
plot_ts_facet("avg_airtmp", hydro_summer_20203, period = 5, title = 20203)

plot_ts_facet("avg_relhum", hydro_summer_10304, period = 5, title = 10304)
plot_ts_facet("avg_relhum", hydro_summer_11502, period = 5, title = 11502)
plot_ts_facet("avg_relhum", hydro_summer_20203, period = 5, title = 20203)

plot_ts_facet("avg_glorad", hydro_summer_10304, period = 5, title = 10304)
plot_ts_facet("avg_glorad", hydro_summer_11502, period = 5, title = 11502)
plot_ts_facet("avg_glorad", hydro_summer_20203, period = 5, title = 20203)
```



### Auto corr
```{r}
library(car)

## drainage
# autocorrelation function, mean acf for drainage (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_acf_drainage_", k, "_", i), data.frame(Names = 1:38))
    for (j in levels(hydro_total$member)) {
      assign(paste0("acf_drainage_", k, "_", i, "_", j), as.data.frame(acf(get(paste0("hydro_", k, "_", i, "_", j))$drainage, plot = FALSE)$acf))
      assign(paste0("df_acf_drainage_", k, "_", i), bind_cols(get(paste0("df_acf_drainage_", k, "_", i)), get(paste0("acf_drainage_", k, "_", i, "_", j))))
    }
    assign(paste0("acf_mean_drainage_", k, "_", i), rowMeans(get(paste0("df_acf_drainage_", k, "_", i))[, 2:11]))
  }
}
# -> Drainage ist in allen Teildatensätzen positiv autokorreliert -> Modell, das diese Autokorrelation berücksichtigt, nötig

# ACF-Plots for drainage (across all members) separated by each waterlevel and season
plot(x = 0:37, y = acf_mean_drainage_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Abflusses für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_drainage_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Abflusses für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_drainage_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Abflusses für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_drainage_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Abflusses für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_drainage_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Abflusses für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_drainage_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Abflusses für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

# Partial autocorrelation function, mean pacf for drainage (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_pacf_drainage_", k, "_", i), data.frame(Names = 1:37))
    for (j in levels(hydro_total$member)) {
      assign(paste0("pacf_drainage_", k, "_", i, "_", j),
             as.data.frame(pacf(get(paste0("hydro_", k, "_", i, "_", j))$drainage, plot = FALSE)$acf))
      assign(paste0("df_pacf_drainage_", k, "_", i),
             bind_cols(get(paste0("df_pacf_drainage_", k, "_", i)), get(paste0("pacf_drainage_", k, "_", i, "_", j))))
    }
    assign(paste0("pacf_mean_drainage_", k, "_", i), rowMeans(get(paste0("df_pacf_drainage_", k, "_", i))[, 2:11]))
  }
}
# 10304, Winter: Autokorrelationskoeffizient bis lag = 2 (durchgängig) signifikant
# 10304, Summer: Autokorrelationskoeffizient bis lag = 3 (durchgängig) signifikant
# 11502, Winter: Autokorrelationskoeffizient bis lag = 5 (nicht durchgängig), bis lag = 3 (durchgängig) signifikant
# 11502, Summer: Autokorrelationskoeffizient bis lag = 3 (nicht durchgängig), bis lag = 1 (durchgängig) signifikant
# 20203, Winter: Autokorrelationskoeffizient bis lag = 5 (nicht durchgängig), bis lag = 3 (durchgängig) signifikant
# 20203, Summer: Autokorrelationskoeffizient bis lag = 5 (nicht durchgängig), bis lag = 3 (durchgängig) signifikant

# PACF-Plots for drainage (across all members) separated by each waterlevel and season
plot(x = 1:37, y = pacf_mean_drainage_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Abflusses für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_drainage_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Abflusses für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_drainage_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Abflusses für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_drainage_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Abflusses für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_drainage_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Abflusses für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_drainage_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Abflusses für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")


## lowlevel
# autocorrelation function, mean acf for lowlevel (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_acf_lowlevel_", k, "_", i), data.frame(Names = 1:38))
    for (j in levels(hydro_total$member)) {
      assign(paste0("acf_lowlevel_", k, "_", i, "_", j), as.data.frame(acf(get(paste0("hydro_", k, "_", i, "_", j))$lowlevel, plot = FALSE)$acf))
      assign(paste0("df_acf_lowlevel_", k, "_", i), bind_cols(get(paste0("df_acf_lowlevel_", k, "_", i)), get(paste0("acf_lowlevel_", k, "_", i, "_", j))))
    }
    assign(paste0("acf_mean_lowlevel_", k, "_", i), rowMeans(get(paste0("df_acf_lowlevel_", k, "_", i))[, 2:11]))
  }
}
# -> Lowlevel ist in allen Teildatensätzen positiv autokorreliert -> Modell, das diese Autokorrelation berücksichtigt, nötig

# ACF-Plots for lowlevel (across all members) separated by each waterlevel and season
plot(x = 0:37, y = acf_mean_lowlevel_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. von Niedrigwasserevents für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_lowlevel_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. von Niedrigwasserevents für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_lowlevel_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. von Niedrigwasserevents für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_lowlevel_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. von Niedrigwasserevents für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_lowlevel_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. von Niedrigwasserevents für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_lowlevel_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. von Niedrigwasserevents für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

# Partial autocorrelation function, mean pacf for lowlevel (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_pacf_lowlevel_", k, "_", i), data.frame(Names = 1:37))
    for (j in levels(hydro_total$member)) {
      assign(paste0("pacf_lowlevel_", k, "_", i, "_", j),
             as.data.frame(pacf(get(paste0("hydro_", k, "_", i, "_", j))$lowlevel, plot = FALSE)$acf))
      assign(paste0("df_pacf_lowlevel_", k, "_", i),
             bind_cols(get(paste0("df_pacf_lowlevel_", k, "_", i)), get(paste0("pacf_lowlevel_", k, "_", i, "_", j))))
    }
    assign(paste0("pacf_mean_lowlevel_", k, "_", i), rowMeans(get(paste0("df_pacf_lowlevel_", k, "_", i))[, 2:11]))
  }
}
# 10304, Winter: Autokorrelationskoeffizient bis lag = 4 (nicht durchgängig), bis lag = 1 (durchgängig) signifikant
# 10304, Summer: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 11502, Winter: Autokorrelationskoeffizient bis lag = 4 (nicht durchgängig), bis lag = 1 (durchgängig) signifikant
# 11502, Summer: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 20203, Winter: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 20203, Summer: Autokorrelationskoeffizient bis lag = 5 (nicht durchgängig), bis lag = 1 (durchgängig) signifikant

# PACF-Plots for lowlevel (across all members) seperated by each waterlevel and season
plot(x = 1:37, y = pacf_mean_lowlevel_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. v. Niedrigwasserevents für Isar Mittenw., Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_lowlevel_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. von Niedrigwasserevents für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_lowlevel_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. von Niedrigwasserevents für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_lowlevel_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. von Niedrigwasserevents für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_lowlevel_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. von Niedrigwasserevents für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_lowlevel_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. von Niedrigwasserevents für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

# DurbinWatsonTest für GAM auf Datensatz hydro_winter_10304_kbe
durbinWatsonTest(gam_1)
# p < 0.05, Wert der Teststatistik: 0.3839 -> Es scheint positive Autokorrelation vorzuliegen

### Driver
## Precip
# autocorrelation function, mean acf for precip (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_acf_precip_", k, "_", i), data.frame(Names = 1:38))
    for (j in levels(hydro_total$member)) {
      assign(paste0("acf_precip_", k, "_", i, "_", j),
             as.data.frame(acf(get(paste0("hydro_", k, "_", i, "_", j))$precip, plot = FALSE)$acf))
      assign(paste0("df_acf_precip_", k, "_", i),
             bind_cols(get(paste0("df_acf_precip_", k, "_", i)), get(paste0("acf_precip_", k, "_", i, "_", j))))
    }
    assign(paste0("acf_mean_precip_", k, "_", i), rowMeans(get(paste0("df_acf_precip_", k, "_", i))[, 2:11]))
  }
}
# -> Precip ist in allen Teildatensätzen leicht positiv autokorreliert -> Modell, das diese Autokorrelation berücksichtigt, wäre von Vorteil

# ACF-Plots for precip (across all members) separated by each waterlevel and season
plot(x = 0:37, y = acf_mean_precip_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Niederschlags für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_precip_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Niederschlags für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_precip_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Niederschlags für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_precip_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Niederschlags für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_precip_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Niederschlags für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_precip_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Niederschlags für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

# partial autocorrelation function, mean pacf for precip (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_pacf_precip_", k, "_", i), data.frame(Names = 1:37))
    for (j in levels(hydro_total$member)) {
      assign(paste0("pacf_precip_", k, "_", i, "_", j),
             as.data.frame(pacf(get(paste0("hydro_", k, "_", i, "_", j))$precip, plot = FALSE)$acf))
      assign(paste0("df_pacf_precip_", k, "_", i),
             bind_cols(get(paste0("df_pacf_precip_", k, "_", i)), get(paste0("pacf_precip_", k, "_", i, "_", j))))
    }
    assign(paste0("pacf_mean_precip_", k, "_", i), rowMeans(get(paste0("df_pacf_precip_", k, "_", i))[, 2:11]))
  }
}
# 10304, Winter: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 10304, Summer: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 11502, Winter: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 11502, Summer: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 20203, Winter: Autokorrelationskoeffizient bis lag = 2 (durchgängig) signifikant
# 20203, Summer: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant

# PACF-Plots for precip (across all members) seperated by each waterlevel and season
plot(x = 1:37, y = pacf_mean_precip_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Niederschlags für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_precip_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Niederschlags für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_precip_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Niederschlags für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_precip_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Niederschlags für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_precip_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Niederschlags für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_precip_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Niederschlags für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")


## Airtmp
# autocorrelation function, mean acf for airtmp (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_acf_airtmp_", k, "_", i), data.frame(Names = 1:38))
    for (j in levels(hydro_total$member)) {
      assign(paste0("acf_airtmp_", k, "_", i, "_", j),
             as.data.frame(acf(get(paste0("hydro_", k, "_", i, "_", j))$airtmp, plot = FALSE)$acf))
      assign(paste0("df_acf_airtmp_", k, "_", i),
             bind_cols(get(paste0("df_acf_airtmp_", k, "_", i)), get(paste0("acf_airtmp_", k, "_", i, "_", j))))
    }
    assign(paste0("acf_mean_airtmp_", k, "_", i), rowMeans(get(paste0("df_acf_airtmp_", k, "_", i))[, 2:11]))
  }
}
# -> Airtmp ist in allen Teildatensätzen positiv autokorreliert -> Modell, das diese Autokorrelation berücksichtigt, nötig

# ACF-Plots for airtmp (across all members) separated by each waterlevel and season
plot(x = 0:37, y = acf_mean_airtmp_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Autokor. zw. Beob. der Temperatur für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_airtmp_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Autokor. zw. Beob. der Temperatur für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_airtmp_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Autokor. zw. Beob. der Temperatur für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_airtmp_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Autokor. zw. Beob. der Temperatur für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_airtmp_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Autokor. zw. Beob. der Temperatur für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_airtmp_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Autokor. zw. Beob. der Temperatur für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

# partial autocorrelation function, mean pacf for airtmp (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_pacf_airtmp_", k, "_", i), data.frame(Names = 1:37))
    for (j in levels(hydro_total$member)) {
      assign(paste0("pacf_airtmp_", k, "_", i, "_", j),
             as.data.frame(pacf(get(paste0("hydro_", k, "_", i, "_", j))$airtmp, plot = FALSE)$acf))
      assign(paste0("df_pacf_airtmp_", k, "_", i),
             bind_cols(get(paste0("df_pacf_airtmp_", k, "_", i)), get(paste0("pacf_airtmp_", k, "_", i, "_", j))))
    }
    assign(paste0("pacf_mean_airtmp_", k, "_", i), rowMeans(get(paste0("df_pacf_airtmp_", k, "_", i))[, 2:11]))
  }
}
# 10304, Winter: Autokorrelationskoeffizient bis lag = 5 (nicht durchgängig), bis lag = 3 (durchgängig) signifikant
# 10304, Summer: Autokorrelationskoeffizient bis lag = 7 (nicht durchgängig), bis lag = 3 (durchgängig) signifikant
# 11502, Winter: Autokorrelationskoeffizient bis lag = 5 (nicht durchgängig), bis lag = 3 (durchgängig) signifikant
# 11502, Summer: Autokorrelationskoeffizient bis lag = 7 (nicht durchgängig), bis lag = 3 (durchgängig) signifikant
# 20203, Winter: Autokorrelationskoeffizient bis lag = 5 (durchgängig) signifikant
# 20203, Summer: Autokorrelationskoeffizient bis lag = 6 (durchgängig) signifikant

# PACF-Plots for airtmp (across all members) seperated by each waterlevel and season
plot(x = 1:37, y = pacf_mean_airtmp_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Temperatur für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_airtmp_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Temperatur für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_airtmp_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Temperatur für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_airtmp_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Temperatur für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_airtmp_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Temperatur für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_airtmp_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Temperatur für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

## Glorad
# autocorrelation function, mean acf for glorad (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_acf_glorad_", k, "_", i), data.frame(Names = 1:38))
    for (j in levels(hydro_total$member)) {
      assign(paste0("acf_glorad_", k, "_", i, "_", j),
             as.data.frame(acf(get(paste0("hydro_", k, "_", i, "_", j))$glorad, plot = FALSE)$acf))
      assign(paste0("df_acf_glorad_", k, "_", i),
             bind_cols(get(paste0("df_acf_glorad_", k, "_", i)), get(paste0("acf_glorad_", k, "_", i, "_", j))))
    }
    assign(paste0("acf_mean_glorad_", k, "_", i), rowMeans(get(paste0("df_acf_glorad_", k, "_", i))[, 2:11]))
  }
}
# -> Glorad ist in allen Teildatensätzen positiv autokorreliert -> Modell, das diese Autokorrelation berücksichtigt, nötig

# ACF-Plots for glorad (across all members) separated by each waterlevel and season
plot(x = 0:37, y = acf_mean_glorad_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der einfall. kurzw. Strahlung für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_glorad_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der einfall. kurzw. Strahlung für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_glorad_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der einfall. kurzw. Strahlung für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_glorad_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der einfall. kurzw. Strahlung für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_glorad_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der einfall. kurzw. Strahlung für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_glorad_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der einfall. kurzw. Strahlung für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

# partial autocorrelation function, mean pacf for glorad (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_pacf_glorad_", k, "_", i), data.frame(Names = 1:37))
    for (j in levels(hydro_total$member)) {
      assign(paste0("pacf_glorad_", k, "_", i, "_", j),
             as.data.frame(pacf(get(paste0("hydro_", k, "_", i, "_", j))$glorad, plot = FALSE)$acf))
      assign(paste0("df_pacf_glorad_", k, "_", i),
             bind_cols(get(paste0("df_pacf_glorad_", k, "_", i)), get(paste0("pacf_glorad_", k, "_", i, "_", j))))
    }
    assign(paste0("pacf_mean_glorad_", k, "_", i), rowMeans(get(paste0("df_pacf_glorad_", k, "_", i))[, 2:11]))
  }
}
# 10304, Winter: Autokorrelationskoeffizient bis lag = 7 (durchgängig) signifikant
# 10304, Summer: Autokorrelationskoeffizient bis lag = 9 (durchgängig) signifikant
# 11502, Winter: Autokorrelationskoeffizient bis lag = 7 (durchgängig) signifikant
# 11502, Summer: Autokorrelationskoeffizient bis lag = 8 (durchgängig) signifikant
# 20203, Winter: Autokorrelationskoeffizient bis lag = 6 (durchgängig) signifikant
# 20203, Summer: Autokorrelationskoeffizient bis lag = 7 (durchgängig) signifikant

# PACF-Plots for glorad (across all members) seperated by each waterlevel and season
plot(x = 1:37, y = pacf_mean_glorad_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der einfall. kurzw. Strahlung für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_glorad_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der einfall. kurzw. Strahlung für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_glorad_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der einfall. kurzw. Strahlung für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_glorad_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der einfall. kurzw. Strahlung für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_glorad_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der einfall. kurzw. Strahlung für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_glorad_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der einfall. kurzw. Strahlung für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")


## Relhum
# autocorrelation function, mean acf for relhum (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_acf_relhum_", k, "_", i), data.frame(Names = 1:38))
    for (j in levels(hydro_total$member)) {
      assign(paste0("acf_relhum_", k, "_", i, "_", j),
             as.data.frame(acf(get(paste0("hydro_", k, "_", i, "_", j))$relhum, plot = FALSE)$acf))
      assign(paste0("df_acf_relhum_", k, "_", i),
             bind_cols(get(paste0("df_acf_relhum_", k, "_", i)), get(paste0("acf_relhum_", k, "_", i, "_", j))))
    }
    assign(paste0("acf_mean_relhum_", k, "_", i), rowMeans(get(paste0("df_acf_relhum_", k, "_", i))[, 2:11]))
  }
}
# -> Relhum ist in allen Teildatensätzen positiv autokorreliert (vor allem für Pegel 20203) -> Modell, das diese Autokorrelation berücksichtigt, nötig

# ACF-Plots for relhum (across all members) separated by each waterlevel and season
plot(x = 0:37, y = acf_mean_relhum_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der relativen Luftfeuchte für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_relhum_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der relativen Luftfeuchte für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_relhum_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der relativen Luftfeuchte für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_relhum_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der relativen Luftfeuchte für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_relhum_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der relativen Luftfeuchte für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_relhum_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der relativen Luftfeuchte für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

# partial autocorrelation function, mean pacf for relhum (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_pacf_relhum_", k, "_", i), data.frame(Names = 1:37))
    for (j in levels(hydro_total$member)) {
      assign(paste0("pacf_relhum_", k, "_", i, "_", j),
             as.data.frame(pacf(get(paste0("hydro_", k, "_", i, "_", j))$relhum, plot = FALSE)$acf))
      assign(paste0("df_pacf_relhum_", k, "_", i),
             bind_cols(get(paste0("df_pacf_relhum_", k, "_", i)), get(paste0("pacf_relhum_", k, "_", i, "_", j))))
    }
    assign(paste0("pacf_mean_relhum_", k, "_", i), rowMeans(get(paste0("df_pacf_relhum_", k, "_", i))[, 2:11]))
  }
}
# 10304, Winter: Autokorrelationskoeffizient bis lag = 3 (durchgängig) signifikant
# 10304, Summer: Autokorrelationskoeffizient bis lag = 3 (durchgängig) signifikant
# 11502, Winter: Autokorrelationskoeffizient bis lag = 3 (durchgängig) signifikant
# 11502, Summer: Autokorrelationskoeffizient bis lag = 3 (durchgängig) signifikant
# 20203, Winter: Autokorrelationskoeffizient bis lag = 5 (durchgängig) signifikant
# 20203, Summer: Autokorrelationskoeffizient bis lag = 4 (nicht durchgängig), bis lag = 1 (durchgängig) signifikant

# PACF-Plots for relhum (across all members) seperated by each waterlevel and season
plot(x = 1:37, y = pacf_mean_relhum_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der relativen Luftfeuchte für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_relhum_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der relativen Luftfeuchte für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_relhum_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der relativen Luftfeuchte für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_relhum_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der relativen Luftfeuchte für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_relhum_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der relativen Luftfeuchte für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_relhum_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der relativen Luftfeuchte für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

## Soilwater
# autocorrelation function, mean acf for soilwater (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_acf_soilwater_", k, "_", i), data.frame(Names = 1:38))
    for (j in levels(hydro_total$member)) {
      assign(paste0("acf_soilwater_", k, "_", i, "_", j),
             as.data.frame(acf(get(paste0("hydro_", k, "_", i, "_", j))$soilwater, plot = FALSE)$acf))
      assign(paste0("df_acf_soilwater_", k, "_", i),
             bind_cols(get(paste0("df_acf_soilwater_", k, "_", i)), get(paste0("acf_soilwater_", k, "_", i, "_", j))))
    }
    assign(paste0("acf_mean_soilwater_", k, "_", i), rowMeans(get(paste0("df_acf_soilwater_", k, "_", i))[, 2:11]))
  }
}
# -> Soilwater ist in allen Teildatensätzen positiv autokorreliert -> Modell, das diese Autokorrelation berücksichtigt, nötig

# ACF-Plots for soilwater (across all members) separated by each waterlevel and season
plot(x = 0:37, y = acf_mean_soilwater_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_soilwater_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_soilwater_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_soilwater_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_soilwater_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_soilwater_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

# partial autocorrelation function, mean pacf for soilwater (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_pacf_soilwater_", k, "_", i), data.frame(Names = 1:37))
    for (j in levels(hydro_total$member)) {
      assign(paste0("pacf_soilwater_", k, "_", i, "_", j),
             as.data.frame(pacf(get(paste0("hydro_", k, "_", i, "_", j))$soilwater, plot = FALSE)$acf))
      assign(paste0("df_pacf_soilwater_", k, "_", i),
             bind_cols(get(paste0("df_pacf_soilwater_", k, "_", i)), get(paste0("pacf_soilwater_", k, "_", i, "_", j))))
    }
    assign(paste0("pacf_mean_soilwater_", k, "_", i), rowMeans(get(paste0("df_pacf_soilwater_", k, "_", i))[, 2:11]))
  }
}
# 10304, Winter: Autokorrelationskoeffizient bis lag = 2 (durchgängig) signifikant
# 10304, Summer: Autokorrelationskoeffizient bis lag = 3 (durchgängig) signifikant
# 11502, Winter: Autokorrelationskoeffizient bis lag = 3 (durchgängig) signifikant
# 11502, Summer: Autokorrelationskoeffizient bis lag = 3 (durchgängig) signifikant
# 20203, Winter: Autokorrelationskoeffizient bis lag = 2 (durchgängig) signifikant
# 20203, Summer: Autokorrelationskoeffizient bis lag = 2 (durchgängig) signifikant

# PACF-Plots for soilwater (across all members) seperated by each waterlevel and season
plot(x = 1:37, y = pacf_mean_soilwater_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_soilwater_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_soilwater_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_soilwater_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_soilwater_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_soilwater_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der oberflächennahen Bodenfeuchte für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")


## Snowstorage
# autocorrelation function, mean acf for snowstorage (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_acf_snowstorage_", k, "_", i), data.frame(Names = 1:38))
    for (j in levels(hydro_total$member)) {
      assign(paste0("acf_snowstorage_", k, "_", i, "_", j),
             as.data.frame(acf(get(paste0("hydro_", k, "_", i, "_", j))$snowstorage, plot = FALSE)$acf))
      assign(paste0("df_acf_snowstorage_", k, "_", i),
             bind_cols(get(paste0("df_acf_snowstorage_", k, "_", i)), get(paste0("acf_snowstorage_", k, "_", i, "_", j))))
    }
    assign(paste0("acf_mean_snowstorage_", k, "_", i), rowMeans(get(paste0("df_acf_snowstorage_", k, "_", i))[, 2:11]))
  }
}
# -> Snowstorage ist in allen Teildatensätzen positiv autokorreliert -> Modell, das diese Autokorrelation berücksichtigt, nötig

# ACF-Plots for snowstorage (across all members) separated by each waterlevel and season
plot(x = 0:37, y = acf_mean_snowstorage_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Schneespeichers für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_snowstorage_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Schneespeichers für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_snowstorage_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Schneespeichers für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_snowstorage_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Schneespeichers für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_snowstorage_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Schneespeichers für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_snowstorage_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Schneespeichers für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

# partial autocorrelation function, mean pacf for snowstorage (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_pacf_snowstorage_", k, "_", i), data.frame(Names = 1:37))
    for (j in levels(hydro_total$member)) {
      assign(paste0("pacf_snowstorage_", k, "_", i, "_", j),
             as.data.frame(pacf(get(paste0("hydro_", k, "_", i, "_", j))$snowstorage, plot = FALSE)$acf))
      assign(paste0("df_pacf_snowstorage_", k, "_", i),
             bind_cols(get(paste0("df_pacf_snowstorage_", k, "_", i)), get(paste0("pacf_snowstorage_", k, "_", i, "_", j))))
    }
    assign(paste0("pacf_mean_snowstorage_", k, "_", i), rowMeans(get(paste0("df_pacf_snowstorage_", k, "_", i))[, 2:11]))
  }
}
# 10304, Winter: Autokorrelationskoeffizient bis lag = 3 (durchgängig) signifikant
# 10304, Summer: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 11502, Winter: Autokorrelationskoeffizient bis lag = 4 (durchgängig) signifikant
# 11502, Summer: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 20203, Winter: Autokorrelationskoeffizient bis lag = 3 (durchgängig) signifikant
# 20203, Summer: Autokorrelationskoeffizient bis lag = 2 (durchgängig) signifikant

# PACF-Plots for snowstorage (across all members) seperated by each waterlevel and season
plot(x = 1:37, y = pacf_mean_snowstorage_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Schneespeichers für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_snowstorage_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Schneespeichers für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_snowstorage_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Schneespeichers für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_snowstorage_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Schneespeichers für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_snowstorage_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Schneespeichers für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_snowstorage_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Schneespeichers für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")


## Groundwaterdepth
# autocorrelation function, mean acf for groundwaterdepth (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_acf_groundwaterdepth_", k, "_", i), data.frame(Names = 1:38))
    for (j in levels(hydro_total$member)) {
      assign(paste0("acf_groundwaterdepth_", k, "_", i, "_", j),
             as.data.frame(acf(get(paste0("hydro_", k, "_", i, "_", j))$groundwaterdepth, plot = FALSE)$acf))
      assign(paste0("df_acf_groundwaterdepth_", k, "_", i),
             bind_cols(get(paste0("df_acf_groundwaterdepth_", k, "_", i)), get(paste0("acf_groundwaterdepth_", k, "_", i, "_", j))))
    }
    assign(paste0("acf_mean_groundwaterdepth_", k, "_", i), rowMeans(get(paste0("df_acf_groundwaterdepth_", k, "_", i))[, 2:11]))
  }
}
# -> Groundwaterdepth ist in allen Teildatensätzen positiv autokorreliert -> Modell, das diese Autokorrelation berücksichtigt, nötig

# ACF-Plots for groundwaterdepth (across all members) separated by each waterlevel and season
plot(x = 0:37, y = acf_mean_groundwaterdepth_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Grundwasserstands für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_groundwaterdepth_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Grundwasserstands für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_groundwaterdepth_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Grundwasserstands für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_groundwaterdepth_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Grundwasserstands für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_groundwaterdepth_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Grundwasserstands für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_groundwaterdepth_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. des Grundwasserstands für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

# partial autocorrelation function, mean pacf for groundwaterdepth (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_pacf_groundwaterdepth_", k, "_", i), data.frame(Names = 1:37))
    for (j in levels(hydro_total$member)) {
      assign(paste0("pacf_groundwaterdepth_", k, "_", i, "_", j),
             as.data.frame(pacf(get(paste0("hydro_", k, "_", i, "_", j))$groundwaterdepth, plot = FALSE)$acf))
      assign(paste0("df_pacf_groundwaterdepth_", k, "_", i),
             bind_cols(get(paste0("df_pacf_groundwaterdepth_", k, "_", i)), get(paste0("pacf_groundwaterdepth_", k, "_", i, "_", j))))
    }
    assign(paste0("pacf_mean_groundwaterdepth_", k, "_", i), rowMeans(get(paste0("df_pacf_groundwaterdepth_", k, "_", i))[, 2:11]))
  }
}
# 10304, Winter: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 10304, Summer: Autokorrelationskoeffizient bis lag = 2 (durchgängig) signifikant
# 11502, Winter: Autokorrelationskoeffizient bis lag = 3 (durchgängig) signifikant
# 11502, Summer: Autokorrelationskoeffizient bis lag = 5 (durchgängig) signifikant
# 20203, Winter: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 20203, Summer: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant

# PACF-Plots for groundwaterdepth (across all members) seperated by each waterlevel and season
plot(x = 1:37, y = pacf_mean_groundwaterdepth_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Grundwasserstands für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_groundwaterdepth_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Grundwasserstands für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_groundwaterdepth_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Grundwasserstands für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_groundwaterdepth_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Grundwasserstands für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_groundwaterdepth_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Grundwasserstands für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_groundwaterdepth_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. des Grundwasserstands für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")


## Infiltration
# autocorrelation function, mean acf for infiltration (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_acf_infiltration_", k, "_", i), data.frame(Names = 1:38))
    for (j in levels(hydro_total$member)) {
      assign(paste0("acf_infiltration_", k, "_", i, "_", j),
             as.data.frame(acf(get(paste0("hydro_", k, "_", i, "_", j))$infiltration, plot = FALSE)$acf))
      assign(paste0("df_acf_infiltration_", k, "_", i),
             bind_cols(get(paste0("df_acf_infiltration_", k, "_", i)), get(paste0("acf_infiltration_", k, "_", i, "_", j))))
    }
    assign(paste0("acf_mean_infiltration_", k, "_", i), rowMeans(get(paste0("df_acf_infiltration_", k, "_", i))[, 2:11]))
  }
}
# -> Infiltration ist in allen Teildatensätzen positiv autokorreliert -> Modell, das diese Autokorrelation berücksichtigt, nötig

# ACF-Plots for infiltration (across all members) separated by each waterlevel and season
plot(x = 0:37, y = acf_mean_infiltration_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der Wasserleitfähigkeit für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_infiltration_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der Wasserleitfähigkeit für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_infiltration_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der Wasserleitfähigkeit für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_infiltration_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der Wasserleitfähigkeit für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_infiltration_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der Wasserleitfähigkeit für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 0:37, y = acf_mean_infiltration_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "ACF", main = "Autokor. zw. Beob. der Wasserleitfähigkeit für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

# partial autocorrelation function, mean pacf for infiltration (across all members) separated by each waterlevel and season
for (k in c("winter", "summer")) {
  for (i in pegel_ids) {
    assign(paste0("df_pacf_infiltration_", k, "_", i), data.frame(Names = 1:37))
    for (j in levels(hydro_total$member)) {
      assign(paste0("pacf_infiltration_", k, "_", i, "_", j),
             as.data.frame(pacf(get(paste0("hydro_", k, "_", i, "_", j))$infiltration, plot = FALSE)$acf))
      assign(paste0("df_pacf_infiltration_", k, "_", i),
             bind_cols(get(paste0("df_pacf_infiltration_", k, "_", i)), get(paste0("pacf_infiltration_", k, "_", i, "_", j))))
    }
    assign(paste0("pacf_mean_infiltration_", k, "_", i), rowMeans(get(paste0("df_pacf_infiltration_", k, "_", i))[, 2:11]))
  }
}
# 10304, Winter: Autokorrelationskoeffizient bis lag = 6 (durchgängig) signifikant
# 10304, Summer: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 11502, Winter: Autokorrelationskoeffizient bis lag = 3 (nicht durchgängig), bis lag = 1 (durchgängig) signifikant
# 11502, Summer: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant
# 20203, Winter: Autokorrelationskoeffizient bis lag = 2 (durchgängig) signifikant
# 20203, Summer: Autokorrelationskoeffizient bis lag = 1 (durchgängig) signifikant

# PACF-Plots for infiltration (across all members) seperated by each waterlevel and season
plot(x = 1:37, y = pacf_mean_infiltration_winter_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Wasserleitfähigkeit für Isar Mittenwald, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_infiltration_summer_10304, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Wasserleitfähigkeit für Isar Mittenwald, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_infiltration_winter_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Wasserleitfähigkeit für Iller Kempten, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_infiltration_summer_11502, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Wasserleitfähigkeit für Iller Kempten, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_infiltration_winter_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Wasserleitfähigkeit für fränk. Saale Salz, Winter") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")
plot(x = 1:37, y = pacf_mean_infiltration_summer_20203, ylim = c(-0.1, 1), type = "h", xlab = "lag", ylab = "PACF", main = "Part. Autokor. zw. Beob. der Wasserleitfähigkeit für fränk. Saale Salz, Sommer") + abline(h = 0) + abline(h = c(-0.05, 0.05), col = "blue", lty = "dashed")

```


```{r}
rolCorPlot(ts_hydro_summer_10304_kbe)
```


## Linear correlation with drainage
```{r}
hydro_winter_10304 <- readRDS(file = "./data/hydro_winter_10304.RDS")
hydro_winter_11502 <- readRDS(file = "./data/hydro_winter_11502.RDS")
hydro_winter_20203 <- readRDS(file = "./data/hydro_winter_20203.RDS")
hydro_summer_10304 <- readRDS(file = "./data/hydro_summer_10304.RDS")
hydro_summer_11502 <- readRDS(file = "./data/hydro_summer_11502.RDS")
hydro_summer_20203 <- readRDS(file = "./data/hydro_summer_20203.RDS")
  
plot_cor_drainage <- function(data, var, ...) {
  ggplot(data, aes(x = get(var), y = drainage)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = "lm", col = "red") +
  labs(...)
}

# Airtmp
cor_airtmp <- numeric(length = 6)
plot_cor_drainage(data = hydro_winter_10304, var = "airtmp", x = "airtmp", title = "Winter 10304 (Isar Mittenwald) correlation airtmp and drainage")
cor_airtmp[1] <- cor(hydro_winter_10304$drainage, hydro_winter_10304$airtmp)
plot_cor_drainage(data = hydro_summer_10304, var = "airtmp", x = "airtmp", title = "Summer 10304 (Isar Mittenwald) correlation airtmp and drainage")
cor_airtmp[2] <- cor(hydro_summer_10304$drainage, hydro_summer_10304$airtmp)

plot_cor_drainage(data = hydro_winter_11502, var = "airtmp", x = "airtmp", title = "Winter 11502 (Iller Kempten) correlation airtmp and drainage")
cor_airtmp[3] <- cor(hydro_winter_11502$drainage, hydro_winter_11502$airtmp)
plot_cor_drainage(data = hydro_summer_11502, var = "airtmp", x = "airtmp", title = "Summer 11502 (Iller Kempten) correlation airtmp and drainage")
cor_airtmp[4] <- cor(hydro_summer_11502$drainage, hydro_summer_11502$airtmp)

plot_cor_drainage(data = hydro_winter_20203, var = "airtmp", x = "airtmp", title = "Winter 20203 (Fränkische Saale Salz) correlation airtmp and drainage")
cor_airtmp[5] <- cor(hydro_winter_20203$drainage, hydro_winter_20203$airtmp)
plot_cor_drainage(data = hydro_summer_20203, var = "airtmp", x = "airtmp", title = "Summer 20203 (Fränkische Saale Salz) correlation airtmp and drainage")
cor_airtmp[6] <- cor(hydro_summer_20203$drainage, hydro_summer_20203$airtmp)

var(cor_airtmp)

# Soilwater
cor_soilwater <- numeric(length = 6)
plot_cor_drainage(data = hydro_winter_10304, var = "soilwater", x = "soilwater", title = "Winter 10304 (Isar Mittenwald) correlation soilwater and drainage")
cor_soilwater[1] <- cor(hydro_winter_10304$drainage, hydro_winter_10304$soilwater)
plot_cor_drainage(data = hydro_summer_10304, var = "soilwater", x = "soilwater", title = "Summer 10304 (Isar Mittenwald) correlation soilwater and drainage")
cor_soilwater[2] <- cor(hydro_summer_10304$drainage, hydro_summer_10304$soilwater)

plot_cor_drainage(data = hydro_winter_11502, var = "soilwater", x = "soilwater", title = "Winter 11502 (Iller Kempten) correlation soilwater and drainage")
cor_soilwater[3] <- cor(hydro_winter_11502$drainage, hydro_winter_11502$soilwater)
plot_cor_drainage(data = hydro_summer_11502, var = "soilwater", x = "soilwater", title = "Summer 11502 (Iller Kempten) correlation soilwater and drainage")
cor_soilwater[4] <- cor(hydro_summer_11502$drainage, hydro_summer_11502$soilwater)

plot_cor_drainage(data = hydro_winter_20203, var = "soilwater", x = "soilwater", title = "Winter 20203 (Fränkische Saale Salz) correlation soilwater and drainage")
cor_soilwater[5] <- cor(hydro_winter_20203$drainage, hydro_winter_20203$soilwater)
plot_cor_drainage(data = hydro_summer_20203, var = "soilwater", x = "soilwater", title = "Summer 20203 (Fränkische Saale Salz) correlation soilwater and drainage")
cor_soilwater[6] <- cor(hydro_summer_20203$drainage, hydro_summer_20203$soilwater)

var(cor_soilwater)


# Precip
cor_precip <- numeric(length = 6)
plot_cor_drainage(data = hydro_winter_10304, var = "precip", x = "precip", title = "Winter 10304 (Isar Mittenwald) correlation precip and drainage")
cor_precip[1] <- cor(hydro_winter_10304$drainage, hydro_winter_10304$precip)
plot_cor_drainage(data = hydro_summer_10304, var = "precip", x = "precip", title = "Summer 10304 (Isar Mittenwald) correlation precip and drainage")
cor_precip[2] <- cor(hydro_summer_10304$drainage, hydro_summer_10304$precip)

plot_cor_drainage(data = hydro_winter_11502, var = "precip", x = "precip", title = "Winter 11502 (Iller Kempten) correlation precip and drainage")
cor_precip[3] <- cor(hydro_winter_11502$drainage, hydro_winter_11502$precip)
plot_cor_drainage(data = hydro_summer_11502, var = "precip", x = "precip", title = "Summer 11502 (Iller Kempten) correlation precip and drainage")
cor_precip[4] <- cor(hydro_summer_11502$drainage, hydro_summer_11502$precip)

plot_cor_drainage(data = hydro_winter_20203, var = "precip", x = "precip", title = "Winter 20203 (Fränkische Saale Salz) correlation precip and drainage")
cor_precip[5] <- cor(hydro_winter_20203$drainage, hydro_winter_20203$precip)
plot_cor_drainage(data = hydro_summer_20203, var = "precip", x = "precip", title = "Summer 20203 (Fränkische Saale Salz) correlation precip and drainage")
cor_precip[6] <- cor(hydro_summer_20203$drainage, hydro_summer_20203$precip)

var(cor_precip)
```

# correlation matrices between drivers
```{r}
hydro_winter_10304_kbe <- readRDS(file = "./data/hydro_winter_10304_kbe.RDS")
hydro_winter_11502_kbe <- readRDS(file = "./data/hydro_winter_11502_kbe.RDS")
hydro_winter_20203_kbe <- readRDS(file = "./data/hydro_winter_20203_kbe.RDS")
hydro_summer_10304_kbe <- readRDS(file = "./data/hydro_summer_10304_kbe.RDS")
hydro_summer_11502_kbe <- readRDS(file = "./data/hydro_summer_11502_kbe.RDS")
hydro_summer_20203_kbe <- readRDS(file = "./data/hydro_summer_20203_kbe.RDS")

subset_hydro_winter_10304_kbe <- subset(hydro_winter_10304_kbe, select = c(avg_precip, avg_airtmp, avg_glorad, avg_relhum, avg_soilwater, avg_snowstorage, groundwaterdepth, avg_infiltration, max_precip, avg_snowstorage_drain))
colnames(subset_hydro_winter_10304_kbe) <- c("Mittlerer Niederschlag", "Mittlere Lufttemperatur", "Mittlere Strahlung", "Mittlere relative Luftfeuchte", "Mittlere Bodenfeuchte", "Mittlerer Schneespeicher", "Grundwasserstand", "Mittlere Versickerung", "Maximaler Niederschlag", "Mittlere Schneeschmelze")

subset_hydro_winter_11502_kbe <- subset(hydro_winter_11502_kbe, select = c(avg_precip, avg_airtmp, avg_glorad, avg_relhum, avg_soilwater, avg_snowstorage, groundwaterdepth, avg_infiltration, max_precip, avg_snowstorage_drain))
colnames(subset_hydro_winter_11502_kbe) <- c("Mittlerer Niederschlag", "Mittlere Lufttemperatur", "Mittlere Strahlung", "Mittlere relative Luftfeuchte", "Mittlere Bodenfeuchte", "Mittlerer Schneespeicher", "Grundwasserstand", "Mittlere Versickerung", "Maximaler Niederschlag", "Mittlere Schneeschmelze")

subset_hydro_winter_20203_kbe <- subset(hydro_winter_20203_kbe, select = c(avg_precip, avg_airtmp, avg_glorad, avg_relhum, avg_soilwater, avg_snowstorage, groundwaterdepth, avg_infiltration, max_precip, avg_snowstorage_drain))
colnames(subset_hydro_winter_20203_kbe) <- c("Mittlerer Niederschlag", "Mittlere Lufttemperatur", "Mittlere Strahlung", "Mittlere relative Luftfeuchte", "Mittlere Bodenfeuchte", "Mittlerer Schneespeicher", "Grundwasserstand", "Mittlere Versickerung", "Maximaler Niederschlag", "Mittlere Schneeschmelze")

subset_hydro_summer_10304_kbe <- subset(hydro_summer_10304_kbe, select = c(avg_precip, avg_airtmp, avg_glorad, avg_relhum, avg_soilwater, avg_snowstorage, groundwaterdepth, avg_infiltration, max_precip, avg_snowstorage_drain))
colnames(subset_hydro_summer_10304_kbe) <- c("Mittlerer Niederschlag", "Mittlere Lufttemperatur", "Mittlere Strahlung", "Mittlere relative Luftfeuchte", "Mittlere Bodenfeuchte", "Mittlerer Schneespeicher", "Grundwasserstand", "Mittlere Versickerung", "Maximaler Niederschlag", "Mittlere Schneeschmelze")

subset_hydro_summer_11502_kbe <- subset(hydro_summer_11502_kbe, select = c(avg_precip, avg_airtmp, avg_glorad, avg_relhum, avg_soilwater, avg_snowstorage, groundwaterdepth, avg_infiltration, max_precip, avg_snowstorage_drain))
colnames(subset_hydro_summer_11502_kbe) <- c("Mittlerer Niederschlag", "Mittlere Lufttemperatur", "Mittlere Strahlung", "Mittlere relative Luftfeuchte", "Mittlere Bodenfeuchte", "Mittlerer Schneespeicher", "Grundwasserstand", "Mittlere Versickerung", "Maximaler Niederschlag", "Mittlere Schneeschmelze")

subset_hydro_summer_20203_kbe <- subset(hydro_summer_20203_kbe, select = c(avg_precip, avg_airtmp, avg_glorad, avg_relhum, avg_soilwater, avg_snowstorage, groundwaterdepth, avg_infiltration, max_precip, avg_snowstorage_drain))
colnames(subset_hydro_summer_20203_kbe) <- c("Mittlerer Niederschlag", "Mittlere Lufttemperatur", "Mittlere Strahlung", "Mittlere relative Luftfeuchte", "Mittlere Bodenfeuchte", "Mittlerer Schneespeicher", "Grundwasserstand", "Mittlere Versickerung", "Maximaler Niederschlag", "Mittlere Schneeschmelze")

saveRDS(subset_hydro_winter_10304, "./data/corrplots/subset_hydro_winter_10304.RDS")
saveRDS(subset_hydro_winter_11502, "./data/corrplots/subset_hydro_winter_11502.RDS")
saveRDS(subset_hydro_winter_20203, "./data/corrplots/subset_hydro_winter_20203.RDS")

saveRDS(subset_hydro_summer_10304, "./data/corrplots/subset_hydro_summer_10304.RDS")
saveRDS(subset_hydro_summer_11502, "./data/corrplots/subset_hydro_summer_11502.RDS")
saveRDS(subset_hydro_summer_20203, "./data/corrplots/subset_hydro_summer_20203.RDS")



library(corrplot)

### pearson
## winter

# 10304
correlationm_pearson_winter_10304_kbe <- cor(subset_hydro_winter_10304_kbe, method = "pearson", use = "complete.obs")
corrplot(correlationm_pearson_winter_10304_kbe, method = "number", type = "lower", tl.col = "black")


# 11502
correlationm_pearson_winter_11502_kbe <- cor(subset_hydro_winter_11502_kbe, method = "pearson", use = "complete.obs")
corrplot(correlationm_pearson_winter_11502_kbe, method = "number", type = "lower", tl.col = "black")

# 20203
correlationm_pearson_winter_20203_kbe <- cor(subset_hydro_winter_20203_kbe, method = "pearson", use = "complete.obs")
corrplot(correlationm_pearson_winter_20203_kbe, method = "number", type = "lower", tl.col = "black")


## summer

# 10304
correlationm_pearson_summer_10304_kbe <- cor(subset_hydro_summer_10304_kbe, method = "pearson", use = "complete.obs")
corrplot(correlationm_pearson_summer_10304_kbe, method = "number", type = "lower", tl.col = "black")

# 11502
correlationm_pearson_summer_11502_kbe <- cor(subset_hydro_summer_11502_kbe, method = "pearson", use = "complete.obs")
corrplot(correlationm_pearson_summer_11502_kbe, method = "number", type = "lower", tl.col = "black")

# 20203
correlationm_pearson_summer_20203_kbe <- cor(subset_hydro_summer_20203_kbe, method = "pearson", use = "complete.obs")
corrplot(correlationm_pearson_summer_20203_kbe, method = "number", type = "lower", tl.col = "black")




### spearman
## winter

# 10304
correlationm_spearman_winter_10304_kbe <- cor(subset_hydro_winter_10304_kbe, method = "spearman", use = "complete.obs")
corrplot(correlationm_spearman_winter_10304_kbe, method = "number", type = "lower", tl.col = "black")

# 11502
correlationm_spearman_winter_11502_kbe <- cor(subset_hydro_winter_11502_kbe, method = "spearman", use = "complete.obs")
corrplot(correlationm_spearman_winter_11502_kbe, method = "number", type = "lower", tl.col = "black")

# 20203
correlationm_spearman_winter_20203_kbe <- cor(subset_hydro_winter_20203_kbe, method = "spearman", use = "complete.obs")
corrplot(correlationm_spearman_winter_20203_kbe, method = "number", type = "lower", tl.col = "black")


## summer

# 10304
correlationm_spearman_summer_10304_kbe <- cor(subset_hydro_summer_10304_kbe, method = "spearman", use = "complete.obs")
corrplot(correlationm_spearman_summer_10304_kbe, method = "number", type = "lower", tl.col = "black")

# 11502
correlationm_spearman_summer_11502_kbe <- cor(subset_hydro_summer_11502_kbe, method = "spearman", use = "complete.obs")
corrplot(correlationm_spearman_summer_11502_kbe, method = "number", type = "lower", tl.col = "black")

# 20203
correlationm_spearman_summer_20203_kbe <- cor(subset_hydro_summer_20203_kbe, method = "spearman", use = "complete.obs")
corrplot(correlationm_spearman_summer_20203_kbe, method = "number", type = "lower", tl.col = "black")

# allgemein:
hydro_summer_kbe <- subset(hydro_summer, member = "kbe")
hydro_winter_kbe <- subset(hydro_winter, member = "kbe")
hydro_total_kbe <- rbind(hydro_summer_kbe, hydro_winter_kbe)
subset_hydro_total_kbe <- subset(hydro_total_kbe, select = c(avg_precip, avg_airtmp, avg_glorad, avg_relhum, avg_soilwater, avg_snowstorage, groundwaterdepth, avg_infiltration, max_precip, avg_snowstorage_drain))
colnames(subset_hydro_total_kbe) <- c("Mittlerer Niederschlag", "Mittlere Lufttemperatur", "Mittlere Strahlung", "Mittlere relative Luftfeuchte", "Mittlere Bodenfeuchte", "Mittlerer Schneespeicher", "Grundwasserstand", "Mittlere Versickerung", "Maximaler Niederschlag", "Mittlere Schneeschmelze")

correlationm_spearman_total_kbe <- cor(subset_hydro_total_kbe, method = "spearman", use = "complete.obs")
corrplot(correlationm_spearman_total_kbe, method = "number", type = "lower", tl.col = "black")
```



# IQA bestimmen
```{r}
avg_soilwater <- c(hydro_winter_20203$avg_soilwater, hydro_winter_10304$avg_soilwater, hydro_winter_11502$avg_soilwater, hydro_summer_10304$avg_soilwater, hydro_summer_11502$avg_soilwater, hydro_summer_20203$avg_soilwater)
quantile(na.omit(avg_soilwater),  c(0.25, 0.75))
# IQA ca. 1.8%, Faktor 10 führt zu 18 0.1% Schritten


avg_snowstorage <- c(hydro_winter_20203$avg_snowstorage, hydro_winter_10304$avg_snowstorage, hydro_winter_11502$avg_snowstorage, hydro_summer_10304$avg_snowstorage, hydro_summer_11502$avg_snowstorage, hydro_summer_20203$avg_snowstorage)
quantile(na.omit(avg_snowstorage),c(0.25, 0.75))
# IQA ca. 36, also führt Faktor 0.5 zu 18 2mm Schritten


avg_airtmp <- c(hydro_winter_20203$avg_airtmp, hydro_winter_10304$avg_airtmp, hydro_winter_11502$avg_airtmp, hydro_summer_10304$avg_airtmp, hydro_summer_11502$avg_airtmp, hydro_summer_20203$avg_airtmp)
quantile(na.omit(avg_airtmp),  c(0.25, 0.75))
# IQA ca. 5.3, also führt Faktor 3.3 zu 16 0.3 grad Schritten


avg_precip <- c(hydro_winter_20203$avg_precip, hydro_winter_10304$avg_precip, hydro_winter_11502$avg_precip, hydro_summer_10304$avg_precip, hydro_summer_11502$avg_precip, hydro_summer_20203$avg_precip)
quantile(na.omit(avg_precip),  c(0.25, 0.75))
# IQA ca. 3.3, also führt Faktor 5 zu 16 0.2 mm/m^2  Schritten



avg_glorad <- c(hydro_winter_20203$avg_glorad, hydro_winter_10304$avg_glorad, hydro_winter_11502$avg_glorad, hydro_summer_10304$avg_glorad, hydro_summer_11502$avg_glorad, hydro_summer_20203$avg_glorad)
quantile(na.omit(avg_glorad),  c(0.25, 0.75))
# IQA ca. 240, also führt Faktor 0.1 zu 24 10 Wh/m^2 Schritten


avg_relhum <- c(hydro_winter_20203$avg_relhum, hydro_winter_10304$avg_relhum, hydro_winter_11502$avg_relhum, hydro_summer_10304$avg_relhum, hydro_summer_11502$avg_relhum, hydro_summer_20203$avg_relhum)
quantile(na.omit(avg_relhum),  c(0.25, 0.75))
# IQA ca. 11, weshalb wir Faktor 1 nutzen


avg_infiltration <- c(hydro_winter_20203$avg_infiltration, hydro_winter_10304$avg_infiltration, hydro_winter_11502$avg_infiltration, hydro_summer_10304$avg_infiltration, hydro_summer_11502$avg_infiltration, hydro_summer_20203$avg_infiltration)
quantile(na.omit(avg_infiltration),  c(0.25, 0.75))
# IQA ca. 0.35, also führt Faktor 50 zu 17 0.02 (Infiltration Einheit)Schritten


groundwaterdepth <- c(hydro_summer_10304$groundwaterdepth, hydro_winter_10304$groundwaterdepth)
quantile(na.omit(groundwaterdepth),  c(0.25, 0.75))
# IQA ca. 0.45, also führt Faktor 33 zu 17 1/30 m schritten


avg_snowstorage_drain_scaled <- c(hydro_summer_10304$avg_snowstorage_drain_scaled, hydro_winter_10304$avg_snowstorage_drain_scaled)
quantile(na.omit(avg_snowstorage_drain_scaled),  c(0.25, 0.75))
# IQA ca. 55, also führt Faktor 0.33333 zu 17 3er Schritten



max_precip <- c(hydro_summer$max_precip, hydro_winter$max_precip_winter)
quantile(na.omit(max_precip),  c(0.25, 0.75))
# IQA ca. 15, also führt Faktor 1 zu 15 1ner Schritten
```