plot-proc-dendro-PacoEzpela.Rmd

---
output:
  pdf_document: default
  html_document: default
---
```{r include=FALSE}
source("lib-dendro.R")

library(ggplot2)
theme_set(theme_bw())
library(ggnewscale)
library(dplyr)
library(lubridate)
library(glue)


### DEFINE GLOBAL VARS ###
PLACE = 'PacoEzpela'
DATA_DIR = glue('processed/{PLACE}-processed')
ENV_DIR = glue('processed/{PLACE}-env-processed')

SELECTED_TMS <- 94252896

Fagus <- c(92223484, 92232411, 92232420, 92232415, 92232419)
Non_Declining_Abies <- c(92222339, 92232424, 92232421, 92232416, 92232423)
Declining_Abies <- c(92232412, 92232418, 92232414, 92232431, 92232417)

### IMPORT DATA ###

# importing dendro data #
list_files <- list.files(file.path(".",DATA_DIR), pattern="*.csv$", full.names=TRUE)
db<-read.all.processed(list_files)

# importing climate data #
db.env <- read.env.proc(file.path(".",ENV_DIR,'proc-env.csv'))
db.env <- db.env[db.env$series == SELECTED_TMS,]
head(db.env)
```


```{r include=FALSE}
### CLEAN & PREPARE DATA ###
db = db %>% mutate(
                  date = date(db$ts),
                  year = year(date),
                  class = case_when(
                    series %in% Declining_Abies ~ factor("D"),
                    series %in% Non_Declining_Abies ~ factor("ND"),
                    series %in% Fagus ~ factor("European beech"),
                    .default = NA
                    )
                  )
```

Filtering by date:

```{r}
ts_start <- "2023-05-01"
ts_end <- "2023-09-30 23:45:00"
db <- db[which(db$ts>=ts_start & db$ts<=ts_end),] 
db.env <- db.env[which(db.env$ts>=ts_start & db.env$ts<=ts_end),] 
```


```{r}
# INSPECT DATA
str(db)
summary(db)
# head(db)
# tail(db)
```


```{r}
## PLOT ALL DENDROS ##
plot_multiple_dendro <- function (data, title, y = data$value)  {
  ggplot(data = data, mapping = aes(x=ts, y=y, col=series))+
  geom_line( )+
  # ggtitle(paste0("Dendro data for sensor series: ",db$series[1], " - ", db$sp[1])) +
  labs(x=expression('Date'),
       y=expression( Delta*"D (um)") ) +
  geom_hline(yintercept=0,lty=2,linewidth=0.2)+
  scale_x_datetime(date_breaks = "1 month", date_labels = "%m-%y") +
  ggtitle(title) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 30, hjust=1))
}
```


```{r}
plot_multiple_dendro(db, "Plot by dendrometer")
ggsave(glue('output/all-dendrometers-{PLACE}.png'), width = 15, height = 10)
```

# Normalization 0-1 to better compare dendros

```{r}
normalized.db <- db %>%
  select (series, ts, value) %>%
  group_by(series) %>%
  mutate( normalized_value = normalize.0_1(value), .keep = 'all' )
```


```{r}
plot_multiple_dendro(normalized.db, "Dendros data normalized to [0-1]", y = normalized.db$normalized_value)
ggsave(glue('output/normalized-all-dendrometers-{PLACE}.png'), width = 15, height = 10)
```


# By species

```{r}
db.D <- db %>% filter (class == "D")
db.ND <- db %>% filter (class == "ND")
db.Fag <- db %>% filter (class == "European beech")
```


```{r}
# plot_multiple_dendro(db.Abi, "Silver fir dendrometers")
```


```{r}
# plot_multiple_dendro(db.Fag, "Dendrómetros European beech")
```

Plot both in one graph:
```{r}
# Define color scales for each class
blue_colors <- c("blue", "cyan", "aquamarine", "cadetblue", "dodgerblue", "darkblue", "turquoise4", "royalblue2")
green_colors <- c("lightgreen", "green", "darkgreen", "darkolivegreen", "olivedrab", "springgreen")
orange_colors <- c("orange", "orangered", "tomato", "darkorange", "coral", "tan")

plotByClass <- 
  ggplot() +
     ggtitle("Dendrometers painted by color scale according to its specie") +
    geom_line(data = db.Fag, aes(x = ts, y = value, col = series)) +
    scale_color_manual(values = green_colors, name = "European beech") +
    new_scale_color() +
    geom_line(data = db.ND, aes(x = ts, y = value, col = series), ) +
    scale_color_manual(values = blue_colors, name = "Non-Declining Abies Alba") +
    new_scale_color() +
    geom_line(data = db.D, aes(x = ts, y = value, col = series), ) +
    scale_color_manual(values = orange_colors, name = "Declining Abies Alba") +
    new_scale_color() +
    geom_hline(yintercept=0,lty=2,linewidth=0.2) +
    labs(x=expression(''),
         y=expression(Delta*" D (um)"))+
    theme_bw() +
    scale_x_datetime(date_breaks = "1 month", date_labels="%b %Y") +
    theme(axis.text.x = element_text(angle = 30, hjust=1))

ggsave(glue('output/all-dendrometers-by-class-{PLACE}.png'),  plot = plotByClass, width = 15, height = 10)
```

