05b_Annual_and_quarterly_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   
---

Overview:   
* Produces summarized data by year (ending in '_a') and quarter (ending in '_q'), ex df_rivers_summ_a
* Reads data from delivered data files (rivers, hydrology)  
* Data produced are written to Data_produced (starting with '05_') for later use   
* HFR: Annual time seies plots for report (Kap 3.1) are produced here, except plots of soft-bottom (pkt 4c) and hard-bottom (pkt 5) variables which are in script 01_ , 

## 0. Libraries
```{r}
library(tidyverse)
library(readxl)
library(broom)
library(lubridate)
library(zoo)
#install.packages("tidyquant")
#library(tidyquant)
#install.packages("cowplot")
#library(cowplot)
# library(pander)
#?geom_smooth
```

## 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", "TrspTot PO4-P", "TrspTot SiO2", "TrspTot SPM")
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()

# rename kolonner
df_rivers <- rename(df_rivers, TotN = TrspTot_TOTN,
                    NO3 = "TrspTot_NO3-N",
                    NH4 = "TrspTot_NH4-N",
                    TotP = TrspTot_TOTP,
                    TOC = TrspTot_TOC,
                    Alkalinity = TrspTot_ALK,
                    Calcium = TrspTot_Ca,
                    Discharge = DisTot,
                    PO4 = "TrspTot_PO4-P",
                    SPM = TrspTot_SPM,
                    Si = TrspTot_SiO2)

print (df_rivers)
str (df_rivers)

