---
title: Vertical profile of WEP
author: Alex Koiter
---

## Load Libraries
```{r}
library(tidyverse)
library(viridis)
library(patchwork)
```

## Read in data
```{r}
conc <- read_csv(here::here("./notebooks/P_concentration.csv")) %>%
  rename("conc" = "ak_content") %>%
  filter(timing == "Before") %>%
  mutate(sample_type = fct_recode(sample_type, "Ah" = "Soil")) %>%
  mutate(sample_type = fct_relevel(sample_type, c("Ah", "Organic", "Litter","Biomass"))) %>%
  mutate(location = fct_relevel(location, c("Upper", "Middle", "Lower")))

mass_data <- read_csv(here::here("./notebooks/mass.csv")) %>%
  filter(timing == "Before") %>%
  mutate(sample_type = fct_relevel(sample_type, c("Biomass", "Litter")))

bd_data <- read_csv(here::here("./notebooks/bulk density.csv")) %>% # Bulk density measurements
  mutate(mass = bd * length_cm/100 * 1) %>%
  mutate(sample_type = fct_recode(sample_type, "Organic" = "LFH"))
```

## Merge data
```{r}
total_data <- conc %>%
  left_join(mass_data) %>%
  left_join(bd_data) %>%
  mutate(p_total = case_when(sample_type == "Biomass" ~ conc * dryweight/1000 /0.25,
                             sample_type == "Litter" ~ conc * dryweight/1000 /0.25,
                             sample_type == "Organic" ~ conc *mass,
                             sample_type == "Ah"  ~ conc *mass)) %>%
  mutate(sample_type = fct_relevel(sample_type, c("Ah", "Organic", "Litter","Biomass"))) %>%
  mutate(location = fct_relevel(location, c("Upper", "Middle", "Lower")))

```

## Summary
```{r}
conc %>%
  group_by(sample_type) %>%
  summarise(mean = mean(conc, na.rm = T),
            median = median(conc, na.rm = T),
            sd = sd(conc, na.rm = T),
            IQR = quantile(conc, 0.75, na.rm = T) -quantile(conc, 0.25, na.rm = T))

conc %>%
  group_by(sample_type, location) %>%
  summarise(mean = mean(conc, na.rm = T),
            median = median(conc, na.rm = T),
            sd = sd(conc, na.rm = T),
            IQR = quantile(conc, 0.75, na.rm = T) -quantile(conc, 0.25, na.rm = T))

total_data %>%
  group_by(sample_type) %>%
  summarise(mean = mean(p_total, na.rm = T),
            median = median(p_total, na.rm = T),
            sd = sd(p_total, na.rm = T),
            IQR = quantile(p_total, 0.75, na.rm = T) -quantile(p_total, 0.25, na.rm = T))

total_data %>%
  group_by(sample_type, location) %>%
  summarise(mean = mean(p_total, na.rm = T),
            median = median(p_total, na.rm = T),
            sd = sd(p_total, na.rm = T),
            IQR = quantile(p_total, 0.75, na.rm = T) -quantile(p_total, 0.25, na.rm = T))
```


## Plots
Bulk Density
```{r}
p5 <- ggplot(data = bd_data, aes(y = sample_type, x = bd, fill = location)) + 
  geom_boxplot() +
  theme_bw(base_size = 12) +
  labs(x = expression(paste("Bulk Density (", kg~m^{-3}, ")")), tag = "b)") +
  scale_fill_viridis_d(begin = 0.3, end = 1) +
  theme(legend.position = "bottom", 
        legend.title = element_blank(),
        axis.title.y = element_blank())
```

biomass and litter weights

```{r}
p4 <- ggplot(data = mass_data, aes(y = sample_type, x = dryweight/1000 /0.25, fill = location)) + 
  geom_boxplot() +
  theme_bw(base_size = 12) +
  labs(x = expression(paste("Mass (", kg~m^{-2}, ")")), y = "P Source", tag = "a)") +
  scale_fill_viridis_d(begin = 0.3, end = 1) +
  theme(legend.position = "bottom", 
        legend.title = element_blank())
```

```{r}
#| label: suppfig-weights-bd
#| fig-cap: "a) Mass of biomass and litter before grazing and mowing and b) the bulk density of the organic layer and 10 cm Ah horizon"
#| fig-alt: "a) Mass of biomass and litter before grazing and mowing and b) the bulk density of the organic layer and 10 cm Ah horizon"
#| 
p4 + p5 + plot_layout(guides = 'collect') & theme(legend.position = 'bottom', legend.title = element_blank())
```


```{r}
p1 <- ggplot(data = conc, aes(y = sample_type, x = conc, fill = location)) + 
  geom_boxplot() + 
  #scale_x_log10() +
  theme_bw(base_size = 12) +
  theme(legend.position = c(0.2, 0.8),
        legend.title = element_blank()) +
  labs(x = expression(paste("WEP Concentration (", mg~kg^{-1}, ")")), y = "P Source", tag = "a)") +
  scale_fill_viridis_d(name = "Location", begin = 0.3, end = 1)
p1

p2 <- ggplot(data = total_data, aes(y = sample_type, x = p_total, fill = location)) + 
  geom_boxplot() + 
  #scale_x_log10() +
  theme_bw(base_size = 12) +
  theme(axis.title.y = element_blank(),
        legend.position = "none") +
  labs(x = expression(paste("WEP Total (", mg~m^{-2}, ")")), tag = "b)") +
  scale_fill_viridis_d(name = "Location", begin = 0.3, end = 1)
p2
```

```{r}
#| label: fig-vertical-WEP
#| fig-cap: "Vertical and longitudinal profiles of a) WEP concentration and b) WEP content in the riparian areas prior to grazing and mowing treatments."
#| fig-alt: "Vertical and longitudinal profiles of a) WEP concentration and b) WEP content in the riparian areas prior to grazing and mowing treatments."
#|  
p3 <- p1+p2 + plot_layout(guides = 'collect') & theme(legend.position = 'bottom', legend.title = element_blank())  

p3
#ggsave(plot = p3, filename = "Figures/Profile P.png", width = 174, height = 100, units = "mm", dpi = 600) 
```