# Seasonal Waddensea Data
# R version 4.0.3 (2020-10-10)
# 04-03-2020
######################################################################################################
# Environmental analysis of Wadden Sea 2019-2020 datasets
######################################################################################################

# Install and load required packages (if not already installed)
required_packages <- c("reshape2", "ggplot2", "ade4", "tidyverse", "ggrepel", "ggdendro",
                       "dendextend", "harrietr", "Hmisc", "colorspace", "oce", "ocedata", "metR",
                       "ggpubr", "rlang")
install_missing_packages <- required_packages[!(required_packages %in% installed.packages()[,"Package"])]
if (length(install_missing_packages) > 0) {
  install.packages(install_missing_packages, repos = "https://cran.rstudio.com/")
}
lapply(required_packages, require, character.only = TRUE)

########### Water column properties ###########

# Data import water column properties
wzdata <- read.csv("/export/lv6/user/tdegroot/Waddensea/wzmetadata.csv", header = TRUE, sep = ",")
alldata <- read.csv("/export/lv6/user/tdegroot/Waddensea/metadata_ws_ns.csv", header = TRUE, sep = ",")
typeof(wzdata)
str(wzdata)
head(wzdata)
summary(wzdata$salinity)
colnames(wzdata)

# Divide data into seasons -- SUBDATA
seasons <- c("winter", "spring", "summer", "autumn")
subdata <- lapply(seasons, function(season) wzdata[wzdata$season == season, ])

# Normalize depth for min value to be 0m or 0.5m
normalize_depth <- function(data) {
  min_values <- aggregate(depth ~ time, data = data, FUN = min)
  data$depth_norm0 <- NA
  for (i in 1:nrow(min_values)) {
    data[data$time == min_values$time[i], "depth_norm0"] <- data[data$time == min_values$time[i], "depth"] - min_values$depth[i]
  }
  return(data)
}

subdata <- lapply(subdata, normalize_depth)

# Function to create profiles
plot_profiles <- function(data, variable, title) {
  ggplot(data = data[!is.na(data$value) & data$variable == variable, ],
         aes(x = time, y = depth, z = value, fill = value)) +
    geom_tile() +
    scale_fill_distiller(name = title, palette = "RdYlBu", direction = -1, na.value = "#00000000", limits = c(14, 24), breaks = seq(14, 24, 3)) +
    scale_y_reverse(breaks = c(0, 2, 4, 6), limits = c(6.5, 0, 5), expand = c(0, 0), sec.axis = dup_axis(name = NULL, labels = NULL)) +
    scale_x_continuous(breaks = c(0, 10, 20, 30, 40, 50), limits = c(-1, 51), expand = c(0, 0), sec.axis = dup_axis(name = NULL, labels = NULL)) +
    labs(title = "", x = "Time (h)", y = "Depth (m)") +
    theme(
      plot.background = element_rect(fill = "transparent", color = "transparent"),
      legend.background = element_rect(fill = "transparent"),
      panel.border = element_rect(fill = "transparent", size = 2),
      panel.background = element_rect(fill = "transparent"),
      panel.grid = element_blank(),
      axis.ticks = element_line(size = 1, color = "grey1"),
      axis.ticks.length = unit(-0.25, "cm"),
      axis.text.x = element_text(margin = unit(c(0.3, 0.3, 0.3, 0.3), "cm")),
      axis.text.y = element_text(margin = unit(c(0.3, 0.3, 0.3, 0.3), "cm")),
      axis.text = element_text(size = 16),
      legend.title = element_text(size = 18, face = 'bold'),
      legend.text = element_text(size = 16),
      axis.title = element_text(size = 18),
      plot.title = element_text(size = 18, face = 'bold'),
      strip.text = element_text(size = 18)
    ) +
    coord_fixed(ratio = 3)
}

temp_profiles <- lapply(subdata, plot_profiles, variable = "temperature", title = "")
dens_profiles <- lapply(subdata, plot_profiles, variable = "density", title = "Density (PSU)")
sal_profiles <- lapply(subdata, plot_profiles, variable = "salinity", title = "Salinity")

temp_all <- ggpubr::ggarrange(plotlist = temp_profiles, ncol = 4, align = "hv", common.legend = TRUE)
dens_all <- ggpubr::ggarrange(plotlist = dens_profiles, ncol = 1, align = "hv", common.legend = FALSE)
sal_all <- ggpubr::ggarrange(plotlist = sal_profiles, ncol = 1, align = "hv", common.legend = FALSE)