# Table of available data for each river
tb <- df_rivers %>% 
  gather("Variable", Value, TotN:Si) %>%
  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(Discharge, na.rm = TRUE), 
            Q10 = quantile(Discharge, 0.1, na.rm = TRUE), 
            Q90 = max(Discharge, 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

#ggsave ("Figures_rapp/Rivers_monthly_Q_Local.png", width = 8, height = 6, dpi=500)

#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(Discharge, na.rm = TRUE), 
            Q10 = quantile(Discharge, 0.1, na.rm = TRUE), 
            Q90 = max(Discharge, 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

#ggsave ("Figures_rapp/Rivers_monthly_Q_Distant.png", width = 8, height = 6, dpi=500)

# gg + scale_y_log10()
```

### d. Summarize by "local rivers"/"distant rivers"    
    * Data ending with _a = annual, ending in _q = quarter
    * set seasons to match hydro and plankton
    * after discussion w/ Øyvind exclude all "Local rivers" and set following years and variables to NA:  Si for 1990-1995, TOC for 1990-1991
    
```{r}
#str (df_rivers)
df <- df_rivers %>%
  mutate(River_type = 
           ifelse(substr(Station_name, 1, 4) %in% c("Nide","Sønd","Stor"), "Local", "Distant"))


df_rivers_summ_a <- df %>% 
  group_by(River_type, Year) %>%
  summarise_at(c("TotN", "NO3", "TotP", "TOC", "Discharge", "PO4", "SPM", "Si"), sum, na.rm = TRUE)


# NB!season 4 for 1989 gets included...
df_rivers_summ_q <- df %>%
  mutate(Quarter = case_when(
    Month %in% 1 ~ 4,
      Month %in% 2:4 ~ 1,
      Month %in% 5:7 ~ 2,
      Month %in% 8:10 ~ 3,
      Month %in% 11:12 ~ 4),
    Year2 = case_when(
      Month == 1 ~ Year - 1,
      Month > 1 ~ Year)
    ) %>%
  group_by(River_type, Year2, Quarter) %>%
  summarise_at(c("TotN", "NO3", "TotP", "TOC", "Discharge", "PO4", "SPM", "Si"), sum, na.rm = TRUE) %>%
  rename(Year = Year2)


# exclude local rivers and set 0 i transport to NA (ie. Si vars from 1990-1995), se comment above
df_rivers_summ_a <- df_rivers_summ_a %>%
  filter(River_type != "Local") %>% 
  droplevels()%>% 
  mutate_all(funs(ifelse(. == 0, NA, .)))

df_rivers_summ_q <- df_rivers_summ_q %>%
  filter(River_type != "Local") %>% 
  droplevels() %>%
  mutate_all(funs(ifelse(. == 0, NA, .)))

# set transport of TOC in 1990 and 1991 to NA (only sampled in Glomma, yields unrealistic low total transports for these years)
is.na (df_rivers_summ_a$TOC) <- 1:2

is.na (df_rivers_summ_q$TOC) <- 1:9


```

### e1. Plot annual data  
*Also testing out plotting moving average using tidyquant::geom_ma  
```{r}
# geom_ma: sma= simple moving average, rolling mean over a period defined by n
#geom_ma(ma_fun = SMA, n = 5, linetype = 1, size = 1) +

#### USE - annual time series plots for rivers (Kap 3.1)

# Distant - several vars
df_rivers_summ_a %>%
  mutate(Discharge = Discharge/1E3) %>%
  gather("Variable", "Value", TotN:Si) %>%
  #mutate(Var=factor(Var, levels = c("Temperature", "Salinity", "O2", "DIN", "PO4", "Si", "Chla", "TotN", "TotP", "POC", "PON", "POP", "TSM", "Secchi")))  %>%
  mutate(Value = Value/1E3) %>%
  filter(River_type %in% "Distant") %>%
  ggplot(aes(Year, Value)) + 
    geom_point() +
    geom_smooth() +
    facet_wrap(~Variable, scales = "free_y", nrow = 4)

#ggsave ("Figures_rapp/Rivers_a_Distant_all.png", width = 8, height = 6, dpi=500)

#### NOT USE
# Local - several vars
df_rivers_summ_a %>%
  #mutate(Discharge = Discharge/1E3) %>%
  gather("Variable", "Value", TotN:Discharge) %>%
  #mutate(Value = Value/1E3) %>%
  filter(River_type %in% "Local") %>%
  ggplot(aes(Year, Value)) + 
    geom_point() +
    geom_smooth() +
    facet_wrap(~Variable, scales = "free_y", )

#ggsave ("Figures_rapp/Rivers_a_Local_all.png", width = 8, height = 6, dpi=500)


str (df_rivers_summ_a)

# Annual time series plot with one variable - both river types
ggplot(df_rivers_summ_a, aes(Year, Discharge/1E6)) +
  geom_smooth(method = "lm") +
  geom_point() + 
  facet_wrap(~River_type, scales = "free_y") +
  labs(x = "Year", y= "Discharge (/1E6)")

ggplot(df_rivers_summ_a, aes(Year, TotN/1E3)) +
  geom_smooth(method = "lm") +
  geom_point() + 
  facet_wrap(~River_type, scales = "free_y") +
  labs(x = "Year", y= "TotN (/1000)")


# Annual time series with one variable - divided by river type
ggplot(df_rivers_summ_a, aes(Year,TotN/1E3, color = River_type)) +
  geom_smooth(method = "lm") +
  geom_point() + 
  facet_grid(River_type~., scales = "free_y") +
  labs(x = "Year", y= "Transport TotN (/1000)", color = "River type")

ggplot(df_rivers_summ_a, aes(Year, NO3/1E3, color = River_type)) +
  geom_smooth(method = "lm") +
  geom_point() + 
  facet_grid(River_type~., scales = "free_y") +
  labs(x = "Year", y= "Transport NO3 (/1000)", color = "River type")

ggplot(df_rivers_summ_a, aes(Year, TotP/1E3, color = River_type)) +
  geom_smooth(method = "lm") +
  geom_point() + 
  facet_grid(River_type~., scales = "free_y") +
  labs(x = "Year", y= "Transport TotP (/1000)", color = "River type")

ggplot(df_rivers_summ_a, aes(Year, TOC/1E3, color = River_type)) +
  geom_smooth(method = "lm") +
  geom_point() + 
  facet_grid(River_type~., scales = "free_y") +
  labs(x = "Year", y= "Transport TOC (/1000)", color = "River type")


# Annual time series plot with several variables for each river type

# cannot have two couluns with river type due to large difference in values
#df_rivers_summ_a %>%
 # gather("Variable", "Value", TotN:Discharge) %>%
  #mutate(Value = Value/1E6) %>%
  #ggplot(aes(Year, Value)) +
  #  geom_ma(ma_fun = SMA, n = 5, linetype = 1, size = 1) + 
   # geom_point() +
    #facet_grid(Variable~River_type, scales = "free_y", labeller = label_both)
    #labs(x = "Year", y= "Transports")

str(df_rivers_summ_a)

# arrange the two plots next to eachother
#ggg <- plot_grid(gg1, gg2, labels = c('(Distant)', '(Local)'), ncol = 2)
# For better results, increase margin on left side of plot  
#ggg  + theme(plot.margin = margin(12, 6, 6, 32))



```

### e2. Plot quarterly data
```{r}
str(df_rivers_summ_q)


### USE seasonal time series for rivers (Appendix)

# Distant - several vars
df_rivers_summ_q %>%
  mutate(Discharge = Discharge/1E3) %>%
  gather("Var", "Value", Discharge, SPM, TOC) %>%
  mutate(Var=factor(Var, levels = c("Discharge", "SPM", "TOC")))  %>%
  mutate(Value = Value/1E3) %>%
    ggplot(aes(Year, Value)) +
      geom_smooth() + 
      geom_point() +
      facet_grid(Var~Quarter, scales = "free_y")

ggsave ("Figures_rapp/Rivers_q_Dis_SPM_TOC.png", width = 8, height = 6, dpi=500)

df_rivers_summ_q %>%
  gather("Var", "Value", NO3, Si, PO4) %>%
  mutate(Var=factor(Var, levels = c("NO3", "PO4", "Si")))  %>%
  mutate(Value = Value/1E3) %>%
    ggplot(aes(Year, Value)) +
      geom_smooth() + 
      geom_point() +
      facet_grid(Var~Quarter, scales = "free_y")

ggsave ("Figures_rapp/Rivers_q_NO3_PO4_Si.png", width = 8, height = 6, dpi=500)

df_rivers_summ_q %>%
  mutate(Discharge = Discharge/1E3) %>%
  gather("Var", "Value", TotN, TotP) %>%
  mutate(Var=factor(Var, levels = c("TotN", "TotP")))  %>%
  mutate(Value = Value/1E3) %>%
    ggplot(aes(Year, Value)) +
      geom_smooth() + 
      geom_point() +
      facet_grid(Var~Quarter, scales = "free_y")

ggsave ("Figures_rapp/Rivers_q_TotN_TotP.png", width = 8, height = 6, dpi=500)

#### NOT USE

## by river type (local, distant)
# plot with onve variable - divided by river type and quarter
ggplot(df_rivers_summ_q, aes(Year, Discharge)) +
  geom_smooth() +
  geom_point() + 
  facet_grid(River_type~Quarter, scales = "free_y")

#ggsave ("Figures_rapp/Rivers_q_Discharge.png", width = 8, height = 6, dpi=500)

ggplot(df_rivers_summ_q, aes(Year, NO3)) +
  geom_smooth() +
  geom_point() + 
  facet_grid(River_type~Quarter, scales = "free_y")

#ggsave ("Figures_rapp/Rivers_q_NO3.png", width = 8, height = 6, dpi=500)

ggplot(df_rivers_summ_q, aes(Year, TotN)) +
  geom_smooth() +
  geom_point() + 
  facet_grid(River_type~Quarter, scales = "free_y")

#ggsave ("Figures_rapp/Rivers_q_TotN.png", width = 8, height = 6, dpi=500)

ggplot(df_rivers_summ_q, aes(Year, TotP)) +
  geom_smooth() +
  geom_point() + 
  facet_grid(River_type~Quarter, scales = "free_y")

#ggsave ("Figures_rapp/Rivers_q_TotP.png", width = 8, height = 6, dpi=500)

ggplot(df_rivers_summ_q, aes(Year, TOC)) +
  geom_smooth() +
  geom_point() + 
  facet_grid(River_type~Quarter, scales = "free_y")

#ggsave ("Figures_rapp/Rivers_q_TOC.png", width = 8, height = 6, dpi=500)


```

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

print (df_rivers_springflood_allyears)

df_rivers_springflood <- df_rivers %>%
  filter(Month %in% 1:6) %>%
  group_by(Station_name, Year) %>%
  mutate(DisTot_max = max(Discharge), 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(Discharge/DisTot_max_mean*100, na.rm = TRUE),
            DisTot_max_month = Month[Discharge == DisTot_max][1],
            DisTot_40perc = Month[Discharge >= 0.40*DisTot_max_mean][1],
            DisTot_60perc = Month[Discharge >= 0.60*DisTot_max_mean][1],
            DisTot_80perc = Month[Discharge >= 0.80*DisTot_max_mean][1])

print (df_rivers_springflood)

```


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

ggsave ("Figures_rapp/Rivers_Maxflood.png", width = 8, height = 6, dpi=500)

```

### g2. Plot of flood timing
```{r}
df_rivers_springflood %>%
  ggplot(aes(Year, DisTot_max_month)) +
  geom_smooth() + geom_point() +
  facet_wrap(~Station_name)

ggsave ("Figures_rapp/Rivers_timing_flood.png", width = 8, height = 6, dpi=500)

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)
```


### h. River group level: Timing and size of spring flood
```{r}
df_rivergroup <- df_rivers %>%
  filter(!Station_name %in% "Otra") %>%
  mutate(River_type = 
           ifelse(substr(Station_name, 1, 4) %in% c("Nide","Sønd","Stor"), "Local", "Distant")) %>%
  group_by(River_type, Year, Month) %>%
  summarize(DisTot = sum(Discharge, na.rm = TRUE))
  
df_rivergroup_springflood_allyears <- df_rivergroup %>%
  group_by(River_type, Year) %>%
  mutate(DisTot_max = max(DisTot[Month %in% 1:6]), na.rm = TRUE) %>%
  group_by(River_type) %>%
  summarize(DisTot_max_mean = mean(DisTot_max, na.rm = TRUE))
df_rivergroup_springflood_allyears

df_rivergroup_springflood <- df_rivergroup %>%
  filter(Month %in% 1:6) %>%
  group_by(River_type, Year) %>%
  mutate(DisTot_max = max(DisTot), na.rm = TRUE) %>%
  group_by(River_type) %>%
  mutate(DisTot_max_mean = mean(DisTot_max, na.rm = TRUE)) %>%
  ungroup() %>%
  group_by(River_type, 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])

