Seasonality-aggregated.Rmd

---
title: "Seasonality-analysis"
author: "Abel Serrano Juste"
date: "`r Sys.Date()`"
output: html_document
---

# INITIALIZATION
SET GLOBAL VARS
```{r}
PATH = dirname(rstudioapi::getSourceEditorContext()$path)
print(PATH)
setwd(PATH)

PLACE = 'Penaflor'
```

Load common libraries and data set
```{r}
source('init-analysis.R')
```

Import other libraries
```{r}
source('correlations.R')
library('correlation')
library('biwavelet')
```


# Calculate mean Seasonality
Mean seasonality for all trees:
```{r}
seasonalities <- calculate_stl_seasonalities(db, unique(db$series))
summary(seasonalities)

seasonalities.agg <- seasonalities %>% group_by(ts) %>% summarise(mean = mean(value), sd = sd(value))
summary(seasonalities.agg)
```

```{r}
ggplot( data = seasonalities.agg, mapping = aes(x=ts, y = mean)) + 
  geom_line() + 
  ggtitle("Mean of seasonalities for all trees")
# if (SAVE) ggsave(glue('output/{PLACE}-mean-seasonalities-all.png'), width = 15, height = 10)
```

# Calculate Seasonality by class

Calculate mean of the seasonality for Qi, PD and P_ND
```{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)
```

```{r}
every_hour_labels = format(lubridate::parse_date_time(hms(hours = 0:23), c('HMS', 'HM')), '%H:%M')
```

## Quercus
Mean seasonality for Quercus Ilex:
```{r}
seasonalities.qi <- calculate_stl_seasonalities(db, Qi)
summary(seasonalities.qi)

seasonalities.qi.agg <- seasonalities.qi %>% group_by(ts) %>% summarise(mean = mean(value), sd = sd(value))
summary(seasonalities.qi.agg)
```

```{r}
ggplot( data = seasonalities.qi, mapping = aes(x=ts, y = value, color = series)) + 
  geom_line(alpha = 0.5) +
  ggtitle(glue('Quercus Ilex seasonalities by dendrometer - {PLACE} '))
# if (SAVE) ggsave(glue('output/{PLACE}-seasonality-Qi.png'), width = 15, height = 10)
```

```{r}
ggplot( data = seasonalities.qi.agg, mapping = aes(x=ts, y = mean)) + 
  geom_line() + 
  ggtitle("Mean of seasonalities for Quercus Ilex")
# if (SAVE) ggsave(glue('output/{PLACE}-mean-seasonality-Qi.png'), width = 15, height = 10)
```

### Aggregating in one day, by removing time-series dimension

removing the time-series dimension, mean of one day for all period:
```{r}
seasons.qi.allperiod = seasonalities.qi.agg %>% mutate(timeOfDay = as_hms(ts)) %>% summarise(meanSeasonalityTime = mean(mean), SE_SeasonalityTime = sd(sd)/sqrt(n()), .by = timeOfDay)
seasons.qi.allperiod
plot_day_seasonality(seasons.qi.allperiod, "Quercus Ilex", PLACE, ALL_PERIOD_ST)
# if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-allperiod-seasonalities-Qi.png'), width = 15, height = 10)
```

removing the time-series dimension, mean of one day for June-July:
```{r}
seasons.qi.jun_jul = seasonalities.qi.agg %>% filter( (ts >= as.Date("2023-06-01") ) & (ts < as.Date("2023-08-01")) ) %>%  mutate(timeOfDay = as_hms(ts)) %>% summarise(meanSeasonalityTime = mean(mean), SE_SeasonalityTime = sd(sd)/sqrt(n()), .by = timeOfDay)
seasons.qi.jun_jul

plot_day_seasonality(seasons.qi.jun_jul, "Quercus Ilex", PLACE, "June to July 2023")
# if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-june-july-seasonalities-Qi.png'), width = 15, height = 10)
```


Now, the two plots above together in one plot:

```{r}
seasons.qi.jun_jul$period <- "June-July"
seasons.qi.allperiod$period <- "All"
seasons.periods.joined <- rbind.data.frame(seasons.qi.allperiod, seasons.qi.jun_jul)
```


