01_Check_data.Rmd

---
title: "Overview and some plots of the data"
output: 
  html_document:
    keep_md: true
    toc: true
    toc_depth: 3
    toc_float: true
    code_folding: hide
    df_print: paged
---

Getting an overview of the data   
* *NOTE:* Tables in .md file edited manually in order to be nicely rendered in Github (mostly for testing/learning...)  

## 0. Libraries
```{r}
library(tidyverse)
library(readxl)
library(pander)
#library(gridExtra)
library(grid)
library (caTools)
```

## 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   

***Helene comments:*** 
* Relevante variabler for plot og test av tidsserie:
* Bruke datasettene med monthly and annual loads, og spesielt variablene Discarge (DisTot) og transport av inorg N (NO3+NO4), PO4, TotN, TotP, TOC og SPM. Alkalinitet og calcium er kun for sur nedbør elvene (Nidelva, Storelva og Gjerstadelva), og kan vi kanskje heller komme tilbake til.
*dropp Otra (som jo er "nedstrøms" fra Arendal)
* aggregering: i sesong (som du jo har startet på allerede), er det også mulig å hente ut mnd for vårflom, slik som mnd for våroppblomstring? Dette vil jo sikkert variere for stasjonene (?), men interessant å se om når vårflommen intreffer har endret seg, og om det evn har noen effekt på spesielt hardbunn.

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; example: Storelva loads
```{r}
df <- read_excel("Datasett/River data (from OKA)/Monthly loads/Storelva_monthly loads.xlsx")  
# df <- read_excel("Datasett/River data (from OKA)/Monthly loads/Storelva_monthly loads.xlsx", skip = 1)  
# summary(df)

df <- as.data.frame(df)[-1,]
# colnames(df) %>% 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[,var] <- as.numeric(df[,var])
df$Time <- with(df, lubridate::ymd(paste(Year, Month, "15")))

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

tb <- df %>% 
  gather("Variable", Value, TrspTot_TOTN:DisTot) %>%
  filter(!is.na(Value)) %>%
  xtabs(~Year + Variable, .)
pandoc.table(tb, style = "rmarkdown")
```

### b. Plot Storelva loads  
```{r}
gg <- df %>%
  gather("Variable", Value,  TrspTot_TOTN:DisTot) %>%
  ggplot(aes(Time, Value)) + 
    geom_line() +
    facet_wrap(~Variable, scales = "free_y")
gg
```

### c. Plot annual flow  
```{r}
df <- read_excel("Datasett/River data (from OKA)/Annual mean flows/Mean annual flow.xlsx")
# df <- read_excel("Datasett/River data (from OKA)/Monthly loads/Storelva_monthly loads.xlsx", skip = 1)  
 summary(df)

df <- as.data.frame(df)[-1,]
colnames(df) <- sub(" ", "_", colnames(df), fixed = TRUE)
# colnames(df) %>% dput()
vars <- c("DisMean")
for (var in vars)
  df[,var] <- as.numeric(df[,var])

# Flow
ggplot(df, aes(Year, DisMean)) + geom_line() + facet_wrap(~Station_name) + labs (y = "Annual mean flow")

```

## 3. Hydrografi  


***Helene comments:*** 
* Dag (du kjenner jo disse dataene også godt): egentlig interessant å se på plot og trender av alle disse variablene over tid...
* aggregering: sesong og årlig, OG i 3 dybdelag: surface 0-10m, intermediate 20-30m, deep 50-75.

### a. Data  
Tables of coverages of depth and seasons  
```{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()

# summary(Df.Arendal)
tb <- xtabs(~Depth + Month, Df.Arendal)
pandoc.table(tb, style = "rmarkdown")

tb <- xtabs(~Year + Month, Df.Arendal)
pandoc.table(tb, style = "rmarkdown")
```

### b. Plot salinity, depth 10 m
```{r}
ggplot(Df.Arendal %>% filter(Depth == 10), aes(Dato, Salt)) +
  geom_line() + geom_point()
```

### c. Salinity, time series by month and depth
```{r}
Df.Arendal %>%
  group_by(Year, Month, Depth) %>%
  summarise_all(mean) %>%
  ggplot(aes(Year, Salt)) + 
    geom_line() + 
    facet_grid(Month~Depth)
```


