library(ggplot2)
library(dplyr)
library(tidyverse)
library(glme)
library(ggpmisc)
#### script that loads all abiotic data, the flux data and scales the data, and also performs different approaches with LMEs
# load following data:
# tachymeter data merged with ch4
setwd("C:/Users/Carolin/Dropbox/Diss/MS/alder_short_communication")
tachy.ch4 <- read.csv("tachy_ch4.csv", header = T, sep = ",", dec = ".")

# load TWI
twi <- read.csv("add_trees_TWI.csv", header =T, sep =";", dec =".")
twi$collar_id.TWI <- names("collar_id")

ch4 <- merge(x = tachy.ch4, y = twi, by = "collar_id", all.x = TRUE)

ch4$CH4.flux <- ch4$CH4.flux*1000

# load CO2
setwd("C:/Users/Carolin/Dropbox/Daten/Fluxes/Treefluxes/Reco")
co2 <- msmts_CO2 <- read.table("tree_data_CO2.txt", header = T, sep ="\t", dec =".")

ch4.2 <- ch4 %>% left_join(co2, by=c("collar_id","meas_date"))

# all numerical data is included
# scale all explanatory variables
ch4.2$hasb2 <- scale(ch4.2$hasb.x,center = T, scale = T)
ch4.2$dbh2 <- scale(ch4.2$dbh.x,center = T, scale = T)
ch4.2$CO2.flux2 <- scale(ch4.2$CO2.flux,center = T, scale = T)
ch4.2$TWI2 <- scale(ch4.2$TWI,center = T, scale = T)
ch4.2$ortho_height2 <- scale(ch4.2$ortho_height,center = T, scale = T)
ch4.2$Rechtswert2 <- scale(ch4.2$Rechtswert,center = T, scale = T)
ch4.2$Hochwert2 <- scale(ch4.2$Hochwert,center = T, scale = T)
ch4.2$CH4.flux2 <- scale(ch4.2$CH4.flux, center = T, scale = T)


#### add other fixed effects ####
ch4.2$leaf.out <- with(ch4.2, ifelse(meas_date =="2019-04-24", "pre leafout",
                                             "post leafout"))

###################################################################

# load environmental data
#### water####
setwd("C:/Users/Carolin/Dropbox/Daten/abiotic_data/waterlevel")
#setwd("C:/Users/koehn/Dropbox/Daten/abiotic_data/waterlevel")

wl <- read.table("1909_groundwater_15_min.txt", header = T, dec = ".", sep = ",",
                 strip.white = T)
wl_AW <- subset(wl, siteID == "AW")
wl_AW$date <- substr(wl_AW$datetime,0,10)
wl_AW$date <- as.Date(wl_AW$date)
# choose the according period, two weeks prior measurement campaign
wl_AW_rel <- subset(wl_AW, date >= "2019-04-10" & date <= "2019-05-28")
wl_AW_rel$datetime <- as.POSIXct(wl_AW_rel$datetime, format = "%Y-%m-%d %H:%M:%S")
#rename
names(wl_AW_rel)[2] <- "datetime_utc"


#### soiltemp and maybe melt ####
setwd("C:/Users/Carolin/Dropbox/Daten/abiotic_data/soiltemp")
#setwd("C:/Users/koehn/Dropbox/Daten/abiotic_data/soiltemp")
st <- read.table("soiltemp_complete.txt", header = T, dec =".", sep=",", strip.white = T)


st$datetime_utc <- as.POSIXct(st$datetime_utc, format = "%Y-%m-%d %H:%M:%S")
st$date <- as.Date(substr(st$datetime_utc, 0,10))

#AW
st_AW <- subset(st, site == "AW")
#choose the correct part of the data
st_AW <- subset(st_AW, datetime_utc >= "2019-04-10 00:00:00" & datetime_utc <= "2019-05-28 23:59:59")
# choose the according period, two weeks prior measurement campaign



st_15 <- st_AW[,c(2,3,5)]
st_15$depth <- 15
st_5 <- st_AW[,c(2,3,4)]
st_5$depth <-5

# aggregate to only have an average soiltemp at 15 cm # 
st_15_agg <- aggregate(soil_t_15 ~ datetime_utc, data = st_15, FUN= "mean" )