```{r}
plot_day_seasonalities_periods(seasons.periods.joined, "Quercus Ilex", PLACE)
# if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-bothperiods-seasonalities-Qi.png'), width = 15, height = 10)
```


```{r}
# Another way to do it but using two different df sources
# ggplot() +
#   ggtitle("Aggregated mean in one day for all period And for period june-july") +
#   geom_line(data = seasons.qi.allperiod, mapping = aes(x=timeOfDay, y = meanSeasonalityTime), col = "blue") + 
#   geom_line(data = seasons.qi.allperiod, aes (x = timeOfDay, y = (meanSeasonalityTime + SE_SeasonalityTime)), col = "red", alpha = 0.5, linetype = "dashed", show.legend = F) +
#   geom_line(data = seasons.qi.allperiod, aes (x = timeOfDay, y = (meanSeasonalityTime - SE_SeasonalityTime)), col = "red", alpha = 0.5, linetype = "dashed", show.legend = F) +
#   geom_line(data = seasons.qi.jun_jul, mapping = aes(x = timeOfDay, y = meanSeasonalityTime), col= "purple") +
#   geom_line(data = seasons.qi.jun_jul, aes (x = timeOfDay, y = (meanSeasonalityTime + SE_SeasonalityTime)), col = "red", alpha = 0.5, linetype = "dashed", show.legend = F) +
#   geom_line(data = seasons.qi.jun_jul, aes (x = timeOfDay, y = (meanSeasonalityTime - SE_SeasonalityTime)), col = "red", alpha = 0.5, linetype = "dashed", show.legend = F) +
#   theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
#   scale_x_time(breaks = seq(0, 85500, by = 3600), labels = every_hour_labels) +
#   labs(x = "Hour of the day", y = "Micrometers of Daily seasonality (um)")
```


## Non-Declining Pines
Mean seasonality for Non-Declining Pines
```{r}
seasonalities.P_ND <- calculate_stl_seasonalities(db, P_ND)
summary(seasonalities.P_ND)

seasonalities.P_ND.agg <- seasonalities.P_ND %>% group_by(ts) %>% summarise(mean = mean(value), sd = sd(value))
summary(seasonalities.P_ND.agg)
```

```{r}
ggplot( data = seasonalities.P_ND, mapping = aes(x=ts, y = value, color = series)) + 
  geom_line(alpha = 0.5) +
  ggtitle(glue('Non-Declining Pinus seasonalities by dendrometer - {PLACE} '))
# if (SAVE) ggsave(glue('output/{PLACE}-seasonality-P_ND.png'), width = 15, height = 10)
```

```{r}
ggplot( data = seasonalities.P_ND.agg, mapping = aes(x=ts, y = mean)) + 
  geom_line() + 
  ggtitle("Mean of seasonalities for Non-Declining Pinus")
# if (SAVE) ggsave(glue('output/{PLACE}-mean-seasonality-P_ND.png'), width = 15, height = 10)
```


### Non-Declining Pines: Aggregating in one day, by removing time-series dimension

removing the time-series dimension, mean of one day for all period:
```{r}
seasons.P_ND.allperiod = seasonalities.P_ND.agg %>% mutate(timeOfDay = as_hms(ts)) %>% summarise(meanSeasonalityTime = mean(mean), SE_SeasonalityTime = sd(sd)/sqrt(n()), .by = timeOfDay)
seasons.P_ND.allperiod
plot_day_seasonality(seasons.P_ND.allperiod, "Non-Declining Pines", PLACE, ALL_PERIOD_ST)
# if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-allperiod-seasonalities-P_ND.png'), width = 15, height = 10)
```

removing the time-series dimension, mean of one day for June-July:
```{r}
seasons.P_ND.jun_jul = seasonalities.P_ND.agg %>% filter( (ts >= as.Date("2023-06-01") ) & (ts < as.Date("2023-08-01")) ) %>%  mutate(timeOfDay = as_hms(ts)) %>% summarise(meanSeasonalityTime = mean(mean), SE_SeasonalityTime = sd(sd)/sqrt(n()), .by = timeOfDay)
seasons.P_ND.jun_jul

plot_day_seasonality(seasons.P_ND.jun_jul, "Non-Declining Pines", PLACE, "June to July 2023")
# if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-june-july-seasonalities-P_ND.png'), width = 15, height = 10)
```