## 4. Soft bottom fauna  

***Helene comments:*** 
* plot og test av trend spesielt relevant for No.species, No.ind, NQI1, H, TOC (alle i df_blot_ind), har laget noen plot i pkt c under. 
* diskuterte litt med Guri om hvordan vi evn skal plotte trender i artsammensetning (artslistene) over tid, dette går jo inn i ordinasjonsanalysen, og Guri mente det kanskje var best (hvertfall for hardbunn) å heller hente ut tidskorrelerte endringer i artsammensetning fra DCA1-aksen. 

### a. Data
* after instruction from Hilde we removed the TOC data from 2014 and sediement fine fraction (63um) for 2012 (measurements errors due to chanage in laboratory). 
```{r}
# species list station B35
df_blot_b35 <- read_excel("Datasett/Bløtbunn/Klimaoverblikk bløtbunn_data til Helene og Dag.xlsx", sheet = "B35_artsliste")
colnames(df_blot_b35)[1] <- "Species"

# species list station B05
df_blot_b05 <- read_excel("Datasett/Bløtbunn/Klimaoverblikk bløtbunn_data til Helene og Dag.xlsx", sheet = "B05_artsliste")
colnames(df_blot_b05)[1] <- "Species"

# number of species and years
cat("Number of species:", nrow(df_blot_b35), "\nNumber of years:", ncol(df_blot_b35), "\n")

# indices and sediment parameters
df_blot_ind <- read_excel("Datasett/Bløtbunn/Klimaoverblikk bløtbunn_data til Helene og Dag.xlsx", sheet = "indekser_sedimentparametere")

# rename variables
df_blot_ind <- rename (df_blot_ind, H = "H'", No.species = "No. species", No.ind = "No. ind.", Fine.Fraction = "%<63um")

# setter TOC i 2014 til NA
df_blot_ind <- df_blot_ind %>%
  mutate(TOC = ifelse(Year == 2014, NA, TOC))

# Setter Fine.Fraction i 2012 til NA  
df_blot_ind <- df_blot_ind %>%
  mutate(Fine.Fraction = ifelse(Year == 2012, NA, Fine.Fraction))

# take mean of individual grabs pr year and station
df_blot_ind <- df_blot_ind %>%
  group_by(STAS, Year) %>%
  summarise_at(vars (H:TOC_norm), funs (median))



```

### b. boxplot of the 10% most abundant species pr station
```{r}

# station B35
df <- df_blot_b35 %>%
  gather("Year", "N", -Species) %>%
  mutate(N = ifelse(is.na(N),0,N)) %>%
  group_by(Species) %>%
  mutate(Median_per_species = median(N))

cat("Percentiles of median abundance per species:\n")
quantile(df$Median_per_species, (0:10)/10, na.rm = TRUE)

gg1 <- df %>%
  filter(Median_per_species >= 2) %>%
  ggplot(aes(Species, N)) + geom_boxplot() + theme(axis.text.x = element_text(angle = 90, hjust = 1))
gg1

# station B05
df <- df_blot_b05 %>%
  gather("Year", "N", -Species) %>%
  mutate(N = ifelse(is.na(N),0,N)) %>%
  group_by(Species) %>%
  mutate(Median_per_species = median(N))

cat("Percentiles of median abundance per species:\n")
quantile(df$Median_per_species, (0:10)/10, na.rm = TRUE)

# 
gg2 <- df %>%
  filter(Median_per_species >= 4) %>%
  ggplot(aes(Species, N)) + geom_boxplot() + theme(axis.text.x = element_text(angle = 90, hjust = 1))
gg2


```