# same for st_5 #
st_5_agg <- aggregate(soil_t_5 ~ datetime_utc, data = st_5, FUN= "mean" )

# join soil temp depths again #
st_final <- merge(st_5_agg, st_15_agg, by = "datetime_utc")

### combine soiltemp and wl ###

st_wl <- merge(st_final, wl_AW_rel, by =  "datetime_utc")

#### sap data ####
#### read in sap flow data ####
setwd("C:/Users/Carolin/Dropbox/Daten/abiotic_data/tree")
sap <- read.table("final_sap.csv", header = T, sep = ",", dec = ".")
# get the time stamp right
#sap$date <- substr(sap$timestamp,1,10)
#sap$date <- as.Date(sap$date, "%Y-%m-%d")
#sap$time <- substr(sap$timestamp,11,18)
# get datetime together for export
#sap$datetime <- paste(sap$date,sap$time)
sap$datetime <- as.POSIXct(sap$datetime, format = "%Y-%m-%d %H:%M:%S")
hist(sap$Umean)


#### check some overall models ####
prelim_plot <- ggplot(ch4.2, aes(x = ortho_height2, y = CH4.flux2)) +
   geom_point() +
   geom_smooth(method = "lm")
prelim_plot
hist(ch4.2$CH4.flux2)
pp <- ggplotly(prelim_plot)

# closer look at lms
basic.lm <- lm(CH4.flux ~ ortho_height2, data = ch4.2)
summary(basic.lm)

plot(basic.lm, which = 2)


# look at the independence of all samples
boxplot(CH4.flux ~ meas_date, data = ch4.2)  # certainly looks like something is going on here

# further watching
colour_plot <- ggplot(ch4.2, aes(x = ortho_height2, y = CH4.flux, colour = hummock.x, 
                                 shape = meas_date)) +
    geom_point(size = 3) +
    theme(legend.position = "none")

colour_plot + theme(legend.position = "top")  

#### multiple analyses based on meas_date ####
split_plot <- ggplot(aes(CO2.flux, CH4.flux), data = ch4.2) + 
   geom_point() + 
   facet_wrap(~ tree_kind.x) + # create a facet for each mountain range
   xlab("CO2 flux") + 
   ylab("CH4 flux")

split_plot

#### fixed effect ####
meas_date_lm <- lm(log(CH4.flux) ~ ortho_height2, data = ch4.2)
summary(meas_date_lm)

# create subsets to do some simple lms
ch4.514 <- subset(ch4, meas_date == "2019-05-14")
ch4.424 <- subset(ch4, meas_date == "2019-04-24")
ch4.528 <- subset(ch4, meas_date == "2019-05-28")


# lms

lm <- lm(CH4.flux ~ hasb, data = ch4.528)
summary(lm)


#### multiple lm plots ####
#### orthometric height ####
ch4_ortho <- ggplot(ch4.2, aes(x=CH4.flux, y = ortho_height, colour = as.factor(meas_date),shape = as.factor(meas_date)))+
  geom_point(size = 4)+
  geom_smooth(method='lm', formula= y~x, se = F)+
  theme_classic()+
  theme(axis.text.x=element_text(angle=0, size=20, colour = "black"))+
  theme(axis.text.y=element_text(angle=0, size=20, colour = "black"))+
  theme(axis.title.y=element_text(angle=90, size=24, colour = "black"))+
  theme(axis.title.x=element_text(angle=0, size=24, colour = "black"))+
  scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+
  #scale_y_continuous(breaks = seq(0,500,100), limits = c(-20,500))+
  #labs(x=expression(paste(CO[2]~'flux [',mg,~m^-2,~h^-1,']')),y=expression(paste(CH[4]~'flux [',mu~g,~m^-2,~h^-1,']')))+
  #labs(x="Height above stem base",y=expression(paste(CH[4]~'flux [',mg,~m^-2,~h^-1,']')))+
  labs(x=expression(paste(CH[4]~'flux [',mu~g,~m^-2,~h^-1,']')),y="Orthometric height [m.a.s.l]")+
  #scale_x_continuous(breaks = seq(0,200,50), limits = c(0,220))+
  #labs(x="TWI",y=expression(paste(CH[4]~'flux [',mu~g,~m^-2,~h^-1,']')))+
  #geom_hline(yintercept=0, colour="grey")+
  theme(legend.position = c(0.8, 0.8), legend.title = element_blank(),
        legend.direction = "vertical")+
  theme(legend.text = element_text(size = 18))+
  theme(panel.border= element_rect(colour = "black", fill=NA, size=1))+
  guides(color = guide_legend(override.aes = list(size=7)))+
  stat_poly_eq(formula = y~x, 
               aes(label = paste(..rr.label.., sep = "~~~")), 
               parse = TRUE, label.x = 0.3, size = 8)