Now, the two plots above together in one plot:

```{r}
seasons.P_ND.jun_jul$period <- "June-July"
seasons.P_ND.allperiod$period <- "All"
seasons.periods.joined <- rbind.data.frame(seasons.P_ND.allperiod, seasons.P_ND.jun_jul)
```


```{r}
plot_day_seasonalities_periods(seasons.periods.joined, "Non-Declining Pines", PLACE)
# if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-bothperiods-seasonalities-P_ND.png'), width = 15, height = 10)
```


## Declining Pines
Mean seasonality for Declining Pines
```{r}
seasonalities.P_D <- calculate_stl_seasonalities(db, P_D)
summary(seasonalities.P_D)

seasonalities.P_D.agg <- seasonalities.P_D %>% group_by(ts) %>% summarise(mean = mean(value), sd = sd(value))
summary(seasonalities.P_D.agg)
```

```{r}
ggplot( data = seasonalities.P_D, mapping = aes(x=ts, y = value, color = series)) + 
  geom_line(alpha = 0.5) +
  ggtitle(glue('Declining Pinus seasonalities by dendrometer - {PLACE} '))
# if (SAVE) ggsave(glue('output/{PLACE}-seasonality-P_D.png'), width = 15, height = 10)
```

```{r}
ggplot( data = seasonalities.P_D.agg, mapping = aes(x=ts, y = mean)) + 
  geom_line() + 
  ggtitle("Mean of seasonalities for Declining Pinus")
# if (SAVE) ggsave(glue('output/{PLACE}-mean-seasonality-P_D.png'), width = 15, height = 10)
```


### Declining Pines: Aggregating in one day, by removing time-series dimension

removing the time-series dimension, mean of one day for all period:
```{r}
seasons.P_D.allperiod = seasonalities.P_D.agg %>% mutate(timeOfDay = as_hms(ts)) %>% summarise(meanSeasonalityTime = mean(mean), SE_SeasonalityTime = sd(sd)/sqrt(n()), .by = timeOfDay)
seasons.P_D.allperiod
plot_day_seasonality(seasons.P_D.allperiod, "Declining Pines", PLACE, ALL_PERIOD_ST)
# if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-allperiod-seasonalities-P_D.png'), width = 15, height = 10)
```

removing the time-series dimension, mean of one day for June-July:
```{r}
seasons.P_D.jun_jul = seasonalities.P_D.agg %>% filter( (ts >= as.Date("2023-06-01") ) & (ts < as.Date("2023-08-01")) ) %>%  mutate(timeOfDay = as_hms(ts)) %>% summarise(meanSeasonalityTime = mean(mean), SE_SeasonalityTime = sd(sd)/sqrt(n()), .by = timeOfDay)
seasons.P_D.jun_jul

plot_day_seasonality(seasons.P_D.jun_jul, "Declining Pines", PLACE, "June to July 2023")
# if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-june-july-seasonalities-P_D.png'), width = 15, height = 10)
```

Now, the two plots above together in one plot:

```{r}
seasons.P_D.jun_jul$period <- "June-July"
seasons.P_D.allperiod$period <- "All"
seasons.periods.joined <- rbind.data.frame(seasons.P_D.allperiod, seasons.P_D.jun_jul)
```


```{r}
plot_day_seasonalities_periods(seasons.periods.joined, "Declining Pines", PLACE)
# if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-bothperiods-seasonalities-P_D.png'), width = 15, height = 10)
```

### Plot aggregated data of each class in one plot

All study period:
```{r}
seasons.allclases.allperiod <- data.frame()
seasons.P_D.allperiod$class = "D"
seasons.P_ND.allperiod$class = "ND"
seasons.qi.allperiod$class = "Quercus"
seasons.allclases.allperiod <- rbind.data.frame(seasons.allclases.allperiod, seasons.P_D.allperiod, seasons.P_ND.allperiod, seasons.qi.allperiod)
plot_day_seasonality_byclass(seasons.allclases.allperiod, "All trees", PLACE, ALL_PERIOD_ST)
if (SAVE) if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-allperiod-seasonalities-alltrees-byclass.png'), width = 15, height = 10)
```

