# Install the readxl,dplyr,tidyr,ggplot2 packages if not already installed
# Load the  packages
library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)

# ---------------------------------Read data function ---------------------------------
read_excel_sheet <- function(file, sheet_name) {
  readxl::read_excel(file, sheet = sheet_name)
}

# --------------------------------- Data Pre-Processing Function---------------------------------
# Handle Missing Data and Data Rounding
handle_missing_data <- function(data) {
  data[data == '---'] <- NA # Replace '---' to NAs
  data <- as.numeric(data) # Convert character to Number
  data[is.na(data)] <- mean(data, na.rm = TRUE) # Replace NAs with Mean
  data <- format(round(as.numeric(data), 2), nsmall = 2) # Round the Data
  data <- as.numeric(data) # Convert Character to Number
  return(data)
}

# --------------------------------- Question:1 Functions---------------------------------

# Function to calculate average seasonal maximum temperatures
calculate_average_seasonal_max <- function(data) {
  average_seasonal_max <- data %>%
    group_by(year) %>%
    summarize(
      spr = mean(spr, na.rm = TRUE),
      sum = mean(sum, na.rm = TRUE),
      aut = mean(aut, na.rm = TRUE),
      win = mean(win, na.rm = TRUE)
    )
  return(average_seasonal_max)
}

# Function to extract every tenth year for plotting
extract_every_tenth_year <- function(data) {
  #every_tenth_year <- data[seq(1, nrow(data), by = 10), ] # From 1884
  every_tenth_year <- data[seq(9, nrow(data), by = 10), ] # From 1892
  return(every_tenth_year)
}

# Function to create a multiple bar graph using ggplot2()
create_multiple_bar_graph <- function(data) {
  plot_data <- gather(data, key = "Season", value = "Temperature", -year)
  
  multiple_bar_plot <- ggplot(plot_data, aes(x = factor(year), y = Temperature, fill = Season)) +
    geom_bar(stat = "identity", position = "dodge", color = "white") +
    labs(title = "Multiple Bar Graph of Average Seasonal Maximum Temperatures",
         x = "Year", y = "Temperature") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
  
  print(multiple_bar_plot)
}

# --------------------------------- Question:2 Functions---------------------------------

# Function to merge and subset dataframes, and create a scatterplot with correlation
analyze_annual_max_temp_vs_rainfall <- function(max_temp_data, rainfall_mm_data) {
  # Merge dataframes by year and include only complete observations
  merged_data <- merge(max_temp_data, rainfall_mm_data, by = "year", all = FALSE)
  
  # Calculate and print the correlation coefficient
  calculate_and_print_correlation(merged_data)
  
  # Create a scatterplot with color and trendline
  plot(
    merged_data$ann.x,
    merged_data$ann.y,
    xlab = "Annual Maximum Temperature",
    ylab = "Annual Rainfall",
    main = "Scatterplot of Annual Maximum Temperature vs. Annual Rainfall",
    col = "blue",    # Set color to blue
    pch = 10,        # Set point shape to solid circle
  )
  # Add a trendline to the scatterplot
  abline(lm(ann.y ~ ann.x, data = merged_data), col = "red", lwd = 2)
}

# Function to calculate and print the correlation coefficient
calculate_and_print_correlation <- function(data) {
  correlation <- cor(data$ann.x, data$ann.y, use = "complete.obs")
  print(paste("Correlation coefficient:", round(correlation, 4)))
}

# --------------------------------- Question:3 Functions---------------------------------

# Function to subset data for specific time periods
subset_temp_data <- function(data, start_year, end_year, month) {
  subset_data <- data[data$year >= start_year & data$year <= end_year, ]
  return(subset_data[[month]])
}

# Function to perform Shapiro test and t-test on two data samples
perform_t_test <- function(data1, data2) {
  # Perform Shapiro test
  shapiro_result1 <- shapiro.test(data1)
  shapiro_result2 <- shapiro.test(data2)
  
  # Print Shapiro test results
  print(paste("Shapiro Test p-value (1890-1910):", shapiro_result1$p.value)) # 0.0522694417128314
  print(paste("Shapiro Test p-value (2000-2020):", shapiro_result2$p.value)) # 0.136542565918646
  
  # Perform t-test
  t_test_result <- t.test(data1, data2)
  
  # Print t-test results
  print(t_test_result)
}

# --------------------------------- Question:4 Function---------------------------------

