plot-growth.Rmd

---
title: "R Notebook"
output: html_notebook
---

```{r}
source('lib-dendro.R')
source('lib-ts-analysis.R')
library(ggplot2)
library(tidyverse)
library(glue)
library(patchwork)
library(imputeTS)
```

```{r}
TreeList <- read.table("TreeList.txt", header=T)

PLACE = 'Penaflor'
DATA_DIR = glue('processed/{PLACE}-processed')
ENV_DIR = glue('processed/{PLACE}-env-processed')


SELECTED_TMS <- case_when(  # select one TMS sensor for this site which has all its data valid
  PLACE == 'Miedes' ~ "94231940",
  PLACE == 'Corbalan' ~ "94231943",
  PLACE == 'Penaflor' ~ "94231950",
  .default = NA
)

SAVE <- T

## STUDY PERIOD ##
ts_start2022 <- "2022-04-01 09:00:00" # from first of April
ts_end2022 <-"2022-11-30 23:45:00" # to 30 November

ts_start2023 <- "2023-04-01 09:00:00" # from first of April
ts_end2023 <-"2023-11-30 23:45:00" # to 30 November

```

Load dataset:
```{r}
list_files <- list.files(file.path(".",DATA_DIR), pattern="*.csv$", full.names=TRUE)
db<-read.all.processed(list_files)
head(db)

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

Declare dendros to analyse:
```{r}
if (PLACE == 'Miedes') {
  SELECTED_DENDROS <- unique(db$series)[unique(db$series) != 92222175] # All but 92222175 
} else {
  SELECTED_DENDROS <- unique(db$series)
}
```

Define classes which each series belong to:
```{r}
Qi = TreeList %>% filter(class == "Quercus", series %in% unique(db$series)) %>% pull(series)
P_D = TreeList %>% filter(class == "D", series %in% unique(db$series)) %>% pull(series)
P_ND = TreeList %>% filter(class == "ND", series %in% unique(db$series)) %>% pull(series)
```

Prepare db
```{r}
db = db %>% mutate(
                  date = date(db$ts),
                  year = year(date),
                  class = case_when(
                    series %in% Qi ~ factor("Quercus"),
                    series %in% P_D ~ factor("D"),
                    series %in% P_ND ~ factor("ND"),
                    .default = NA
                    )
                  ) %>%
            filter(!is.na(class))
db
```

Set ggplot label options, themes & vars (Note that some of the below overrides lib-ts-analysis.R)
```{r}
every_hour_labels = format(lubridate::parse_date_time(hms(hours = 0:23), c('HMS', 'HM')), '%H:%M')
labels = c("Quercus", "Declining Pines", "Non-Declining Pines")

n.per.class <- db %>% group_by(class, year) %>% summarise (n = n_distinct(series)) %>% pivot_wider(names_from = class, values_from = n)

c_labels_2022 = c(glue("{labels[1]} ({n.per.class[1,'Quercus']})"), glue("{labels[2]} ({n.per.class[1,'D']})"), glue("{labels[3]} ({n.per.class[1,'ND']})"))
c_labels_2023 = c(glue("{labels[1]} ({n.per.class[2,'Quercus']})"), glue("{labels[2]} ({n.per.class[2,'D']})"), glue("{labels[3]} ({n.per.class[2,'ND']})"))
c_values = c(Quercus = "purple", D = "darkred", ND = "darkorange")
```

Making bigger text:
```{r}
theme_set( theme_bw() +
          theme(
            plot.title = element_text(size = 20, hjust = 0.5, face = "bold"),
            text=element_text(size=20, color = "gray9"),
            legend.title=element_text(size=16, color = "gray9"),
            legend.text=element_text(size=16, color = "gray9"),
            axis.title=element_text(size=20),
            axis.text=element_text(size=20, color = "gray9")
          ))