June to July:
```{r}
seasons.allclases.jun_jul <- data.frame()
seasons.P_D.jun_jul$class = "D"
seasons.P_ND.jun_jul$class = "ND"
seasons.qi.jun_jul$class = "Quercus"
seasons.allclases.jun_jul <- rbind.data.frame(seasons.allclases.jun_jul, seasons.P_D.jun_jul, seasons.P_ND.jun_jul, seasons.qi.jun_jul)
plot_day_seasonality_byclass(seasons.allclases.jun_jul, "All trees", PLACE, "June to July 2023")
if (SAVE) {if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-june-july-seasonalities-alltrees-byclass.png'), width = 15, height = 10)}
```
# Climate data

Importing environmental data and keeping the one sensor with all valid data
```{r}
db.env <- read.env.proc(file.path(PATH,ENV_DIR,'proc-env.csv'))
db.env <- db.env[db.env$series == SELECTED_TMS,]
```

Filter env data to period of study
```{r}
db.env <- db.env[which(db.env$ts>=ts_start & db.env$ts<=ts_end),]
summary(db.env)
```

Calculations of useful aggregation data for climate:
  * soil moisture: max VWC, range VWC, mean of VWC (misleading)
  * temp: Mean of daily temp, minimum daily temp and maximum daily temp, range temp and interquartile temp (Q3-Q1). We will use soil temperature as it's used on many other studies and it isn't so much influenced by sun irradiation and night/day differences.

```{r}
clim.daily <- db.env %>% mutate (date = date(ts)) %>% group_by(date) %>% summarise(max.temp = max(surface.temp), min.temp = min(surface.temp), mean.temp = mean(surface.temp), sd.temp = sd(surface.temp), range.vwc = max(vwc) - min(vwc), max.vwc = max(vwc), mean.vwc = mean(vwc), range.temp = max.temp - min.temp, interquartil.temp = as.numeric(quantile(surface.temp, prob=c(.75)) - quantile(surface.temp, prob=c(.25)) ) )
summary(clim.daily)
```

# Correlations: Seasonality ~ climate

Now, let's see if we can see a relation within Seasonality and climate.

## All trees
Let's add up all trees seasonality to a matrix with env data and explore the correlations:
```{r}
matrix <- inner_join(seasonalities.P_D.agg, db.env, by='ts') %>% select(soil.temp, surface.temp, air.temp, vwc, mean) %>% rename(seasonal.value = mean)
matrix
```

```{r}
library("correlation")

results <- correlation(matrix)
results
```

```{r}
library(see)

results %>%
  summary(redundant = TRUE) %>%
  plot()
```
Now using spearman method:
```{r}
results <- correlation(matrix, method = "spearman")

results %>%
  summary(redundant = TRUE) %>%
  plot()
```
Doing cross-correlation:
```{r}
ccf(matrix$seasonal.value, matrix$surface.temp, lag.max = 102)

find_Max_CCF(matrix$seasonal.value, matrix$surface.temp)
```

## By class

Let's test correlations within temp and seasonality for each different class:

### Quercus Ilex
```{r}
matrix <- inner_join(seasonalities.qi.agg, db.env, by='ts') %>% select(soil.temp, surface.temp, air.temp, vwc, mean) %>% rename(seasonal.value = mean)
matrix

results <- correlation(matrix, method = "pearson")

results %>%
  summary(redundant = TRUE) %>%
  plot()
```

Doing cross-correlation:
```{r}
ccf(matrix$seasonal.value, matrix$surface.temp, lag.max = 102)

find_Max_CCF(matrix$seasonal.value, matrix$surface.temp)
```

We found that after 6 instances of 15 mins (1h and half), there is 0.70 inverse correlation within surface temperature and Quercus ilex maximum swelling.

Let's plot this to visually see it for June and July Months (Daily mean seasonality + Daily temperature fluctuations)
```{r}
temp.jun_jul <- db.env %>% filter( (ts >= as.Date("2023-06-01") ) & (ts < as.Date("2023-08-01")) ) %>%  mutate(timeOfDay = as_hms(ts)) %>% summarise(meanTemp = mean(surface.temp), se_temp = (sd(surface.temp) / sqrt(n()) ), .by = timeOfDay)
temp.jun_jul

plot_day_seasonality_and_temp(seasons.qi.jun_jul, temp.jun_jul, "Quercus Ilex", PLACE, "June to July 2023")
# if (SAVE) if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-june-july-seasonalities&Temp-Qi.png'), width = 15, height = 10)
```

