05_Annual_data_all.Rmd

---
title: "05 Collect annual data of all types"
output: 
  html_document:
    keep_md: true
    toc: true
    toc_depth: 3
    toc_float: true
    code_folding: hide
    df_print: paged   
---

__Test file__ for producing summarized data by year and quarter   
Actual data used (in folder 'Data_produced') are made in __script 05b__    


## 0. Libraries
```{r}
library(tidyverse)
library(readxl)
library(broom)
library(lubridate)
# library(pander)
```

## 1. Folders and files    
Content of folders (see code)  
```{r}
dir("Datasett")
dir("Datasett/River data (from OKA)")
dir("Datasett/hydrografi")
dir("Datasett/Bløtbunn")
dir("Datasett/hardbunn_kopi")
dir("Datasett/Plankton")
```

## 2. River data   
Content of folders (see code)
```{r}
dir("Datasett/River data (from OKA)/Annual mean flows")
dir("Datasett/River data (from OKA)/Monthly loads")
dir("Datasett/River data (from OKA)/Concentrations (individual samples)")
dir("Datasett/River data (from OKA)/Monthly flow-weighted concentrations")
```

### a. Data of monthly river loads and total discharge  
```{r}
df1 <- read_excel("Datasett/River data (from OKA)/Monthly loads/Storelva_monthly loads.xlsx")  
df2 <- read_excel("Datasett/River data (from OKA)/Monthly loads/Gjerstadelva_Nidelva_monthly loads.xlsx")  
df3 <- read_excel("Datasett/River data (from OKA)/Monthly loads/RIDx5_monthly loads.xlsx")  

# head(df1, 3)
# head(df2, 3)
# head(df3, 3)

# colnames(df1) %>% dput()
# colnames(df2) %>% dput()
# colnames(df3) %>% dput()   RID rivers also include PO4, SiO2 and SPM

df_rivers <- bind_rows(df1[-1,], df2[-1,], df3[-1,])

# colnames(df_rivers) %>% dput()
vars <- c("TrspTot TOTN", "TrspTot NO3-N", "TrspTot NH4-N", "TrspTot TOTP",           
          "TrspTot TOC", "TrspTot ALK", "TrspTot Ca", "DisTot")
for (var in vars)
  df_rivers[[var]] <- as.numeric(df_rivers[[var]])
df_rivers$Time <- with(df_rivers, lubridate::ymd(paste(Year, Month, "15")))

# Add "_" in column names (TrspTot Ca -> TrspTot_Ca)
colnames(df_rivers) <- sub(" ", "_", colnames(df_rivers), fixed = TRUE)

# endre rekkefølge på elver fra nord til sør for ggplot
df_rivers$Station_name <- factor(df_rivers$Station_name, levels = c("Glomma ved Sarpsfoss", "Drammenselva", "Numedalslågen", "Skienselva", "Søndeledelva v. Søndeleddammen", "Storelva v/ Nes verk", "Nidelva ovenf. Rygene", "Otra"))

# endre navn på elver
sel <- levels(df_rivers$Station_name) == "Glomma ved Sarpsfoss"; sum(sel)
levels(df_rivers$Station_name)[sel] <- "Glomma"

sel <- levels(df_rivers$Station_name) == "Søndeledelva v. Søndeleddammen"; sum(sel)
levels(df_rivers$Station_name)[sel] <- "Gjerstadelva"

sel <- levels(df_rivers$Station_name) == "Storelva v/ Nes verk"; sum(sel)
levels(df_rivers$Station_name)[sel] <- "Storelva"

sel <- levels(df_rivers$Station_name) == "Nidelva ovenf. Rygene"; sum(sel)
levels(df_rivers$Station_name)[sel] <- "Nidelva"

levels(df_rivers$Station_name)

# Dropp Otra fra plot og analyser (nedstrøms)
df_rivers <- df_rivers %>% 
  filter(Station_name != "Otra") %>% 
  droplevels()

#print (df_rivers)

# Table of available data for each river
tb <- df_rivers %>% 
  gather("Variable", Value, TrspTot_TOTN:DisTot) %>%
  filter(!is.na(Value)) %>%
  xtabs(~Station_name + Variable, .)
tb

```

### b. Local rivers, plot monthly mean discharge by station
```{r}
gg <- df_rivers %>%
  filter(substr(Station_name, 1, 4) %in% c("Nide","Gjer","Stor")) %>%
  group_by(Station_name, Month) %>%
  summarise(Mean = mean(DisTot, na.rm = TRUE), 
            Q10 = quantile(DisTot, 0.1, na.rm = TRUE), 
            Q90 = max(DisTot, 0.9, na.rm = TRUE)) %>%
  ggplot(., aes(Month, Mean)) + 
    geom_ribbon(aes(ymin = Q10, ymax = Q90), fill = "lightgreen") +
    geom_line() + geom_point() +
    facet_wrap(~Station_name)
gg
# gg + scale_y_log10()
```