```

# TWD
```{r}
db.twd <- db %>% summarise( mean_twd = mean(twd, na.rm = T), twd_se = sd(twd, na.rm = T) / sqrt(n()), .by = c(date, class), n = as.integer(n() / 96)) %>% mutate(doy = yday(date))
db.twd2022 <- db.twd %>% filter(year(date) == 2022)
db.twd2023 <- db.twd %>% filter(year(date) == 2023)
db.twd
```

```{r}
db.twd.agg.yearly <- db %>% mutate(year = year(date)) %>% summarise( mean_twd = mean(twd, na.rm = T), n = n_distinct(series), twd_se = sd(twd, na.rm = T) / sqrt(n), .by = c(series, class, year))
boxplot(mean_twd ~ class + year, data = db.twd.agg.yearly)
```

```{r}
db.twd.agg.2022 <- db.twd.agg.yearly %>% filter(year == 2022)
db.twd.agg.2023 <- db.twd.agg.yearly %>% filter(year == 2023)
boxplot(mean_twd ~ class, data = db.twd.agg.2022, title = '2022')
boxplot(mean_twd ~ class, data = db.twd.agg.2023)
```

```{r}
aov.model <- aov(mean_twd ~ class + year, data = db.twd.agg.yearly)
summary(aov.model)
aov.model <- aov(mean_twd ~ class, data = db.twd.agg.yearly)
summary(aov.model)
aov.model <- aov(mean_twd ~ class, data = db.twd.agg.2022)
summary(aov.model)
aov.model <- aov(mean_twd ~ class, data = db.twd.agg.2023)
summary(aov.model)
```

Boxplot:
```{r}
plot.twd2022 <- db.twd2022 %>%
  ggplot () +
  geom_boxplot(aes (x = class, y = mean_twd, color = class)) +
  scale_color_manual(labels = c_labels_2022, values=c_values) +
  scale_x_discrete(labels=c("Q" = labels[1], "D" = labels[2], "ND" = labels[3])) +
  labs(title = "2022", x = "Class", y = expression(paste("Tree Water Deficit (", mu, "m)")))

plot.twd2023 <- db.twd2023 %>%
  ggplot () +
  geom_boxplot(aes (x = class, y = mean_twd, color = class)) +
   scale_color_manual(labels = c_labels_2023, values=c_values) +
  labs(title = "2023", x = "Class", y = expression(paste("Tree Water Deficit (", mu, "m)"))) +
  scale_x_discrete(labels=c("Q" = labels[1], "D" = labels[2],
                              "ND" = labels[3]))

combined <- (plot.twd2022 | plot.twd2023) + plot_layout(guides = "collect", axis_titles = "collect") + plot_annotation(title = glue('Boxplot of TWD dispersion - {PLACE}')) & theme(legend.position = "bottom")
combined
```

```{r}
aov.model <- aov(mean_twd ~ class, data = db.twd2022)
summary(aov.model)
TukeyHSD(aov.model)
```

```{r}
aov.model <- aov(mean_twd ~ class, data = db.twd2023)
summary(aov.model)
TukeyHSD(aov.model)
```

```{r}
db.twd <- db.twd %>% mutate( year = year(date))
str(db.twd)
aov.model <- aov(mean_twd ~ class + year, data = db.twd)
summary(aov.model)
```

```{r}
if (SAVE) ggsave(glue('output/boxplot-TWD-byClass-{PLACE}.png'), combined, width = 14, height = 6)
```

Continuous TWD evolution along with VWC:

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


```{r}
plot.twd2023 <- db.twd %>%
  ggplot () +
  geom_line(aes (x = date, y = mean_twd, color = class), show.legend = T) +
  geom_line(aes (x = date, y = (mean_twd + twd_se), color = class), alpha = 0.5, linetype = "dashed", show.legend = F) +
  geom_line(aes (x = date, y = (mean_twd - twd_se), color = class), alpha = 0.5, linetype = "dashed", show.legend = F) +
  scale_color_manual(labels = c_labels_2022, values=c_values) +
  labs(title = "2023", x = "Date", y = expression(paste("Tree Water Deficit (", mu, "m)"))) +
  scale_x_date(limits = as.Date(c("2023-01-01", "2023-12-31")), breaks = seq(as.Date("2023-01-01"), as.Date("2023-12-31"), by = "1 month"), date_labels = "%b-%y") +
  scale_y_continuous(sec.axis = sec_axis(trans = ~ . /10, name = "Volumetric Water Content (%)")) +
   geom_line(data = vwc, aes(x = date, y = max*1000, linetype = "Volumetric Water Content (%)"), col = "blue", alpha = 0.8, show.legend	= c(colour = FALSE, linetype = T)) +
  scale_linetype_manual(NULL, values = 4)