# Function to perform ANOVA and post-hoc tests on annual numbers of wet days
perform_anova_and_posthoc <- function(data, start_year, end_year, decade_breaks) {
  # Create a new variable 'Decade' to represent each complete decade
  data$Decade <- cut(data$year, breaks = decade_breaks, labels = FALSE)
  
  # Perform ANOVA
  anova_result <- aov(ann ~ as.factor(Decade), data = data)
  
  # Print ANOVA results
  print(summary(anova_result))
  
  #To identify where the differences lie,we perform post-hoc tests. 
  # Perform post-hoc tests (Tukey's HSD)
  posthoc_result <- TukeyHSD(anova_result)
  
  # Print post-hoc test results
  print(posthoc_result)
  
  # Create a box plot for visual representation
  boxplot(ann ~ as.factor(Decade), data = data,
          main = "Box Plot of Annual Numbers of Wet Days by Decade",
          xlab = "Decade",
          ylab = "Annual Numbers of Wet Days",
          col = "lightblue",
          border = "black")
}

# ******************------------------  Main Script ------------------******************
# ******************------------------  Main Script ------------------******************
# ******************------------------  Main Script ------------------******************

# Store the file name
walesclimate <- "Wales Climate Data.xlsx"
# Get data from different sheets
max_temp_data <- read_excel_sheet(walesclimate, "Max Temp")
min_temp_data <- read_excel_sheet(walesclimate, "Min Temp")
sunshine_data <- read_excel_sheet(walesclimate, "Sunshine")
rainfall_mm_data <- read_excel_sheet(walesclimate, "Rainfall")
rainfall_days_data <- read_excel_sheet(walesclimate, "Rain Days")

# Data Pre-Processing
max_temp_data$win <- handle_missing_data(max_temp_data$win)
min_temp_data$win <- handle_missing_data(min_temp_data$win)
sunshine_data$win <- handle_missing_data(sunshine_data$win)
rainfall_mm_data$win <- handle_missing_data(rainfall_mm_data$win)
rainfall_days_data$win <- handle_missing_data(rainfall_days_data$win)

#--------------------------------------------------------------------------------------------------------
# *** Q1. Produce an appropriate graphic to show how the average seasonal maximum temperatures 
# *** have changed from 1892 to 2022. You should plot the values for every tenth year rather than 
# *** for every year.

# Calculate average seasonal maximum temperatures
average_seasonal_max <- calculate_average_seasonal_max(max_temp_data)
  
# Extract every tenth year for plotting
every_tenth_year <- extract_every_tenth_year(average_seasonal_max)
  
# Create a multiple bar graph using ggplot2
create_multiple_bar_graph(every_tenth_year)

#--------------------------------------------------------------------------------------------------------
# *** Q2. Calculate the correlation between the total annual rainfall and the annual maximum temperature.

# Print the length of both dataframe
print(length(max_temp_data$ann)) #139
print(length(rainfall_mm_data$ann)) #187
# Both dataframe has different number of records

# Analyze the relationship between annual maximum temperature and rainfall
analyze_annual_max_temp_vs_rainfall(max_temp_data, rainfall_mm_data) # Correlation coefficient: 0.0725"

#--------------------------------------------------------------------------------------------------------
# *** Q3. Perform a t-test to identify whether there is a significant difference in minimum temperatures 
# *** between August 1890-1910 and August 2000-2020.

# Subset the data for the two time periods
min_temp_1890_1910 <- subset_temp_data(min_temp_data, 1890, 1910, "aug")
min_temp_2000_2020 <- subset_temp_data(min_temp_data, 2000, 2020, "aug")
  
# Boxplot for minimum temperature of August between 1890 and 1910
boxplot(min_temp_1890_1910, main = "Minimum Temperature of August between 1890 and 1910")
# No Outliers. No Need To apply filter on data

# Boxplot for minimum temperature of August between 2000 and 2020
boxplot(min_temp_2000_2020, main = "Minimum Temperature of August between 2000 and 2020")
# No Outliers. No Need To apply filter on data
# Perform Shapiro test and t-test
perform_t_test(min_temp_1890_1910, min_temp_2000_2020)
  

#--------------------------------------------------------------------------------------------------------
# *** Q4. Perform an ANOVA test to identify if there are statistically significant differences in the 
# *** annual numbers of wet days in each complete decade starting in 1891 and finishing in 2020. 
# *** If such a difference is found, perform the appropriate ad-hoc tests to identify where the 
# *** differences lie.

# Define decade breaks for ANOVA
decade_breaks <- seq(1890, 2020, by = 10)
  
# Perform ANOVA and post-hoc tests
perform_anova_and_posthoc(rainfall_days_data, 1891, 2020, decade_breaks)
# The ANOVA result suggests that there is a marginally significant difference in 
# The mean annual numbers of wet days across the different decades (p-value = 0.0917).