### c. Distant rivers, plot monthly mean discharge by station  
Including Otra
```{r}
gg <- df_rivers %>%
  filter(!substr(Station_name, 1, 4) %in% c("Nide","Gjer","Stor")) %>%
  group_by(Station_name, Month) %>%
  summarise(Mean = mean(DisTot, na.rm = TRUE), 
            Q10 = quantile(DisTot, 0.1, na.rm = TRUE), 
            Q90 = max(DisTot, 0.9, na.rm = TRUE)) %>%
  ggplot(., aes(Month, Mean)) + 
    geom_ribbon(aes(ymin = Q10, ymax = Q90), fill = "lightgreen") +
    geom_line() + geom_point() +
    facet_wrap(~Station_name)
gg
# gg + scale_y_log10()
```

### d. Summarize by "local rivers"/"distant rivers" and quarter  
    * Exclude Otra
```{r}
df_rivers_summ <- df_rivers %>%
  filter(!Station_name %in% "Otra") %>%
  mutate(River_type = ifelse(substr(Station_name, 1, 4) %in% c("Nide","Gjer","Stor"), "Local", "Distant")) %>%
  mutate(Quarter = case_when(
    Month %in% 1:3 ~ 1,
    Month %in% 4:6 ~ 2,
    Month %in% 7:9 ~ 3,
    Month %in% 10:12 ~ 4
  )) %>%
  group_by(River_type, Year, Quarter) %>%
  summarise_at(c("TrspTot_TOTN", "TrspTot_NO3-N", "TrspTot_TOTP", "TrspTot_TOC", "DisTot"), mean, na.rm = TRUE)

```

### e. Plot in order to check data
```{r}
ggplot(df_rivers_summ, aes(Year, DisTot, color = River_type)) +
  geom_smooth() +
  geom_point() + 
  facet_grid(River_type~Quarter, scales = "free_y", labeller = label_both)

```

### f. Timing and size of spring flood
```{r}
df_rivers_springflood_allyears <- df_rivers %>%
  group_by(Station_name, Year) %>%
  mutate(DisTot_max = max(DisTot[Month %in% 1:6]), na.rm = TRUE) %>%
  group_by(Station_name) %>%
  summarize(DisTot_max_mean = mean(DisTot_max, na.rm = TRUE))
df_rivers_springflood_allyears

df_rivers_springflood <- df_rivers %>%
  filter(Month %in% 1:6) %>%
  group_by(Station_name, Year) %>%
  mutate(DisTot_max = max(DisTot), na.rm = TRUE) %>%
  group_by(Station_name) %>%
  mutate(DisTot_max_mean = mean(DisTot_max, na.rm = TRUE)) %>%
  ungroup() %>%
  group_by(Station_name, Year) %>%
  summarize(DisTot_max_rel = max(DisTot/DisTot_max_mean*100, na.rm = TRUE),
            DisTot_max_month = Month[DisTot == DisTot_max][1],
            DisTot_40perc = Month[DisTot >= 0.40*DisTot_max_mean][1],
            DisTot_60perc = Month[DisTot >= 0.60*DisTot_max_mean][1],
            DisTot_80perc = Month[DisTot >= 0.80*DisTot_max_mean][1])
```

### g. Save both
```{r}
# write.csv(df_rivers_summ, "Data_produced/05_df_rivers_summ.csv", row.names = FALSE, quote = FALSE)
# write.csv(df_rivers_springflood, "Data_produced/05_df_rivers_springflood.csv", 
#           row.names = FALSE, quote = FALSE)
```


### h1. Plot of max flood
```{r}
ggplot(df_rivers_springflood, aes(Year, DisTot_max_rel)) +
  geom_smooth() + geom_point() +
  facet_wrap(~Station_name)
```

### h2. Plot of flood timing
```{r}
df_rivers_springflood %>%
  gather("Parameter", "Month", DisTot_max_month, DisTot_40perc, DisTot_60perc, DisTot_80perc) %>%
  ggplot(aes(Year, Month, group = Parameter, color = Parameter)) +
  geom_smooth(method = "lm") + geom_point() +
  facet_wrap(~Station_name)
```

## 3. Hydrological data
### a. Read data
```{r}
load("Datasett/Hydrografi/Arendal_allvars_1990_2016.Rdata")
Df.Arendal$Month <- Df.Arendal$Dato %>% as.character() %>% substr(6,7) %>% as.numeric()
Df.Arendal$Year <- Df.Arendal$Dato %>% as.character() %>% substr(1,4) %>% as.numeric()
Df.Arendal$Time <- ymd_hms(paste(Df.Arendal$Dato, "00:00:00"))   # R's time format
```