plot.twd2022 <- db.twd %>%
  ggplot () +
  geom_line(aes (x = date, y = mean_twd, color = class), show.legend = T) +
  geom_line(aes (x = date, y = (mean_twd + twd_se), color = class), alpha = 0.5, linetype = "dashed", show.legend = F) +
  geom_line(aes (x = date, y = (mean_twd - twd_se), color = class), alpha = 0.5, linetype = "dashed", show.legend = F) +
  scale_color_manual(labels = c_labels_2023, values=c_values) +
  labs(title = "2022", x = "Date", y = expression(paste("Tree Water Deficit (", mu, "m)"))) +
  scale_x_date(limits = as.Date(c("2022-01-01", "2022-12-31")), breaks = seq(as.Date("2022-01-01"), as.Date("2022-12-31"), by = "1 month"), date_labels = "%b-%y") +
  scale_y_continuous(sec.axis = sec_axis(trans = ~ . /10, name = "Volumetric Water Content (%)")) +
   geom_line(data = vwc, aes(x = date, y = max*1000, linetype = "Volumetric Water Content (%)"), col = "blue", alpha = 0.8, show.legend	= c(colour = FALSE, linetype = T)) +
  scale_linetype_manual(NULL, values = 4)

combined <- plot.twd2022 / plot.twd2023 + 
  plot_layout(guides = "collect", axis_titles = "collect") + 
  plot_annotation(title = glue('TWD evolution by classes along with Soil Moisture - {PLACE}')) &
  theme(legend.position = "bottom")
combined
```
```{r}
if (SAVE) ggsave(glue('output/evolution-TWD-byClass-{PLACE}.png'), combined, width = 14, height = 6)
```

Same but now along with temp min and temp max:

```{r}
temp <- db.env %>% mutate(date = date(db.env$ts)) %>% summarise(tmax = max(surface.temp), tmin = min(surface.temp), .by = date)
temp
```

```{r warning=FALSE}
plot.twd2023 <- db.twd %>%
  ggplot () +
  geom_line(aes (x = date, y = mean_twd, color = class), show.legend = T) +
  geom_line(aes (x = date, y = (mean_twd + twd_se), color = class), alpha = 0.5, linetype = "dashed", show.legend = F) +
  geom_line(aes (x = date, y = (mean_twd - twd_se), color = class), alpha = 0.5, linetype = "dashed", show.legend = F) +
  scale_color_manual(labels = c_labels_2022, values=c_values) +
  labs(title = "2023", x = "Date", y = expression(paste("Tree Water Deficit (", mu, "m)")) , alpha = "Climate") +
  scale_x_date(limits = as.Date(c("2023-01-01", "2023-12-31")), breaks = seq(as.Date("2023-01-01"), as.Date("2023-12-31"), by = "1 month"), date_labels = "%b-%y") +
  scale_y_continuous(sec.axis = sec_axis(trans = ~ . /20, name = "Temperature (ºC)")) +
  geom_linerange(data = temp, aes(x = date, ymin = tmin*20, ymax = tmax*20, alpha = "temp"), col = "firebrick1") +
  scale_alpha_manual(labels = c("Range from minimum to maximum temperature (ºC)"), values = c(temp = 0.4)) # to set legend text

plot.twd2022 <- db.twd %>%
  ggplot () +
  geom_line(aes (x = date, y = mean_twd, color = class), show.legend = T) +
  geom_line(aes (x = date, y = (mean_twd + twd_se), color = class), alpha = 0.5, linetype = "dashed", show.legend = F) +
  geom_line(aes (x = date, y = (mean_twd - twd_se), color = class), alpha = 0.5, linetype = "dashed", show.legend = F) +
  scale_color_manual(labels = c_labels_2023, values=c_values) +
  labs(title = "2022", x = "Date", y = expression(paste("Tree Water Deficit (", mu, "m)")), alpha = "Climate") +
  scale_x_date(limits = as.Date(c("2022-01-01", "2022-12-31")), breaks = seq(as.Date("2022-01-01"), as.Date("2022-12-31"), by = "1 month"), date_labels = "%b-%y") +
  scale_y_continuous(sec.axis = sec_axis(trans = ~ . /20, name = "Temperature (ºC)")) +
  geom_linerange(data = temp, aes(x = date, ymin = tmin*20, ymax = tmax*20, alpha = "temp"), col = "firebrick1") +
  scale_alpha_manual(labels = c("Range from minimum to maximum temperature (ºC)"), values = c(temp = 0.4)) # to set legend text