print(ch4_ortho)

#### co2 and ch4 scatter plot ####
co2_ch4 <- ggplot(ch4.2, aes(x=CH4.flux, y = CO2.flux, colour = as.factor(meas_date), shape = as.factor(meas_date)))+
  geom_point(size = 4)+
  geom_smooth(method='lm', formula= y~x, se = F)+
  theme_classic()+
  theme(axis.text.x=element_text(angle=0, size=20, colour = "black"))+
  theme(axis.text.y=element_text(angle=0, size=20, colour = "black"))+
  theme(axis.title.y=element_text(angle=90, size=24, colour = "black"))+
  theme(axis.title.x=element_text(angle=0, size=24, colour = "black"))+
  scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+
  #scale_y_continuous(breaks = seq(0,500,100), limits = c(-20,500))+
  labs(x=expression(paste(CH[4]~'flux [',mu~g,~m^-2,~h^-1,']')),
       y=expression(paste(CO[2]~'flux [',mg,~m^-2,~h^-1,']')))+
  
  #labs(x="Height above stem base",y=expression(paste(CH[4]~'flux [',mg,~m^-2,~h^-1,']')))+
  #labs(x="Height above stem base [m]",y=expression(paste(CH[4]~'flux [',mu~g,~m^-2,~h^-1,']')))+
  #scale_x_continuous(breaks = seq(0,200,50), limits = c(0,220))+
  #labs(x="TWI",y=expression(paste(CH[4]~'flux [',mu~g,~m^-2,~h^-1,']')))+
  geom_hline(yintercept=0, colour="grey")+
  theme(legend.position = c(0.8, 0.4), legend.title = element_blank(),
        legend.direction = "vertical")+
  theme(legend.text = element_text(size = 18))+
  theme(panel.border= element_rect(colour = "black", fill=NA, size=1))+
  guides(color = guide_legend(override.aes = list(size=7)))+
  stat_poly_eq(formula = y~x, 
               aes(label = paste(..rr.label.., sep = "~~~")), 
               parse = TRUE, label.x = 0.3, size = 8)
print(co2_ch4)


#### scatter hasb ####
ch4_hasb <- ggplot(ch4.2, aes(x=CH4.flux, y = hasb.x, colour = as.factor(meas_date), shape = as.factor(meas_date)))+
  geom_point(size = 4)+
  geom_smooth(method='lm', formula= y~x, se = F)+
  theme_classic()+
  theme(axis.text.x=element_text(angle=0, size=20, colour = "black"))+
  theme(axis.text.y=element_text(angle=0, size=20, colour = "black"))+
  theme(axis.title.y=element_text(angle=90, size=24, colour = "black"))+
  theme(axis.title.x=element_text(angle=0, size=24, colour = "black"))+
  scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+
  #scale_y_continuous(breaks = seq(0,500,100), limits = c(-20,500))+
  #labs(x=expression(paste(CO[2]~'flux [',mg,~m^-2,~h^-1,']')),y=expression(paste(CH[4]~'flux [',mu~g,~m^-2,~h^-1,']')))+
  #labs(x="Height above stem base",y=expression(paste(CH[4]~'flux [',mg,~m^-2,~h^-1,']')))+
  labs(x=expression(paste(CH[4]~'flux [',mu~g,~m^-2,~h^-1,']')),y="Height above stem base [cm]")+
  #scale_x_continuous(breaks = seq(0,200,50), limits = c(0,220))+
  #labs(x="TWI",y=expression(paste(CH[4]~'flux [',mu~g,~m^-2,~h^-1,']')))+
  #geom_hline(yintercept=0, colour="grey")+
  theme(legend.position = c(0.8, 0.8), legend.title = element_blank(),
        legend.direction = "vertical")+
  theme(legend.text = element_text(size = 18))+
  theme(panel.border= element_rect(colour = "black", fill=NA, size=1))+
  guides(color = guide_legend(override.aes = list(size=7)))+
  stat_poly_eq(formula = y~x, 
               aes(label = paste(..rr.label.., sep = "~~~")), 
               parse = TRUE, label.x = 0.3, size = 8)