### c. Plots of Indices and sediment parameters vs time
```{r}
str(df_blot_ind)

# individual variables
#ggplot(df_blot_ind, aes(Year, H, NQI1)) + geom_line() + facet_wrap(~STAS) 


# several variables - stations separate
df_blot_ind %>%
  gather("Variable", "Value", H:TOC) %>%
  filter(STAS %in% "B05") %>%
  mutate(Variable=factor(Variable, levels = c("No.species", "No.ind", "H", "NQI1", "TOC", "Fine.Fraction")))  %>%
  ggplot(aes(Year, Value)) +
    geom_point() +
    geom_smooth() +
    facet_wrap(~Variable, scales = "free_y", nrow = 3)
    #labs(x = "Year", y= "Transports")
 
ggsave ("Figures_rapp/Soft_B05_all.png", width = 8, height = 6, dpi=500)

df_blot_ind %>%
  gather("Variable", "Value", H:TOC) %>%
  filter(STAS %in% "B35") %>%
  mutate(Variable=factor(Variable, levels = c("No.species", "No.ind", "H", "NQI1", "TOC", "Fine.Fraction")))  %>%
  ggplot(aes(Year, Value)) +
    geom_point() +
    geom_smooth() +
    facet_wrap(~Variable, scales = "free_y", nrow = 3)
    #labs(x = "Year", y= "Transports")
 
ggsave ("Figures_rapp/Soft_B35_all.png", width = 8, height = 6, dpi=500)

```


### d. Save
```{r}
write.csv(df_blot_ind, "Data_produced/01_df_blot_ind.csv", row.names = FALSE, quote = FALSE)

```

```{r soft bottom indices for publication}
# Stations together for publication
df_blot_ind %>%   
  rename(`Shannon-Wiener` = H) %>% 
  mutate(Station = ifelse(STAS == "B05", "BR1", "BT44")) %>% 
  rename(Species = `No.species`, Individuals = `No.ind`, `Fine fraction` = `Fine.Fraction`) %>% 
  gather("Variable", "Value", `Shannon-Wiener`:TOC) %>%
  mutate(Variable = factor(Variable, 
                           levels = c("Species", "Individuals", "Shannon-Wiener", 
                                      "NQI1", "TOC", "Fine fraction")))  %>%
  ggplot(aes(Year, Value, color = Station, fill = Station)) +
    geom_point(size = 0.9) +
    geom_smooth(size = 0.7, aes(linetype = Station))+
    facet_wrap(~Variable, scales = "free_y", ncol = 2) +
    labs(x = "", y= "") +
  theme_bw()

ggsave ("Figures_publ/Softbottom_indicies.png", width = 8, height = 6, dpi=500)

```

 
## 5. Hard-bottom flora and fauna 


***GSA comments:*** 
* Vi kutter 416 og 417, siden tidsserien her var relativt kort. *(Data er fortsatt med i parametere utenom ordinasjonsaksene)*
* Vi kutter strandsone og konsentrerer oss bare om transektdata også i ordinasjonene. Endret fra tidligere versjon av HBanalysesett.csv
* I ordinasjonene er det både alger og dyr
* Jeg skal forsøke å finne en god måte å plotte artsscorer sammen med punktene i DCA-plot slik at vi kan tolke ordinasjonsresultater også med tanke på hvilke arter som utgjør endringene i sammfunnsstruktur over tid (*finnes nå i ordinasjons-notebooken (02_)*).
* La til figuren fra hardbunnsskriptet over endringer i nedre voksedyp for de ni indeksartene (makroalger) over tid.
* Har lagt til antall arter i gruppene dyr, rød-, grønn- og brunalger over tid
** NB: Vi må snakke med Janne om hvor mye utbytting av personell over tid kan ha å si for resultatene!**


### a. Data (already chewed a bit on)
```{r}
# dir("Datasett/hardbunn_kopi")
# readLines("Datasett/hardbunn_kopi/HBanalysesett.csv", 2)
df <- readr::read_csv2("Datasett/Hardbunn_KOPI/HBanalysesett.csv")
```

### b. Plot DCA
```{r}
df %>%
  gather("DCA_axis", "Value", DCA1:DCA4) %>%
  ggplot(aes(Year, Value)) +
    geom_line() +
    facet_grid(Site~DCA_axis)
```

### c. Plot NMDS  
```{r}
df %>%
  gather("NMDS_axis", "Value", NMDS1:NMDS2) %>%
  ggplot(aes(Year, Value)) +
    geom_line() +
    facet_grid(Site~NMDS_axis)
```

### d. Plot DCAstrand
```{r}
df %>%
  gather("DCA_axis", "Value", DCA1strand:DCA4strand) %>%
  ggplot(aes(Year, Value)) +
    geom_line() +
    facet_grid(Site~DCA_axis)
```