combined <- plot.twd2022 / plot.twd2023 + 
  plot_layout(guides = "collect", axis_titles = "collect") + 
  plot_annotation(title = glue('TWD evolution by classes along with Daily Temperature - {PLACE}')) &
  theme(legend.position = "bottom")
combined
```

```{r}
if (SAVE) ggsave(glue('output/evolution-Temp-byClass-{PLACE}.png'), combined, width = 14, height = 6)
```


# GRO
```{r}
db.gro <- db %>% summarise( gro = mean(gro_yr, na.rm = T), gro_se = sd(gro_yr, na.rm = T) / sqrt(n()), .by = c(date, class)) %>% mutate(doy = yday(date))
db.gro2022 <- db.gro %>% filter(year(date) == 2022)
db.gro2023 <- db.gro %>% filter(year(date) == 2023)
db.gro
```

```{r}
plotgro2022 <- db.gro2022 %>% 
  ggplot() +
  ggtitle('2022') +
  geom_line(aes (x = doy, y = gro, color = class)) +
  geom_line(aes (x = doy, y = (gro + gro_se), color = class), alpha = 0.5, linetype = "dashed") +
  geom_line(aes (x = doy, y = (gro - gro_se), color = class), alpha = 0.5, linetype = "dashed") +
  labs(x = "Day of the year", y = expression(paste("Acumulated growth (",mu, "m)")) ) +
  scale_color_manual(labels = c_labels_2022, values=c_values) +
  scale_x_continuous(breaks = seq(0, 360, by = 30), limits = c(0, 366))

plotgro2023 <- db.gro2023 %>% ggplot () +
  ggtitle('2023') +
  geom_line(aes (x = doy, y = gro, color = class)) +
  geom_line(aes (x = doy, y = (gro + gro_se), color = class), alpha = 0.5, linetype = "dashed") +
  geom_line(aes (x = doy, y = (gro - gro_se), color = class), alpha = 0.5, linetype = "dashed") +
  labs(x = "Day of the year", y = expression(paste("Acumulated growth (",mu, "m)")) ) +
  scale_color_manual(labels = c_labels_2023, values=c_values) +
  scale_x_continuous(breaks = seq(0, 360, by = 30), limits = c(0, 366))

combined <- plotgro2022 / plotgro2023 + plot_layout(guides = "collect", axis_titles = "collect") + plot_annotation(title = glue('Acumulated Growth - {PLACE}')) & theme(legend.position = "bottom")
combined
```

```{r}
if (SAVE) ggsave(glue('output/GRO-byClass-{PLACE}.png'), combined, width = 14, height = 6)
```

# Growing rate
Using Growing rate instead of GRO from treenetproc.

First, we will use the trend only to get more accurate values.

For that, we need to fill the NA values and then use a time-series decomposition function.
```{r}
db.selected <- db %>% filter(db$series %in% SELECTED_DENDROS)

statsNA(db.selected$value)

db.selected$value <- db.selected$value %>% na_interpolation(option = "spline")
```

Then we can extract the trend to have more accurate results
```{r}
db.trends <- decompose_ts_stl(db.selected, SELECTED_DENDROS)
db.trends <- db.trends %>% select(ts, trend, series) %>% rename(value = trend) %>%
  mutate (class = case_when(
                    series %in% Qi ~ factor("Quercus"),
                    series %in% P_D ~ factor("D"),
                    series %in% P_ND ~ factor("ND"),
                    .default = NA
                    )
)
str(db.trends)
```

Now we calculate the Growing rate for every dendrometer and we compute the aggregated mean per classes:
```{r}
db.gr <- data.frame()
for (dendro.no in SELECTED_DENDROS) {
  # from database of trends, filter data of the time series of this dendro:
  dat = db.trends[db.trends$series == dendro.no,] %>% select(ts, value, class)
  aux <- data.frame(
     series = as.factor(dendro.no),
     class = dat$class[1],
     calc.growth.rate(dat) # calculates Growing rate for the given data
  ) 
  db.gr <- rbind.data.frame(db.gr, aux) # append in a output / result dataframe.
}

