plot-proc-dendro-Corbalan.Rmd

```{r}
source("lib-dendro.R")

library(ggplot2)
library(dplyr)
library(lubridate)

### DEFINE GLOBAL VARS ###
PATH = '/home/akronix/workspace/dendro';
setwd(PATH)

DATA_DIR = 'processed/Corbalan-processed'
ENV_DIR = 'processed/Corbalan-env-processed'

SELECTED_TMS <- "94231943"

### IMPORT DENDRO DATA ###

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

# Set initial and final date for analysis
# ts_start<-"2022-04-01 09:00:00" # 2 days after installation
# ts_end<-"2023-09-27 08:00:00" # last timestamp of downloaded data
# db <- reset.initial.values(db, ts_start, ts_end)

### CLEAN & PREPARE DATA ###

# INSPECT DATA
str(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('output/all-dendrometers-Corbalan.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('output/normalized-all-dendrometers-Corbalan.png', width = 15, height = 10)
```


# Mean + strip region of SE
```{r}
dbagg = db %>% group_by(ts) %>% dplyr::summarise(mean = mean(value, na.rm = T), sd = sd(value, na.rm = T) )

#library(Rmisc) # for summarySE()
#dbagg <- summarySE(db, measurevar = "value", groupvars = c("ts"), na.rm = TRUE)
head(dbagg)
tail(dbagg)
```

```{r}
plot_cat_full_year<-
  ggplot(data = dbagg,  aes(x=ts, y=mean)) +
  ggtitle("Mean of data for Corbalan dendrometers with ± standard deviation strips") +
  geom_ribbon(aes(ymin=mean-sd, ymax=mean+sd), fill='lightgreen', alpha=0.3, show.legend	= FALSE, linetype = 0) +
  geom_line( aes (linetype = ""), col='darkgreen', show.legend	= F) +
  labs(x=expression(''),
       y=expression(Delta*" D (um)"))+
  theme_bw() +
  geom_hline(yintercept=0,lty=2,linewidth=0.2)+
  #facet_grid(class~.,scales = "free_y")+
  scale_x_datetime(date_breaks = "1 month", date_labels="%b %Y") +
  theme(axis.text.x = element_text(angle = 30, hjust=1))
plot_cat_full_year
ggsave('output/mean-Corbalan.png', width = 15, height = 10)
```


# By species

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

Qi = TreeList %>% filter(class == "Quercus") %>% pull(series)
P_D = TreeList %>% filter(class == "D") %>% pull(series)
P_ND = TreeList %>% filter(class == "ND") %>% pull(series)

db.Qi = db %>% filter (series %in% Qi)
db.D = db %>% filter (series %in% P_D)
db.ND = db %>% filter (series %in% P_ND)
```


```{r}
library(ggnewscale)

# Define color scales for each class
green_colors <- c("lightgreen", "green", "darkgreen", "darkolivegreen", "olivedrab", "springgreen")
blue_colors <- c("blue", "cyan", "aquamarine", "cadetblue", "dodgerblue", "darkblue")
orange_colors <- c("orange", "orangered", "tomato", "darkorange", "coral", "tan")

ggplot() +
   ggtitle("Dendrometers painted by color scale according to its class") +
  geom_line(data = db.D, aes(x = ts, y = value, col = series)) +
  scale_color_manual(values = orange_colors, name = "Declining trees") +
  new_scale_color() +
  geom_line(data = db.ND, aes(x = ts, y = value, col = series), ) +
  scale_color_manual(values = blue_colors, name = "Non-Declining trees") +
  new_scale_color() +
  geom_line(data = db.Qi, aes(x = ts, y = value, col = series)) +
  scale_color_manual(values = green_colors, name = "Quercus Ilex") +
  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('output/all-dendrometers-by-class-Corbalan.png', width = 15, height = 10)