```

### g1. Plot of max flood
```{r}
ggplot(df_rivergroup_springflood, aes(Year, DisTot_max_rel)) +
  geom_smooth() + geom_point() +
  facet_wrap(~River_type)
```

### g2. Plot of flood timing
```{r}
df_rivergroup_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_jitter(width = 0, height = 0.15) +
  facet_wrap(~River_type)
```


### h. Save all
```{r}
write.csv(df_rivers_summ_a, "Data_produced/05_df_rivers_summ_a.csv", row.names = FALSE, quote = FALSE)
write.csv(df_rivers_summ_q, "Data_produced/05_df_rivers_summ_q.csv", row.names = FALSE, quote = FALSE)
#write.csv(df_rivergroup_springflood, "Data_produced/05_df_rivergroup_springflood.csv", 
         # row.names = FALSE, quote = FALSE)
```


## 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 February, as with plankton data (section 4c)  
    * For quarters, also see script 04, plot in section 3b  
```{r}
df <- Df.Arendal %>%
  mutate(
    Depth = case_when(
      Depth %in% c(0,5,10) ~ "Surface",
      Depth %in% c(20,30) ~ "Intermediate",
      Depth %in% c(50,75) ~ "Deep")
    )

str (df)

# rename kolonner
df <- rename(df, Temperature = Temperatur,
                    Salinity = Salt,
                    Density = Tetthet,
                    Si = SiO4,
                    Chla = Klorofyll,
                    Secchi = Siktdyp)