### b. Summarize by depth bins and quarter    
    * Depth bins = 0-10, 10-30, 30-50
    * Quarters starting with March
```{r}
df_hydro_summ <- Df.Arendal %>%
  mutate(
    Quarter = case_when(
      Month %in% 1:2 ~ 1,
      Month %in% 3:5 ~ 2,
      Month %in% 6:8 ~ 3,
      Month %in% 9:11 ~ 4,
      Month %in% 12 ~ 1),
    Year2 = case_when(
      Month == 12 ~ Year + 1,
      Month < 12 ~ Year),
    Depth = case_when(
      Depth %in% c(0,5,10) ~ "Surface",
      Depth %in% c(20,30) ~ "Intermediate",
      Depth %in% c(50,75) ~ "Deep")
    ) %>%
  group_by(Year2, Quarter, Depth) %>%
  summarize_at(vars(Temperatur:Siktdyp), mean, na.rm = TRUE) %>%
  rename(Year = Year2)

df_hydro_summ$Depth <- factor(df_hydro_summ$Depth, 
                                 levels = c("Surface", "Intermediate", "Deep"))

head (df_hydro_summ)
print (df_hydro_summ)

str(df_hydro_summ)


# plot variable by quartes vr year
gg1 <- df_hydro_summ %>% 
  filter(Depth %in% "Surface") %>% 
  ggplot  (aes(Year, Temperatur)) +
  geom_smooth() + geom_point() +
   facet_wrap(~Quarter)

gg1


```

### c. Save
```{r}
# write.csv(df_hydro_summ, "Data_produced/05_df_hydro_summ.csv", row.names = FALSE, quote = FALSE)
```

### d. Plot temperature
```{r}
ggplot(df_hydro_summ, aes(Year, Temperatur, color = Depth)) +
  geom_smooth(method = "lm") + geom_point() +
  facet_grid(.~Quarter, labeller = label_both)
```

### e. Plot some nutrients etc.
```{r}
df_hydro_summ %>%
  gather("Var", "Concentration", NO2_NO3, TotP, TotN, TSM) %>%
  ggplot(aes(Year, Concentration, color = Depth)) +
    geom_smooth(method = "lm") + geom_point() +
    facet_grid(Var~Quarter, scales = "free_y", labeller = label_both)
```

## 4. Plankton  
Summarize main groups only  
    * Will add ordination scores (DCA) later
    
### a. Read plankton data
```{r}
df_plank <- read_excel("Datasett/Plankton/Planteplankton Arendal.xlsx") # range = "A1:V471"
df_plank$Year <- lubridate::year(df_plank$Dato)
df_plank$Month <- lubridate::month(df_plank$Dato)
```

### b. Plankton: Select by depth  
0-30 or 5 m  
```{r}
xtabs(~Dyp, df_plank)

# Select
sel <- df_plank$Dyp %in% c("0-30 m", "5 m", "5m"); 
df_plank <- df_plank[sel,]

# Stats
cat("Select", sum(sel), "lines\n")
cat(mean(sel)*100, "% of the data")
```

### c. Summarize data  
Use quarters starting with February  
```{r}
df_plank_summ <- df_plank %>%
  mutate(
    Quarter = case_when(
      Month %in% 1:2 ~ 1,
      Month %in% 3:5 ~ 2,
      Month %in% 6:8 ~ 3,
      Month %in% 9:11 ~ 4,
      Month %in% 12 ~ 1),
    Year2 = case_when(
      Month == 12 ~ Year + 1,
      Month < 12 ~ Year),
    Total = Kiselalger + Dinoflagellater + Flagellater
    ) %>%
  group_by(Year2, Quarter) %>%
  summarize_at(.vars = vars(Kiselalger:Flagellater, Total), 
               .funs = funs(med = median, max = max)
               ) %>%
  rename(Year = Year2)

```

### d. Save
```{r}
# write.csv(df_plank_summ, "Data_produced/05_df_plank_summ.csv", row.names = FALSE, quote = FALSE)
```

### e. Plot medians
```{r}
df_plank_summ %>%
  gather("Group", "Median", Kiselalger_med:Total_med) %>%
  mutate(Quarter = paste("Quarter", Quarter)) %>%
  ggplot(aes(Year, Median/1E6)) +
    geom_smooth(method = "lm") + geom_point() +
    facet_grid(Group~Quarter, scales = "free_y")
```
### f. Plot maxima
```{r}
df_plank_summ %>%
  gather("Group", "Maximum", Kiselalger_max:Total_max) %>%
  mutate(Quarter = paste("Quarter", Quarter)) %>%
  ggplot(aes(Year, Maximum/1E6)) +
    geom_smooth(method = "lm") + geom_point() +
    facet_grid(Group~Quarter, scales = "free_y")
```