### Non-Declining Pines
```{r}
matrix <- inner_join(seasonalities.P_ND.agg, db.env, by='ts') %>% select(soil.temp, surface.temp, air.temp, vwc, mean) %>% rename(seasonal.value = mean)
matrix

results <- correlation(matrix, method = "spearman")

results %>%
  summary(redundant = TRUE) %>%
  plot()
```

Doing cross-correlation:
```{r}
ccf(matrix$seasonal.value, matrix$surface.temp, lag.max = 102)

find_Max_CCF(matrix$seasonal.value, matrix$surface.temp)
```

Let's plot this to visually see it for June and July Months (Daily mean seasonality + Daily temperature fluctuations)
```{r}
plot_day_seasonality_and_temp(seasons.P_ND.jun_jul, temp.jun_jul, "Non-Declining Pines", PLACE, "June to July 2023")
# if (SAVE) if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-june-july-seasonalities&Temp-P_ND.png'), width = 15, height = 10)
```

### Declining Pines
```{r}
matrix <- inner_join(seasonalities.P_D.agg, db.env, by='ts') %>% select(soil.temp, surface.temp, air.temp, vwc, mean) %>% rename(seasonal.value = mean)
matrix

results <- correlation(matrix, method = "spearman")

results %>%
  summary(redundant = TRUE) %>%
  plot()
```

Doing cross-correlation:
```{r}
ccf(matrix$seasonal.value, matrix$surface.temp, lag.max = 102)

find_Max_CCF(matrix$seasonal.value, matrix$surface.temp)
```

Let's plot this to visually see it for June and July Months (Daily mean seasonality + Daily temperature fluctuations)
```{r}
plot_day_seasonality_and_temp(seasons.P_D.jun_jul, temp.jun_jul, "Declining Pines", PLACE, "June to July 2023")
# if (SAVE) ggsave(glue('output/{PLACE}_aggoneday-june-july-seasonalities&Temp-P_D.png'), width = 15, height = 10)
```

# Exploring wavelets

Here we test if there's another correlation with a longer frequency than the daily frequency between temperature and daily seasonality. We could think that, in warm periods (days, weeks, months,...) like the summer season we have a correlation with an specific seasonality.
I would expect to have a year correlation too (every summer we have similar tree), but unfortunately we don't have enough data to test this yet.

Next steps take a long while so,
```{r}
stop('Wavelet analysis below. Too time-consuming and not really relevant.')
```

## For all trees: surface.temp ~ seasonality

```{r}
surface.temp = db.env %>% 
           # select data at hourly intervals by
           # create a new variable named "minutes"
           mutate (minutes = minute (ts)) %>% # Add variable minutes
           # filter "oclock data"
           filter (minutes == "0") %>% # Filter hourly data
           # convert datetime to a numeric variable
           # to make it more intuitive, we will express time as days,
           # and make time 0 = the first instance (corresponds with ts_start)
           mutate (time.days = (as.numeric(ts)- as.numeric(db.env[1,"ts"]))/(60*60*24)) %>%
           # select the variables of interest (time and air.temp)
           select (time.days, surface.temp)
           
any(is.na(surface.temp))
summary(surface.temp)

seasonality = seasonalities.agg %>% 
           # select data at hourly intervals by
           # create a new variable named "minutes"
           mutate (minutes = minute (ts), value = mean) %>% # Add variable minutes
           # filter "oclock data"
           filter (minutes == "0") %>% # Filter hourly data
           # convert datetime to a numeric variable
           # to make it more intuitive, we will express time as days,
           # and make time 0 = the first instance (corresponds with ts_start)
           mutate (time.days = (as.numeric(ts)- as.numeric(seasonalities.agg[1,"ts"]))/(60*60*24)) %>%
           # select the variables of interest (time and air.temp)
           select (time.days, value)
seasonality <- data.frame(seasonality)

any(is.na(seasonality))
summary(seasonality)

# wavelet analysis (this will take a while)
wavelet = wtc(surface.temp, seasonality,
              max.scale = 32,
              #display the progress bar
              quiet = F)

# wavelet plot
par(oma = c(0, 0, 0, 1), mar = c(5, 4, 4, 5) + 0.1)
plot(wavelet,
     #legend colors
     plot.cb = TRUE,
     #lag phase (lines)
     plot.phase = TRUE, 
     ylab = "Frequency (days)",
     xlab = "Time (days)") 
```