df_hydro_summ_a <- df %>%
  group_by(Year, Depth) %>%
  summarize_at(vars(Temperature:Secchi), mean, na.rm = TRUE)

df_hydro_summ_q <- df %>%
  mutate(
    Quarter = case_when(
      Month %in% 1 ~ 4,
      Month %in% 2:4 ~ 1,
      Month %in% 5:7 ~ 2,
      Month %in% 8:10 ~ 3,
      Month %in% 11:12 ~ 4),
    Year2 = case_when(
      Month == 1 ~ Year - 1,
      Month > 1 ~ Year)
    ) %>%
  group_by(Year2, Quarter, Depth) %>%
  summarize_at(vars(Temperature:Secchi), mean, na.rm = TRUE) %>%
  rename(Year = Year2)

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

# summere til DIN
df_hydro_summ_a$DIN <- df_hydro_summ_a$NO2_NO3 + df_hydro_summ_a$NH4
df_hydro_summ_q$DIN <- df_hydro_summ_q$NO2_NO3 + df_hydro_summ_q$NH4

```

### b2. Gap in TSM - plot
```{r}
ggplot(df_hydro_summ_a, aes(Year, TSM)) +
  geom_line() + geom_point() +
  facet_grid(.~Depth)

ggplot(df_hydro_summ_q, aes(Year, TSM)) +
  geom_line() + geom_point() +
  facet_grid(Quarter~Depth)