print(ch4_hasb)


#### bind the two scattrplots together ####
plot_grid <- plot_grid(ch4_ortho, co2_ch4,ch4_hasb,
                       align = "hv", axis = "b", nrow = 2,
                       labels = c("(a)","(b)","(c)"))


setwd("C:/Users/Carolin/Dropbox/Diss/MS/alder_short_communication/figures")
#setwd("C:/Users/koehn/Dropbox/Diss/MS/alder_short_communication/figures")

ggsave("20250115_scatter.png", plot = plot_grid, device = "png",
       scale = 1, width = 360, height = 300, units = "mm",
       dpi = 500)



# load the interpolated water data from "models with interpolated waterlevels.R"
# and merge by collar id and meas_date
ch4.exp <- ch4.2 %>%
  select(collar_id,meas_date,CH4.flux,hasb2,dbh2,hummock.x,ortho_height2,TWI2,CO2.flux2, leaf.out) %>%
  inner_join(dat_4242 %>% select(collar_id,meas_date,interpolated_waterlevel_dist4_424), by = c("collar_id","meas_date"))
# adding an ID
x <- 1
ch4.exp$counter <-  x + (1:nrow(ch4.exp)) + 1

ch4.exp2 <- ch4.2 %>%
  select(collar_id,meas_date,CH4.flux,hasb2,dbh2,hummock.x,ortho_height2,TWI2,CO2.flux2, leaf.out) %>%
  inner_join(dat_5142 %>% select(collar_id,meas_date,interpolated_waterlevel_dist4_514), by = c("collar_id","meas_date"))
x <- 30
ch4.exp2$counter <-  x + (1:nrow(ch4.exp2)) + 1

ch4.exp3 <- ch4.2 %>%
  select(collar_id,meas_date,CH4.flux,hasb2,dbh2,hummock.x,ortho_height2,TWI2,CO2.flux2, leaf.out) %>%
  inner_join(dat_5282 %>% select(collar_id,meas_date,interpolated_waterlevel_dist4_528), by = c("collar_id","meas_date"))
x <- 60
ch4.exp3$counter <-  x + (1:nrow(ch4.exp3)) + 1

# rbind these merged dataframes

ch4.exp_fin <- bind_rows(ch4.exp, ch4.exp2, ch4.exp3)

# add a variable if tree was most likely inundated
ch4.exp_fin$inundated[ch4.exp_fin$meas_date == "2019-05-14"] <- "yes"

#### final step get the soil temperature data together with the flux data ####
# loaded the soil temperature data from the 140 script
# st_final
# clculate daily average value from both depths from st_final
st_final$date <- as.Date(substr(st_final$datetime_utc,0,10))
st_final$date <- as.factor(st_final$date)
st_final$date <- as.Date(st_final$date)

setDT(st_final)

# Calculate daily averages across multiple columns
df_avg <- aggregate(cbind(soil_t_5, soil_t_15) ~ date, data = st_final, FUN = mean, na.rm = TRUE)

# add the values via date in ch4.exp_fin
ch4.exp_fin$st_5[ch4.exp_fin$meas_date == "2019-04-24"] <- 11
ch4.exp_fin$st_15[ch4.exp_fin$meas_date == "2019-04-24"] <- 9.9


ch4.exp_fin$st_5[ch4.exp_fin$meas_date == "2019-05-14"] <- 9.9
ch4.exp_fin$st_15[ch4.exp_fin$meas_date == "2019-05-14"] <- 9.6