```



```{r}
# Dataframe with all data and a new column "class"
# db.by.class = db %>% select(series, ts, value) %>% mutate(class =  case_when(
#                                       series %in% Qi ~ "Quercus",
#                                       series %in% P_D ~ "Declining",
#                                       series %in% P_ND ~ "Non-Declining"
#                                       )
# )
```

One plot of dendros by class:
```{r}
plot_multiple_dendro(db.Qi, "Quercus Ilex - Corbalan", y = db.Qi$value)
ggsave('output/Qi-Corbalan.png', width = 15, height = 10)
```

```{r}
plot_multiple_dendro(db.ND, "Non-Declining - Corbalan", y = db.ND$value)
ggsave('output/ND-Corbalan.png', width = 15, height = 10)
```

```{r}
plot_multiple_dendro(db.D, "Declining - Corbalan", y = db.D$value)
ggsave('output/D-Corbalan.png', width = 15, height = 10)
```


```{r}
db.Qi.agg <- db.Qi %>% group_by(ts) %>% dplyr::summarise(mean = mean(value, na.rm = T), sd = sd(value, na.rm = T) )
head(db.Qi.agg)
tail(db.Qi.agg)
```


```{r}
#library(mdthemes)
plot_mean_sd <- function (data, species, output_fn = FALSE) {
  plot <-
    ggplot(data = data,  aes(x=ts, y=mean)) +
    ggtitle(paste0("Mean of data for ", species, " with ± standard deviation strips")) +
    geom_ribbon(aes(ymin=mean-sd, ymax=mean+sd), fill='lightgreen', alpha=0.3, show.legend	= FALSE, linetype = 0) +
    geom_line( aes (linetype = species), col='darkgreen', show.legend	= T) +
    labs(x=expression(''),
         y=expression(Delta*" D (um)")) +
    theme_bw() +
    #md_theme_bw()+
    geom_hline(yintercept=0,lty=2,linewidth=0.2) +
    #facet_grid(class~.,scales = "free_y")+
    scale_x_datetime(date_breaks = "1 month", date_labels="%b %Y") +
    theme(axis.text.x = element_text(angle = 30, hjust=1))
        
  plot(plot)
  if (output_fn != F){
    ggsave(output_fn, width = 15, height = 10)  
  }
}
```


```{r}
plot_mean_sd(db.Qi.agg, "Quercus Ilex", output_fn = "output/mean-Qi-Corbalan.png")
```


```{r}
db.D.agg <- db.D %>% group_by(ts) %>% dplyr::summarise(mean = mean(value, na.rm = T), sd = sd(value, na.rm = T) )
head(db.D.agg)
tail(db.D.agg)
```


```{r}
plot_mean_sd(db.D.agg, "Declining Pinus", output_fn = "output/mean-declining-Corbalan.png" )
```

```{r}
db.ND.agg <- db.ND %>% group_by(ts) %>% dplyr::summarise(mean = mean(value, na.rm = T), sd = sd(value, na.rm = T) )
head(db.ND.agg)
tail(db.ND.agg)
```


```{r}
plot_mean_sd(db.ND.agg, "Non-Declining Pinus", output_fn = "output/mean-ND-Corbalan.png" )
```

Join all aggregated df in one
```{r}
db.Qi.agg$class = factor("Quercus Ilex")
db.D.agg$class = factor("Declining Pinus")
db.ND.agg$class = factor("Non-Declining Pinus")
db.by.class = rbind.data.frame(db.Qi.agg, db.D.agg, db.ND.agg)
str(db.by.class)
```

Plot mean by class in one plot
```{r}
plot_multiple_dendro_by_class <-
  ggplot(data = db.by.class, mapping = aes(x=ts, y=mean, col = class))+
  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("Means by class") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 30, hjust=1))
plot_multiple_dendro_by_class
ggsave('output/mean-byclass-Corbalan.png', width = 15, height = 10)
```

```{r}
# library(plotly)
# plot_ly(
#   data = db.by.class,
#    x = ~ts,
#    y = ~mean,
#    type = "scatter",
#    mode = "lines",
#    color = ~class)
# +   colors = c("blue")
```


# Environmental data from TMS sensors

```{r}
# importing processed environmental data #
list_files <- list.files(file.path(PATH,ENV_DIR), pattern="*.csv$", full.names=TRUE)
db.env<-read.all.env.processed(list_files)
str(db.env)
```

Now calculate the mean for vwc and temp:
```{r}
db.env.agg = db.env %>% group_by(ts) %>% dplyr::summarise(temp = mean(surface.temp, na.rm = T), vwc = mean(vwc, na.rm = T))
db.env.agg
```

Let's keep environmental data from only one dendro whose data is all valid:
```{r}
db.env <- db.env[db.env$series == SELECTED_TMS,]
```

And we plot it:

```{r}
plot_temp_and_humidity <-
  ggplot(data = db.env, aes(x=ts, y=surface.temp)) +
    ggtitle("Temperature and humidity for Corbalan (2022-2023)") +
    geom_line(aes(colour = surface.temp)) +
    theme_bw() +
    scale_colour_gradient(low = "light blue", high = "red") +
    scale_x_datetime(date_breaks = "1 month", date_labels="%b %Y") +
    theme(axis.text.x = element_text(angle = 30, hjust=1)) +
    scale_y_continuous("Temperature (ºC)", breaks=seq(0,50,10), sec.axis = dup_axis(name = "Volumetric Water Content (%)" ))+
    geom_line(aes(x = ts, y = vwc * 100,  linetype = "Volumetric Water Content (%)"), col="blue", show.legend	= TRUE) +
    scale_linetype_manual(NULL, values = 1) +
    labs(x=expression(''), colour = "Temperature (ºC)")
plot_temp_and_humidity
ggsave('output/VWCytempAgg-Corbalan.png', width = 15, height = 10)
```

```{r}
# library(plotly)
# plot_ly(
#   data = db.env,
#   x = ~ts,
#   y = ~vwc*100,
#   type = "scatter",
#   mode = "lines",
#   # color = ~state,
#   colors = c("blue")
#   )
```