```

### b3. Gap in TSM - interpolate data    
Using na.approx from the zoo package  
```{r}

for (dep in unique(df_hydro_summ_a$Depth)){
  sel <- df_hydro_summ_a$Depth == dep
  x <- zoo(df_hydro_summ_a$TSM[sel])
  df_hydro_summ_a$TSM[sel]  <- as.numeric(na.approx(x))
}


for (dep in unique(df_hydro_summ_q$Depth)){
  for (qu in 1:4){
    sel <- with(df_hydro_summ_q, Depth == dep & Quarter == qu)
    x <- zoo(df_hydro_summ_q$TSM[sel])
    df_hydro_summ_q$TSM[sel]  <- as.numeric(na.approx(x))
  }
}

```

### b3. Gap in TSM - plot interpolation
```{r}
df_hydro_summ_a %>% 
  mutate(Interpol = Year %in% 2012:2013) %>%
  ggplot(aes(Year, TSM)) +
  geom_line() + geom_point(aes(color = Interpol)) +
  scale_color_manual(values = c("black","red")) +
  facet_grid(.~Depth)

ggsave ("Figures_rapp/Hydro_interp_TSM.png", width = 8, height = 6, dpi=500)

df_hydro_summ_q %>% 
  mutate(Interpol = Year %in% 2012:2013) %>%
  ggplot(aes(Year, TSM)) +
  geom_line() + geom_point(aes(color = Interpol)) +
  scale_color_manual(values = c("black","red")) +
  facet_grid(Quarter~Depth)
```

### c. Save
```{r}
write.csv(df_hydro_summ_a, "Data_produced/05_df_hydro_summ_a.csv", row.names = FALSE, quote = FALSE)
write.csv(df_hydro_summ_q, "Data_produced/05_df_hydro_summ_q.csv", row.names = FALSE, quote = FALSE)
```

### d1. Plot annual data
```{r}

str(df_hydro_summ_a)

# annual time series - single variable
ggplot(df_hydro_summ_a, aes(Year, Temperatur, color = Depth)) +
  geom_smooth() + 
  geom_point()

# annual time series - several variables
df_hydro_summ_a %>%
  gather("Var", "Concentration", Temperature, Salinity, O2) %>%
  mutate(Var=factor(Var, levels = c("Temperature", "Salinity", "O2")))  %>%
  ggplot(aes(Year, Concentration, color = Depth)) +
    geom_smooth() +
    geom_point() +
    facet_wrap(~Var, scales = "free_y")