db.gr
```

Let's restrict the output to the study period:
```{r}
db.gr2022 <- db.gr %>% filter(date >= ts_start2022 & date <= ts_end2022)
db.gr2023 <- db.gr %>% filter(date >= ts_start2023 & date <= ts_end2023)
```

Now we calculate means for each class:
```{r}
db.gr2022 <- db.gr2022 %>% summarise( mean_gr = mean(rate, na.rm = T), gr_se = sd(rate, na.rm = T) / sqrt(n()), .by = c(date, class)) %>% mutate(doy = yday(date))
db.gr2023 <- db.gr2023 %>% summarise( mean_gr = mean(rate, na.rm = T), gr_se = sd(rate, na.rm = T) / sqrt(n()), .by = c(date, class)) %>% mutate(doy = yday(date))
```

```{r}
unique(db.gr2022$doy)
```


Finally, we plot it:

```{r}
plotgr2022 <- db.gr2022 %>% 
  ggplot() +
  ggtitle('2022') +
  geom_line(aes (x = doy, y = mean_gr, color = class)) +
  geom_line(aes (x = doy, y = (mean_gr + gr_se), color = class), alpha = 0.5, linetype = "dashed") +
  geom_line(aes (x = doy, y = (mean_gr - gr_se), color = class), alpha = 0.5, linetype = "dashed") +
  labs(x = "Day of the year", y = bquote("Growing rate (" ~ mu * "m" ~ "·" ~ day ^-1 ~ ")" ) ) +
  scale_color_manual(name="Tree class", labels = c_labels, values = c_class_values, 
                   guide = guide_legend(order = 1)) +
  scale_x_continuous(breaks = seq(90, 360, by = 30), limits = c(91, 334), expand = expansion( mult = 0.01))

plotgr2023 <- db.gr2023 %>%
  ggplot () +
  ggtitle('2023') +
  geom_line(aes (x = doy, y = mean_gr, color = class)) +
  geom_line(aes (x = doy, y = (mean_gr + gr_se), color = class), alpha = 0.5, linetype = "dashed") +
  geom_line(aes (x = doy, y = (mean_gr - gr_se), color = class), alpha = 0.5, linetype = "dashed") +
  labs(x = "Day of the year",  y = bquote("Growing rate (" ~ mu * "m" ~ "·" ~ day ^-1 ~ ")" ) ) +
   scale_color_manual(name="Tree class", labels = c_labels, values = c_class_values, 
                   guide = guide_legend(order = 1)) +
  scale_x_continuous(breaks = seq(90, 360, by = 30), limits = c(91, 334), expand = expansion( mult = 0.01))

combined <- plotgr2022 * plotgr2023 + plot_layout(guides = "collect", axis_titles = "collect") & theme(legend.position = "bottom")
combined
```

```{r}
if (SAVE) ggsave(glue('output/GrowthRate-byClass-{PLACE}.png'), combined, width = 4200, height = 3000, units = "px")
```

And smoothed to get better the higher level trend with a smooth line:
```{r}
# plotgr2022 <- db.gr2022 %>% 
#   ggplot() +
#   ggtitle('2022') +
#   geom_smooth(span = 0.1, aes (x = doy, y = mean_gr, color = class), se=T) +
#   labs(x = "Day of the year", y = bquote("Growing rate (" ~ um ~ "·" ~ day ^-1 ~ ")" ) ) +
#   scale_color_manual(labels = c_labels, values=c(Quercus = "green", D = "darkred", ND = "darkorange")) +
#   scale_x_continuous(breaks = seq(0, 360, by = 30), limits = c(0, 366))
# 
# plotgr2023 <- db.gr2023 %>%
#   ggplot () +
#   ggtitle('2023') +
#   geom_smooth(span = 0.1, aes (x = doy, y = mean_gr, color = class), se=T) +
#   labs(x = "Day of the year",  y = bquote("Growing rate (" ~ um ~ "·" ~ day ^-1 ~ ")" ) ) +
#   scale_color_manual(labels = c_labels, values=c(Quercus = "green", D = "darkred", ND = "darkorange")) +
#   scale_x_continuous(breaks = seq(0, 360, by = 30), limits = c(0, 366))
# 
# combined <- plotgr2022 / plotgr2023 + plot_layout(guides = "collect", axis_titles = "collect") + plot_annotation(title = glue('Growing rate per day - {PLACE}')) & theme(legend.position = "bottom")
# combined
```