# Load necessary libraries
library(readxl)
library(ggplot2)
library(stats)
library(tidyr)

# Load the data from the Excel file
max_temp_data <- read_excel("Desktop/WalesClimateData.xlsx", sheet = "MaxTemp")
min_temp_data <- read_excel("Desktop/WalesClimateData.xlsx", sheet = "MinTemp")
rainfall_data <- read_excel("Desktop/WalesClimateData.xlsx", sheet = "Rainfall")

# 1: Plot average seasonal maximum temperatures from 1892 to 2022 for summer

# Aggregate data to get average summer temperatures
average_sum_temp <- aggregate(sum ~ year, data = max_temp_data, mean)

# Select decades for plotting
decade_years <- seq(1892, 2022, by = 10)
selected_years <- average_sum_temp$year %in% decade_years
selected_data <- average_sum_temp[selected_years, ]

# Create the plot
ggplot(selected_data, aes(x = year, y = sum)) +
  geom_line() +
  labs(
    title = "Average Seasonal Summer Temperatures (1892-2022)",
    x = "Year",
    y = "Average Summer Temperature"
  ) +
  theme_minimal()

# 2. Calculate correlation between total annual rainfall and annual maximum temperature

# Ensure common years for both datasets
common_years <- intersect(rainfall_data$year, max_temp_data$year)

# Calculate correlation
correlation <- cor(
  rainfall_data$ann[rainfall_data$year %in% common_years],
  max_temp_data$ann[max_temp_data$year %in% common_years]
)
print(correlation)

# 3: Perform t-test for summer minimum temperatures between two periods

# Extract data for t-test
t_test_result <- t.test(
  min_temp_data$aug[min_temp_data$year >= 1890 & min_temp_data$year <= 1910],
  min_temp_data$aug[min_temp_data$year >= 2000 & min_temp_data$year <= 2020]
)
print(t_test_result)

# 4. Perform ANOVA test for annual numbers of wet days in each complete decade

# Check variability within each group
rainfall_data$decade <- cut(rainfall_data$year, breaks = seq(1890, 2020, by = 10))
variability_check <- rainfall_data %>%
  group_by(decade) %>%
  summarize(variability = sd(ann)) %>%
  arrange(desc(variability))

# Print the variability check
print(variability_check)

# Perform ANOVA only if there is enough variability
if (max(variability_check$variability) > 0) {
  rainfall_data$decade <- as.factor(rainfall_data$decade)  # Convert to factor
  anova_result <- oneway.test(ann ~ decade, data = rainfall_data)
  
  # Print ANOVA summary
  print(anova_result)
  
  # Post hoc tests if ANOVA is significant
  if (anova_result$p.value < 0.05) {
    post_hoc <- TukeyHSD(anova_result)
    print(post_hoc)
  }
} else {
  print("Not enough variability within groups for ANOVA.")
}