### e. Plot species - LGD - Nedre voksedyp
```{r, fig.width=9}
df %>%
  gather("Species", "Value", `LGD.Halidrys siliquosa`:`LGD.Rhodomela confervoides`) %>%
  mutate(Species = sub("LDG.", "", Species)) %>%
  ggplot(aes(Year, Value)) +
    geom_line() +
    facet_grid(Site~Species)
```
#### Fjerna to stasjoner, snudd y-akse og plotta sammen. Bedre til visuell presentasjon.
```{r LGD index species}

# Guri: lagt inn dette som eksempel, tatt med Loess smoother som i de andre plottene, og endre aksetittel til engelsk. Kan vi broppe LGD for an hvert artsnavn? etterhvert også bytte stasjonsnavn til de "nye" (407=HT113 og 410=HR104)
LGDfig <- df[df$Site %in% c(407, 410), ] %>%
  gather("Species", "Value", `LGD.Halidrys siliquosa`:`LGD.Rhodomela confervoides`) %>%
  mutate(Species = sub("LGD.", "", Species)) %>%
  ggplot(aes(Year, Value, color = as.character(Site))) +
    geom_point() +
    geom_smooth() +
    facet_wrap(~ Species, nrow=3) +
    scale_y_reverse(limits=c(30,0)) +
    scale_color_hue(labels = c("HT113", "HR104")) +
    labs(x = "", y= "Lower growth depth (m)", color = "Station")

LGDfig

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


df[df$Site %in% c(407, 410), ] %>%
  gather("Species", "Value", `LGD.Halidrys siliquosa`:`LGD.Rhodomela confervoides`) %>%
  mutate(Species = sub("LGD.", "", Species)) %>%
  ggplot(aes(Year, Value)) +
    geom_line() +
    facet_grid(Site~Species)

```

MSMDI EQR

```{r Hard-bottom diversity}

MSMDIfig <- df[df$Site %in% c(407, 410), ] %>%
  ggplot(aes(Year, EQR, color = as.character(Site))) +
    geom_point() +
    geom_smooth() +
#    ylim(0, 1) +
    scale_color_hue(labels = c("HT113", "HR104")) +
    labs(x = "", y= "MSMDI-index (EQR-value)", color = "Station")

ggsave(MSMDIfig, "Figures_rapp/MSMDI_EQR.png", width = 8, height = 6, dpi=500)
```
```{r, fig.height = 8, fig.width=8}
grid.arrange(LGDfig, MSMDIfig, nrow=2)
```


### f. Plot number of species within different groups
***GSA comments:***
* Probably makes more sense to plot porportions over time though!
```{r Hard-bottom species groups}

# str(df)

# Guri: Tar også med dette plottet som eksempel

df[df$Site %in% c(407, 410), ] %>%
  gather("Group", "Value", 36:40) %>%
  ggplot(aes(Year, Value, color = as.character(Site))) +
    geom_point() +
    geom_smooth()+
    facet_wrap(~ Group, scales = "free_y") +
    scale_color_hue(labels = c("HT113", "HR104")) +
    labs(x = "", y= "Number of species", color = "Station")

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

plot(df$Green/(df$Brown) ~ df$Year)

```

### g. Plot diversity indices

```{r Hard-bottom diversity}
df[df$Site %in% c(407, 410), ] %>%
  gather("Group", "Value", `H`:`J`) %>%
  ggplot(aes(Year, Value, color = as.character(Site))) +
    geom_point() +
    geom_smooth()+
    facet_wrap(~ Group, scales = "free_y") +
    scale_color_hue(labels = c("HT113", "HR104")) +
    labs(x = "", y= "", color = "Station")

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

```

### h. Plot numbers and diversity indices - for publication

