# 1. Preparing data
# cleaning the environment, loading packages
rm(list=ls(all=TRUE))

par.sav <- par(no.readonly = T)

library(readxl)    # reading excel files
library(chron)     # for handling dates
library(lubridate) # for calculating day of the year
library(lomb)      # for fitting a Lomb-Scargle periodogram
library(zoo)       # for calculating rolling mean

# preparing tick information for the plots
begin.years <- dates(c("01/01/2010", "01/01/2011",
                       "01/01/2012", "01/01/2013",
                       "01/01/2014", "01/01/2015",
                       "01/01/2016", "01/01/2017",
                       "01/01/2018", "01/01/2019",
                       "01/01/2020", "01/01/2021",
                       "01/01/2022", "01/01/2023",
                       "01/01/2024", "01/01/2025",
                       "01/01/2026"))  

begin.years.long <- dates(c("01/01/10",
                            "01/01/14",
                            "01/01/18",
                            "01/01/22",
                            "01/01/26"))

begin.years.short <- setdiff(begin.years, begin.years.long)
## short and long refer to the tick length

#######################################################
# 2 Reading daily temperature data June 01, 2009 to 
#   December 31, 2025

datT <- read_xlsx("temperature.xlsx")
datT$date <- as.numeric(chron(as.character(datT$Zeitstempel),
                          format=c(dates="d.m.y"))
                    - dates("12/31/2009"))
  ## changing date to days since 31.12.2009

resTrend <- lm (Wert ~ date, data = datT)
summary(resTrend)


#####################################################
# 3a plotting temperature data  
#    using all days
par(mar = c(5, 6, 1, 1), mai = c(1, 1, 0.2, 0.2),
    lwd = 2, las = 1,cex = 0.9, 
    cex.axis = 1, cex.lab = 1.5)

plot(datT$date, datT$Wert, axes = F,
     xlim = c(-300, 5850),
     ylab = "Temperature [°C]",
     xlab = " ",
     pch = 21, col = "darkgrey", bg = "darkgrey")
yticks <- c(-10, 0, 10, 20)
axis(1, at = (begin.years.long - dates("12/31/2009")), 
     labels = begin.years.long, lwd = 2, tck = -0.07)
axis(1, at = (begin.years.short - dates("12/31/2009")), 
     labels = rep("", length(begin.years.short)), lwd = 2,
     tck = -0.04)
axis(2, at = yticks, labels = yticks, lwd =2)
box(lwd =2)
abline(v=(begin.years - dates("12/31/2009")), col = "grey")

rollLine <- rollmean(cbind(datT$date, datT$Wert), k= 3)
lines(rollLine[,1], rollLine[,2], col = "red", lwd = 2)
par(par.sav)

#####################################################
# 3b plotting temperature data  
#    using only days with counts

par(mar = c(5, 6, 1, 1), mai = c(1, 1, 0.2, 0.2),
    lwd = 2, las = 1,cex = 0.9, 
    cex.axis = 1, cex.lab = 1.5)

# finding temperature data of days with counts
dataDate <- names(read_xlsx("transectCountsDeps2010_2025.xlsx"))
dataDate <- dataDate[-c(1 : 3)]
dataDate <- as.numeric(chron(as.character(dataDate),
                              format=c(dates="d.m.y"))
                        - dates("12/31/2009"))
scr <- match(dataDate, datT$date)
datTCounts <- datT[scr, ]

plot(datTCounts$date, datTCounts$Wert, axes = F,
     xlim = c(-0, 5850),
     ylab = "Temperature [°C]",
     xlab = " ",
     pch = 21, col = "darkgrey", bg = "darkgrey")
yticks <- c(-10, 0, 10, 20)
axis(1, at = (begin.years.long - dates("12/31/2009")), 
     labels = begin.years.long, lwd = 2, tck = -0.07)
axis(1, at = (begin.years.short - dates("12/31/2009")), 
     labels = rep("", length(begin.years.short)), lwd = 2,
     tck = -0.04)
axis(2, at = yticks, labels = yticks, lwd =2)
box(lwd =2)
abline(v=(begin.years - dates("12/31/2009")), col = "grey")

rollLine <- rollmean(cbind(datTCounts$date, datTCounts$Wert), k= 3)
lines(rollLine[,1], rollLine[,2], col = "red", lwd = 2)
par(par.sav)

#################################################
# 3 Calculating Lomb-Scargle-Periodogram with all 
#   temperature data

data.lomb <- data.frame(date = datT$date,
                        temp = datT$Wert)