## By class: surface.temp ~ seasonality

For Quercus:
```{r}
seasonality = seasonalities.qi.agg %>% 
           # select data at hourly intervals by
           # create a new variable named "minutes"
           mutate (minutes = minute (ts), value = mean) %>% # Add variable minutes
           # filter "oclock data"
           filter (minutes == "0") %>% # Filter hourly data
           # convert datetime to a numeric variable
           # to make it more intuitive, we will express time as days,
           # and make time 0 = the first instance (corresponds with ts_start)
           mutate (time.days = (as.numeric(ts)- as.numeric(seasonalities.agg[1,"ts"]))/(60*60*24)) %>%
           # select the variables of interest (time and air.temp)
           select (time.days, value)
seasonality <- data.frame(seasonality)

any(is.na(seasonality))
summary(seasonality)

# wavelet analysis (this will take a while)
wavelet = wtc(surface.temp, seasonality,
              max.scale = 32,
              #display the progress bar
              quiet = F)

# wavelet plot
par(oma = c(0, 0, 0, 1), mar = c(5, 4, 4, 5) + 0.1)
plot(wavelet,
     #legend colors
     plot.cb = TRUE,
     #lag phase (lines)
     plot.phase = TRUE, 
     ylab = "Frequency (days)",
     xlab = "Time (days)") 
```

For Non-Declining pines:
```{r}
seasonality = seasonalities.P_ND.agg %>% 
           # select data at hourly intervals by
           # create a new variable named "minutes"
           mutate (minutes = minute (ts), value = mean) %>% # Add variable minutes
           # filter "oclock data"
           filter (minutes == "0") %>% # Filter hourly data
           # convert datetime to a numeric variable
           # to make it more intuitive, we will express time as days,
           # and make time 0 = the first instance (corresponds with ts_start)
           mutate (time.days = (as.numeric(ts)- as.numeric(seasonalities.agg[1,"ts"]))/(60*60*24)) %>%
           # select the variables of interest (time and air.temp)
           select (time.days, value)
seasonality <- data.frame(seasonality)

any(is.na(seasonality))
summary(seasonality)

# wavelet analysis (this will take a while)
wavelet = wtc(surface.temp, seasonality,
              max.scale = 32,
              #display the progress bar
              quiet = F)

# wavelet plot
par(oma = c(0, 0, 0, 1), mar = c(5, 4, 4, 5) + 0.1)
plot(wavelet,
     #legend colors
     plot.cb = TRUE,
     #lag phase (lines)
     plot.phase = TRUE, 
     ylab = "Frequency (days)",
     xlab = "Time (days)") 
```

For Declining pines:
```{r}
seasonality = seasonalities.P_D.agg %>% 
           # select data at hourly intervals by
           # create a new variable named "minutes"
           mutate (minutes = minute (ts), value = mean) %>% # Add variable minutes
           # filter "oclock data"
           filter (minutes == "0") %>% # Filter hourly data
           # convert datetime to a numeric variable
           # to make it more intuitive, we will express time as days,
           # and make time 0 = the first instance (corresponds with ts_start)
           mutate (time.days = (as.numeric(ts)- as.numeric(seasonalities.agg[1,"ts"]))/(60*60*24)) %>%
           # select the variables of interest (time and air.temp)
           select (time.days, value)
seasonality <- data.frame(seasonality)

any(is.na(seasonality))
summary(seasonality)

# wavelet analysis (this will take a while)
wavelet = wtc(surface.temp, seasonality,
              max.scale = 32,
              #display the progress bar
              quiet = F)

# wavelet plot
par(oma = c(0, 0, 0, 1), mar = c(5, 4, 4, 5) + 0.1)
plot(wavelet,
     #legend colors
     plot.cb = TRUE,
     #lag phase (lines)
     plot.phase = TRUE, 
     ylab = "Frequency (days)",
     xlab = "Time (days)") 
```