```{r Hard-bottom for publ}
df[df$Site %in% c(407, 410), ] %>% 
  dplyr::select(Year, Site, Animals, antMakroalger:`EQR`) %>% 
  rename(Macroalgae = antMakroalger,
         MSMDI = EQR,
         `Shannon-Wiener` = H,
         `Pielou's evenness` = J,
         `Species` = S) %>% 
  gather("Group", "Value", Animals:MSMDI) %>%
  mutate(Station = ifelse(Site == 407, "HT113", "HR104")) %>% 
  mutate(Group = factor(Group, levels = c("Species", "Animals", "Macroalgae", 
                                          "Shannon-Wiener", "Pielou's evenness", "MSMDI"))) %>% 
  ggplot(aes(Year, Value, color = Station, fill = Station)) +
    geom_point(size = 0.9) +
    geom_smooth(size = 0.7, aes(linetype = Station))+
    facet_wrap(~ Group, scales = "free", ncol = 2) +
    labs(x = "", y= "") +
  theme_bw()

ggsave ("Figures_publ/Hardbottom_indicies.png", width = 8, height = 6, dpi=500)

```


## 6. Phytoplankton

***Helene comments:*** 
* du kjenner jo disse dataene også godt. Vi har jo bare dataene aggregert til klasser og grupper (altså ikke fulle artslister). du har jo allerede plottet tidsserier for gruppene (det holder kanskje?), og teste trend over tid i disse?
* aggregering: sesong, og hente ut artsammensetning ved maks chla på våren til ordinasjon. Plottene du hadde gjort tidligere viser jo dette med at høstoppblomstringer ser ut til å ha forsvunnet. 

### a. Data   
* Note: deleted the empty second row manually, and one empty column 
* Also copied the data (just the range with data, A1:V471) to a new sheet
* NB! filter on depth to only include obs at 5m after 2000 (se below)
```{r}
df <- read_excel("Datasett/Plankton/Planteplankton Arendal.xlsx") # range = "A1:V471"
df$Year <- lubridate::year(df$Dato)
df$Month <- lubridate::month(df$Dato)
xtabs(~Year + Month, df)

# They used depth profiles (0-30m) before 2000, after 2000 they mostly sampled at 5m, but there are some dates where there are profiles after 2000 as well. See more info in data description.
df <- df %>% 
  select(Dato:Month) %>%
  filter(Dyp == "0-30 m" | Dyp == "5 m" | Dyp == "5m")

str (df)
summary (df)

```

### b. Plot time series of main groups

```{r, fig.height=3.5}
df %>%
  gather("Group", "Value", Kiselalger:Flagellater) %>%
  ggplot(aes(Dato, Value)) +
    geom_line() +
    facet_wrap(~Group, scale = "free_y")
```

### c. Plot time series of main groups by month
```{r}
df %>%
  gather("Group", "Value", Kiselalger:Flagellater) %>%
  group_by(Year, Month, Group) %>%
  summarize(Value = mean(Value)) %>%
  filter(Month %in% c(2,4,6,8,10)) %>%
   ggplot(aes(Year, Value)) +
    geom_line() +
    facet_grid(Group~Month, scale = "free_y")
```

### d. Plot time series of 6 main taxa (by median abundance) by month
```{r, fig.height=6}
mean_abundance <- apply(df %>% select(Cyanophycea:Chlorophycea), 2, mean, na.rm = TRUE) %>% sort() %>% rev()
median_abundance <- apply(df %>% select(Cyanophycea:Chlorophycea), 2, median, na.rm = TRUE) %>% sort() %>% rev()
taxa <- names(median_abundance) %>% head(6)
df %>%
  gather("Taxon", "Value", Cyanophycea:Chlorophycea) %>%
  group_by(Year, Month, Taxon) %>%
  summarize(Value = mean(Value)) %>%
  filter(Month %in% c(2,4,6,8,10) & Taxon %in% taxa) %>%
  ggplot(aes(Year, Value)) +
    geom_line() +
    facet_grid(Taxon~Month, scale = "free_y")
```

### e. Plot time series of 5 lesser taxa (by median abundance) by month
```{r, fig.height=6}
taxa <- names(median_abundance) %>% tail(5)
df %>%
  gather("Taxon", "Value", Cyanophycea:Chlorophycea) %>%
  group_by(Year, Month, Taxon) %>%
  summarize(Value = mean(Value)) %>%
  filter(Month %in% c(2,4,6,8,10) & Taxon %in% taxa) %>%
  ggplot(aes(Year, Value)) +
    geom_line() +
    facet_grid(Taxon~Month, scale = "free_y")
```