temp_all
dens_all
sal_all

ggpubr::ggarrange(temp_all, sal_all, dens_all, ncol = 3)
ggpubr::ggarrange(temp_all, ncol = 1)

########## Analysis CH4 related parameters water column and sediments #######################################

# Data import and format
wzdatach4 <- read.csv("/export/lv6/user/tdegroot/Waddensea/datach4wz.csv", header = TRUE, sep = ",")
wzdatach4$ch4 <- as.numeric(wzdatach4$ch4)
wzdatach4$depth <- as.numeric(wzdatach4$depth)
wzdatach4$depth <- factor(wzdatach4$depth, levels = c(-1, -3), labels = c(-1, -3), ordered = TRUE)

# Function to create CH4 plots
create_ch4_plot <- function(data, title_suffix) {
  ggplot(data = data, aes(x = time, y = ch4, color = depth)) +
    facet_grid(season ~ variable, scales = "free_y") +
    scale_color_manual(values = c("coral", "seagreen")) +
    geom_point(size = 2) +
    geom_path(size = 1) +
    scale_x_continuous(breaks = seq(0, 50, 10), limits = c(0, 50), expand = c(0, 0)) +
    scale_y_continuous(breaks = seq(0, 150, 30), limits = c(0, 160), expand = c(0, 0)) +
    xlab("Time (h)") + ylab("nM/day") +
    ggtitle(paste("CH4 Conc. -", title_suffix)) +
    common_theme()
}

# Function to create mox and k plots
create_mox_k_plot <- function(data, parameter) {
  ggplot(data = data, aes(x = time, y = !!rlang::sym(parameter), color = depth)) +
    facet_grid(season ~ ., scales = "fixed") +
    scale_color_manual(values = c("indianred", "steelblue")) +
    geom_path(size = 1) +
    geom_point(size = 2) +
    geom_pointrange(aes(ymin = !!rlang::sym(parameter) - se, ymax = !!rlang::sym(parameter) + se), size = 0.3) +
    scale_y_continuous(breaks = seq(0, 0.15, 0.05), limits = c(0, 0.15), expand = c(0, 0)) +
    scale_x_continuous(breaks = seq(0, 50, 10), limits = c(0, 50), expand = c(0, 0)) +
    xlab("Time (h)") + ylab("nM/day") +
    ggtitle(paste(parameter, "- All Seasons")) +
    common_theme()
}

# Function to create sediment plots
create_sediment_plot <- function(data, title_suffix) {
  ggplot(data = data, aes(x = value, y = depth, color = season)) +
    facet_grid(season ~ ., scales = "free_x") +
    geom_point(size = 3, alpha = 0.7) +
    geom_path(size = 1, alpha = 0.7) +
    scale_x_continuous(breaks = seq(0, 20, 4), limits = c(0, 20), expand = c(0, 0)) +
    scale_y_continuous(breaks = seq(-20, 0, 5), limits = c(-20, 0.5), expand = c(0, 0)) +
    xlab("Conc. (nmol/mL)") + ylab("Depth (cm)") +
    ggtitle(paste("Sediment Conc. -", title_suffix)) +
    common_theme()
}

# Create subsets for different seasons
winter <- wzdatach4[wzdatach4$season == "winter", ]
spring <- wzdatach4[wzdatach4$season == "spring", ]
summer <- wzdatach4[wzdatach4$season == "summer", ]
autumn <- wzdatach4[wzdatach4$season == "autumn", ]

# Create CH4 plots for all seasons
ch4_all <- create_ch4_plot(wzdatach4, "All Seasons")
ch4_cold <- create_ch4_plot(winter, "Cold Seasons")
ch4_warm <- create_ch4_plot(wzdatach4[wzdatach4$season %in% c("spring", "summer"), ], "Warm Seasons")

# Combine CH4 plots
ch4_combined <- ggpubr::ggarrange(ch4_cold, ch4_warm, ncol = 1, common.legend = TRUE)

# Create mox and k plots
mox_plot <- create_mox_k_plot(wzdatach4, "mox")
k_plot <- create_mox_k_plot(wzdatach4, "k")

# Combine mox and k plots
mox_k_combined <- ggpubr::ggarrange(mox_plot, k_plot, ncol = 1, common.legend = TRUE)