Plot means for each specie in one graph:
```{r}
db.mean.Fag <- db.Fag %>% group_by(ts) %>% dplyr::summarise(n = n(), mean = mean(value, na.rm = T), sd = sd(value, na.rm = T), se = sd / sqrt(n) ) %>% mutate(class = "European beech")
db.mean.D <- db.D %>% group_by(ts) %>% dplyr::summarise(n = n(), mean = mean(value, na.rm = T), sd = sd(value, na.rm = T), se = sd / sqrt(n) ) %>% mutate(class = "Declining Silver fir")
db.mean.ND <- db.ND %>% group_by(ts) %>% dplyr::summarise(n = n(), mean = mean(value, na.rm = T), sd = sd(value, na.rm = T), se = sd / sqrt(n) ) %>% mutate(class = "Non-Declining Silver fir")

db.means <- rbind.data.frame(db.mean.Fag, db.mean.D, db.mean.ND)

plotMeans <- 
  ggplot(data = db.means, mapping = aes(x=ts, y=mean, col = class))+
  ggtitle("Mean of the dendrometer values for the three groups of study") +
  scale_color_manual(values=c("European beech" = "green","Non-Declining Silver fir" = "blue", "Declining Silver fir" = "orange")) +
  scale_fill_manual(values=c("European beech" = "green","Non-Declining Silver fir" = "blue", "Declining Silver fir" = "orange")) +
  geom_line( ) +
  geom_ribbon(aes(ymin=mean-se, ymax=mean+se, fill=class), alpha=0.2, show.legend	= FALSE, linetype = 0) +
  labs(x=expression('Date'),
       y=expression( Delta*"D (um)") ) +
  geom_hline(yintercept=0,lty=2,linewidth=0.2)+
  scale_x_datetime(date_breaks = "1 month", date_labels="%b %Y") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 30, hjust=1))
  
plotMeans
ggsave(glue('output/means-by-class-{PLACE}.png'), plot = plotMeans, width = 15, height = 10)
```

# Tree Water Deficit (TWD)

Calculate means for TWD by sp:
```{r}
db.twd <- db %>% summarise(mean_twd = mean(twd, na.rm = T), n = as.integer(n() / 96), twd_se = sd(twd, na.rm = T) / sqrt(n), .by = c(date, class)) %>% mutate(doy = yday(date))
```


Calculate VWC:

```{r}
vwc <- db.env %>% mutate(date = date(db.env$ts)) %>% summarise(max = max(vwc), .by = date)
# vwc
```


```{r warning=FALSE}
plot.twd <- db.twd %>%
  ggplot () +
  geom_line(aes (x = date, y = mean_twd, color = class)) +
  geom_ribbon(aes (x=date, ymin=mean_twd-twd_se, ymax=mean_twd+twd_se, fill=class), alpha=0.3, show.legend	= FALSE, linetype = 0) +
  scale_color_manual(values=c("European beech" = "green","ND" = "blue", "D" = "orange"), labels = c("European beech" = "European beech","ND" = "Non-Declining Silver fir", "D" = "Declining Silver fir")) +
  scale_fill_manual(values=c("European beech" = "green","ND" = "blue", "D" = "orange")) +
  labs(title = "Means of Tree Water Deficit for the three groups of study", x = "Date", y = expression(paste("Tree Water Deficit (", mu, "m)")))+
  scale_x_date(date_breaks = "1 month", date_labels="%b %Y")+
  scale_y_continuous(limits=c(0, 150), sec.axis = sec_axis(trans = ~ . / 1.5, name = "Volumetric Water Content (%)")) +
  geom_line(data = vwc, aes(x = date, y = max*100, linetype = "Volumetric Water Content (%)"), col = "gray28", alpha = 0.9, show.legend	= c(colour = FALSE, linetype = T)) +
  scale_linetype_manual(NULL, values = 4) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 30, hjust=1))
plot.twd
ggsave(glue('output/twd-evolution-by-class-{PLACE}.png'), plot = plot.twd, width = 15, height = 10)
```

```{r}
boxplot(mean_twd ~ class, data = db.twd)
```