ggsave ("Figures_rapp/Hydro_a_T_S_O.png", width = 8, height = 6, dpi=500)

df_hydro_summ_a %>%
  gather("Var", "Concentration", DIN, PO4, Si) %>%
  #mutate(Var=factor(Var, levels = c("DIN", "PO4", "Si")))  %>%
  ggplot(aes(Year, Concentration, color = Depth)) +
    geom_smooth() +
    geom_point() +
    facet_wrap(~Var, scales = "free_y")

ggsave ("Figures_rapp/Hydro_a_InorgNut.png", width = 8, height = 6, dpi=500)

#HFR: rev siktedyp
df_hydro_summ_a %>%
  gather("Var", "Concentration", Chla, TotP, TotN, Secchi) %>%
  mutate(Var=factor(Var, levels = c("Chla", "TotN", "TotP", "Secchi")))  %>%
  ggplot(aes(Year, Concentration, color = Depth)) +
    geom_smooth() +
    geom_point() +
    facet_wrap(~Var, scales = "free_y")

ggsave ("Figures_rapp/Hydro_a_Chla_TotNP.png", width = 8, height = 6, dpi=500)

df_hydro_summ_a %>%
  gather("Var", "Concentration", POC, PON, POP, TSM) %>%
  #mutate(Var=factor(Var, levels = c("Temperatur", "Salt", "O2")))  %>%
  ggplot(aes(Year, Concentration, color = Depth)) +
    geom_smooth() +
    geom_point() +
    facet_wrap(~Var, scales = "free_y")

ggsave ("Figures_rapp/Hydro_a_POM.png", width = 8, height = 6, dpi=500)

# "samleplot" med alle var

df_hydro_summ_a %>%
  gather("Var", "Concentration", Temperature, Salinity, O2, PO4, Si, Chla:Secchi) %>%
  mutate(Var=factor(Var, levels = c("Temperature", "Salinity", "O2", "DIN", "PO4", "Si", "Chla", "TotN", "TotP", "POC", "PON", "POP", "TSM", "Secchi")))  %>%
  ggplot(aes(Year, Concentration, color = Depth)) +
    geom_smooth() +
    geom_point() +
    facet_wrap(~Var, scales = "free_y", ncol = 3)

```

### d2. Plot quarterly data
```{r}
#str (df_hydro_summ_q)

# heller bruke sma?
#til geom_smooth: method = "lm"

# timeseries - divided by quarters - single var
ggplot(df_hydro_summ_q, aes(Year, Temperatur, color = Depth)) +
  geom_smooth() + 
  geom_point() +
  facet_grid(.~Quarter, labeller = label_both)

# several vars
df_hydro_summ_q %>%
  gather("Var", "Concentration", Temperature, Salinity, O2) %>%
   mutate(Var=factor(Var, levels = c("Temperature", "Salinity", "O2")))  %>%
  ggplot(aes(Year, Concentration, color = Depth)) +
    geom_smooth() + 
    geom_point() +
    facet_grid(Var~Quarter, scales = "free_y")

ggsave ("Figures_rapp/Hydro_q_T_S_O.png", width = 8, height = 6, dpi=500)

df_hydro_summ_q %>%
  gather("Var", "Concentration", DIN, PO4, Si) %>%
   #mutate(Var=factor(Var, levels = c("DIN", "PO4", "Si")))  %>%
  ggplot(aes(Year, Concentration, color = Depth)) +
    geom_smooth() + 
    geom_point() +
    facet_grid(Var~Quarter, scales = "free_y")

ggsave ("Figures_rapp/Hydro_q_InorgNut.png", width = 8, height = 6, dpi=500)

df_hydro_summ_q %>%
  gather("Var", "Concentration", Chla, TotP, TotN, Secchi) %>%
   mutate(Var=factor(Var, levels = c("Chla", "TotN", "TotP", "Secchi")))  %>%
  ggplot(aes(Year, Concentration, color = Depth)) +
    geom_smooth() + 
    geom_point() +
    facet_grid(Var~Quarter, scales = "free_y")