res.lomb <-lsp(data.lomb, type = "period",
               from = 20, to = 1000, ofac = 10, plot = F)

par(mar = c(5, 6, 1, 1), mai = c(1, 1, 0.2, 0.2),
    lwd = 2, las = 1, cex = 0.9, 
    cex.axis = 1, cex.lab = 1.5)

plot(res.lomb$scanned, 
     res.lomb$power,
     type = "l",
     xlab="Period [days]", ylab="Normalized power - temperature",
     xlim = c(0, 1000), ylim = c(0, 0.8),
     axes = F)
ticksx<-c(0, 364, 500, 1000)
ticksy<-c(0, 0.2, 0.4, 0.6, 0.8)
axis(1, at=ticksx, labels=ticksx, lwd.ticks=2)
axis(2, at=ticksy, labels=ticksy, lwd.ticks=2)
abline(h = res.lomb$sig.level, col = "red")
text(150,res.lomb$sig.level + 0.08,
     labels = " p < 0.01",
     col = "red", cex = 1)
box()

# plotting periodogram for days with counts 
# into the periodogram with all temperature data 
data.lombR <- data.frame(date = datTCounts$date,
                         temp = datTCounts$Wert)
res.lombR <-lsp(data.lombR, type = "period", 
                from = 20, to = 1000, ofac = 5, plot = F)
par(fig = c(0.6, 1, 0.4, 0.8), mar=c(2, 5, 2, 2), 
    cex.lab = 1, new=TRUE)
plot(res.lombR$scanned, 
     res.lombR$power,
     type = "l",
     xlab="", ylab="Dates with counts",
     xlim = c(0, 1000), ylim = c(0, 0.8),
     col = "blue",
     axes = F)
axis(1, at=364)
axis(2, at = c(0, 0.4, 0.8))
par(par.sav)

############################################
# 4 Reading rainfall data

  ## processing identical to temperature data

datP <- read_xlsx("rainfall.xlsx")

datP$date <- as.numeric(chron(as.character(datP$Zeitstempel),
                             format=c(dates="d.m.y"))
                       - dates("12/31/2009"))

#########################################################
# 5a  Plotting precipitation data using all days

par(mar = c(5, 6, 1, 1), mai = c(1, 1, 0.2, 0.2),
    lwd = 2, las = 1,cex = 0.9, 
    cex.axis = 1, cex.lab = 1.5)

plot(datP$date, datP$Wert, axes = F,
     xlim = c(-300, 5850),
     ylab = "Precipitation [mm]",
     xlab = " Date - starting January 01, 2010",
     pch = 21, col = "darkgrey", bg = "darkgrey")
yticks <- c (0, 10, 20, 30, 40, 50)
axis(1, at = (begin.years.long - dates("12/31/2009")), 
     labels = begin.years.long, lwd = 2, tck = -0.07)
axis(1, at = (begin.years.short - dates("12/31/2009")), 
     labels = rep("", length(begin.years.short)), lwd = 2,
     tck = -0.04)
axis(2, at = yticks, labels = yticks, lwd =2)
box(lwd =2)
abline(v=(begin.years - dates("12/31/2009")), col = "grey")
rollLine <- rollmean(cbind(datP$date, datP$Wert), k= 3)
lines(rollLine[,1], rollLine[,2], col = "red", lwd = 2)
par(par.sav)

############################################
# 5b Periodogram for rainfall data

data.lomb <- data.frame(date = datP$date,
                        precip = datP$Wert)
res.lomb <-lsp(data.lomb, type = "period",
               plot = F)

par(mar = c(5, 6, 1, 1), mai = c(1, 1, 0.2, 0.2),
    lwd = 2, las = 1,cex = 0.9, 
    cex.axis = 1, cex.lab = 1.5)

plot(res.lomb$scanned, 
     res.lomb$power,
     type = "l",
     xlab="Period [days]", 
     ylab = expression(paste("Normalized power x ", 10^-3,sep="")),
     xlim = c(0, 1000), ylim = c(0, 0.006),
     axes = F)
ticksx<-c(0, 182, 364, 500, 1000)
ticksy<-c(0, 0.002, 0.004, 0.004, 0.006)
axis(1, at=ticksx, labels=ticksx, lwd.ticks=2)
axis(2, at=ticksy, labels=ticksy*1000, lwd.ticks=2)
abline(h = res.lomb$sig.level, col = "red")
text(600,res.lomb$sig.level + 0.0005,
     labels = " p < 0.01",
     col = "red", cex = 1)
box()
par(par.sav)