# Create sediment plots for all seasons
sediment_all <- create_sediment_plot(sed, "All Seasons")
sediment_winter <- create_sediment_plot(winter, "Winter")
sediment_spring <- create_sediment_plot(spring, "Spring")
sediment_summer <- create_sediment_plot(summer, "Summer")
sediment_autumn <- create_sediment_plot(autumn, "Autumn")

# Combine sediment plots
sediment_combined <- ggpubr::ggarrange(sediment_winter, sediment_spring, sediment_summer, sediment_autumn,
                                       ncol = 1, align = "hv", common.legend = FALSE)

############### Methane flux atmosphere analysis ########################

# Data import and format "Probe" column
wzflux <- read.csv("/export/lv6/user/tdegroot/Waddensea/wzflux.csv", header = TRUE, sep = ",")
data <- reshape2::melt(data = wzflux[, c(1:5)], id.vars = c("season", "time"))

# Create plots for flux, wind, and u10
create_flux_wind_u10_plot <- function(data, variable, col_var) {
  ggplot(data = data[data$variable == variable, ], aes(x = time, y = value, color = col_var)) +
    facet_grid(season ~ variable) +
    scale_color_manual(values = c(col_var)) +
    geom_line(size = 1, alpha = 0.75) +
    geom_point(size = 3, alpha = 0.5) +
    scale_x_continuous(breaks = c(0, 10, 20, 30, 40, 50), limits = c(-1, 51), expand = c(0, 0)) +
    scale_y_continuous(breaks = c(0, 100, 200, 300, 400, 500), limits = c(0, 525), expand = c(0, 0)) +
    xlab("Time (h)") + ylab("Value") +
    common_theme()
}

flux <- create_flux_wind_u10_plot(data, "flux", "darkgoldenrod")
wind <- create_flux_wind_u10_plot(data, "wind", "seagreen3")
u10 <- create_flux_wind_u10_plot(data, "u10", "darkmagenta")

# Combine all graphs in one figure
fluxALL <- ggpubr::ggarrange(create_flux_wind_u10_plot(win, "flux", "darkgoldenrod"),
                             create_flux_wind_u10_plot(spr, "flux", "darkgoldenrod"),
                             create_flux_wind_u10_plot(sum, "flux", "darkgoldenrod"),
                             create_flux_wind_u10_plot(aut, "flux", "darkgoldenrod"),
                             ncol = 1, align = "hv", common.legend = FALSE)

windALL <- ggpubr::ggarrange(create_flux_wind_u10_plot(win, "wind", "seagreen3"),
                             create_flux_wind_u10_plot(spr, "wind", "seagreen3"),
                             create_flux_wind_u10_plot(sum, "wind", "seagreen3"),
                             create_flux_wind_u10_plot(aut, "wind", "seagreen3"),
                             ncol = 1, align = "hv", common.legend = FALSE)

ggpubr::ggarrange(fluxALL, windALL, align = "hv", ncol = 2)

# Read environmental data for PCA
env <- wzdatach4[, c("season", "depth", "ch4", "kvalue", "mox", "dens", "temp", "sal")]
penv <- env[, c("ch4", "kvalue", "mox", "dens", "temp", "sal")]
sc.envp <- scale(penv, center = TRUE, scale = TRUE)

# Identify and replace NAs
sc.envp[is.na(sc.envp)] <- apply(sc.envp, 2, function(vec) mean(vec, na.rm = TRUE))

# PCA analysis
pcaenvp <- PCA(sc.envp, scale.unit = TRUE, ncp = 5, graph = TRUE)

# Plots for PCA
scree_plot <- fviz_eig(pcaenvp, addlabels = TRUE, ylim = c(0, 50))
variables_plot <- fviz_pca_var(pcaenvp, col.var = "cos2",
                               gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
                               repel = TRUE)
biplot_plot <- fviz_pca_biplot(pcaenvp, col.ind = env$season, palette = "Spectral", 
                               addEllipses = TRUE, label = "var", col.var = "black", 
                               repel = TRUE, legend.title = "Seasons")

# Print the plots
print(ch4_all)
print(ch4_combined)
print(mox_k_combined)
print(sediment_all)
print(sediment_combined)
print(flux)
print(wind)
print(u10)
print(scree_plot)
print(variables_plot)
print(biplot_plot)