ggsave ("Figures_rapp/Hydro_q_Chla_TotNP.png", width = 8, height = 6, dpi=500)

df_hydro_summ_q %>%
  gather("Var", "Concentration", POC, PON, POP, TSM) %>%
  #mutate(Var=factor(Var, levels = c("Klorofyll", "TotN", "TotP", "O2")))  %>%
  ggplot(aes(Year, Concentration, color = Depth)) +
    geom_smooth() + 
    geom_point() +
    facet_grid(Var~Quarter, scales = "free_y")

ggsave ("Figures_rapp/Hydro_q_POM.png", width = 8, height = 6, dpi=500)

```

## 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  
As hydrological data (section 3b), use quarters starting with February (see script 04, plot in section 3b) 
```{r}
df_plank_summ_a <- df_plank %>%
  mutate(Total = Kiselalger + Dinoflagellater + Flagellater) %>%
  group_by(Year) %>%
  summarize_at(.vars = vars(Kiselalger:Flagellater, Total), 
               .funs = funs(med = mean, max = max)
               )

df_plank_summ_q <- df_plank %>%
  mutate(
    Quarter = case_when(
      Month %in% 1 ~ 4,
      Month %in% 2:4 ~ 1,
      Month %in% 5:7 ~ 2,
      Month %in% 8:10 ~ 3,
      Month %in% 11:12 ~ 4),
    Year2 = case_when(
      Month == 1 ~ Year - 1,
      Month > 1 ~ Year),
    Total = Kiselalger + Dinoflagellater + Flagellater
    ) %>%
  group_by(Year2, Quarter) %>%
  summarize_at(.vars = vars(Kiselalger:Flagellater, Total), 
               .funs = funs(med = mean, max = max)
               ) %>%
  rename(Year = Year2)

# rename kolonner
df_plank_summ_a <- rename(df_plank_summ_a, Diatoms = Kiselalger_med,
                    Dinoflagellates = Dinoflagellater_med,
                    Flagellates = Flagellater_med
                    )

df_plank_summ_q <- rename(df_plank_summ_q, Diatoms = Kiselalger_med,
                    Dinoflagellates = Dinoflagellater_med,
                    Flagellates = Flagellater_med
                    )

```

### d. Save
```{r}
write.csv(df_plank_summ_a, "Data_produced/05_df_plank_summ_a.csv", row.names = FALSE, quote = FALSE)
write.csv(df_plank_summ_q, "Data_produced/05_df_plank_summ_q.csv", row.names = FALSE, quote = FALSE)
```


### e1. Plot medians by year
```{r}

str(df_plank_summ_a)

df_plank_summ_a %>%
  gather("Group", "Value", Diatoms:Flagellates) %>%
  ggplot(aes(Year, Value)) +
    geom_smooth() + geom_point() +
    facet_grid(Group~., scales = "free_y")

ggsave ("Figures_rapp/Plankton_a_Groups.png", width = 8, height = 6, dpi=500)

```


### e2. Plot medians by quarter
```{r}
df_plank_summ_q %>%
  gather("Group", "Value", Diatoms:Flagellates) %>%
  #mutate(Quarter = paste("Quarter", Quarter)) %>%
  ggplot(aes(Year, Value), ylim(c(0, NA))) +
    geom_smooth() + geom_point() +
    facet_grid(Group~Quarter, scales = "free_y")

ggsave ("Figures_rapp/Plankton_q_Groups.png", width = 8, height = 6, dpi=500)
```
### f. Plot maxima
```{r}
df_plank_summ_q %>%
  gather("Group", "Maximum", Diatoms:Flagellates) %>%
  mutate(Quarter = paste("Quarter", Quarter)) %>%
  ggplot(aes(Year, Maximum/1E6)) +
    geom_smooth(method = "lm") + geom_point() +
    facet_grid(Group~Quarter, scales = "free_y")
```