ch4.exp_fin$st_5[ch4.exp_fin$meas_date == "2019-05-28"] <- 12.8
ch4.exp_fin$st_15[ch4.exp_fin$meas_date == "2019-05-28"] <- 12.4


#### mixed effect models #### -> random effect
library(lme4)


# should the variables be random or fixed effects?!
# -> rather fixed effects: firstly because we believe they can have an impact
# and secondly because it doesnt make sense to use to factoral variables as random 
# effects
# add categorical column emitter
ch4.exp_fin$emitter <- with(ch4.exp_fin, ifelse(collar_id =="AWT-001808", "high emitter",
                           ifelse(collar_id =="AWT-001965", "high emitter",
                                  ifelse(collar_id =="AWT-001724", "high emitter", 
                                         ifelse(collar_id =="AWT-001738", "high emitter",
                                                ifelse(collar_id =="AWT-001968", "high emitter",
                                                       ifelse(collar_id =="AWT-001810", "high emitter", "low emitter")))))))


mod<-lmer(CH4.flux ~ (1|collar_id)+ortho_height2+CO2.flux2+TWI2+
            leaf.out+hasb2+inundated,data=ch4.exp_fin)
summary(mod)$coefficients
summary(mod)

anova(mod)

predicted <- predict(mod)
ch4.2 <- ch4.2[-3,]
df <- data.frame(Observed = ch4.exp_fin$CH4.flux,
                 Predicted = predicted)
ggplot(df, aes(x = Observed, y = Predicted)) +
  geom_point() +
  geom_abline(intercept = 3, slope = 1,
              linetype = "dashed") +
  xlab("Observed CH4 flux") +
  ylab("Predicted CH4 flux")
# get coefficients
coeffs <- coef(summary(mod))
p <- pnorm(abs(coeffs[, "t value"]), lower.tail = FALSE) * 2
cbind(coeffs, "p value" = round(p,5))

qqnorm(resid(mod))
qqline(resid(mod))  # points fall nicely onto the line - good!

# using the lmerTest to do formal testing on the lme
library(lmerTest)
library(ggeffects)
library(performance)

summary(mod)

r2(mod)

##################################################################
# delete some rows due to missing TWI2 values
ch4.3 <- ch4.2[-(13:15),]

#### plot fitted vs real
plot(predict(mod),ch4.3$CH4.flux2,
     xlab="predicted",ylab="actual")
abline(a=0,b=1)

#### create an ordered plot ####
ch4$collar_id <- as.factor(ch4$collar_id)
ch4.avg <- aggregate(CH4.flux ~ collar_id, data = ch4, FUN = mean)
ch4.avg$CH4.flux <- ch4.avg$CH4.flux*1000

ch4.ordered <- ch4.avg[order(ch4.avg$CH4.flux), ]
ch4.ordered$names <- ch4.ordered$collar_id

# now plot it with bar plot
library(ggrepel)
ch4.ordered$label <- ifelse(ch4.ordered$CH4.flux > 40, as.character(ch4.ordered$names), NA)

ordered <- ggplot(ch4.ordered, aes(x = reorder(collar_id, -CH4.flux), y = CH4.flux)) +
  geom_hline(yintercept = 0, col = "red")+
  geom_bar(stat = "identity", fill = "slategrey", colour = "black") +
  theme_classic(base_size = 15)+
  theme(
    axis.text.x = element_blank(),  # Suppress x-axis text
    axis.ticks.x = element_blank(),  # Suppress x-axis title
  ) +
  labs(x = "Collar ID",y=expression(paste(CH[4]~'flux [',mu~g,~m^-2,~h^-1,']')))+
  geom_text(aes(label = label, vjust = -.25, hjust = -.1), size = 4)
print(ordered)
#### save plot in figures ####
setwd("C:/Users/Carolin/Dropbox/Diss/MS/alder_short_communication/figures")
#setwd("C:/Users/koehn/Dropbox/Diss/MS/alder_short_communication/figures")

ggsave("20241217_ordered_plot.png", plot = ordered, device = "png",
       scale = 1, width = 180, height = 120, units = "mm",
       dpi = 500)