library(readxl)
library(dplyr)
library(tidyr)
library(lubridate)
library(ggplot2)


try({
  if (requireNamespace("rstudioapi", quietly = TRUE)) {
    this_file <- rstudioapi::getActiveDocumentContext()$path
    if (nzchar(this_file)) {
      setwd(dirname(this_file))
    }
  }
}, silent = TRUE)


data_file <- "Port Talbot Hourly AQ.xlsx"


PortTalbotAQ_raw <- read_excel(data_file, skip = 10, col_types = "text")


PortTalbotAQ <- PortTalbotAQ_raw %>%
  rename(
    date = `Date`,
    time = `Time`,
    
    o3  = `Ozone`,
    no  = `Nitric oxide`,
    no2 = `Nitrogen dioxide`,
    nox = `Nitrogen oxides as nitrogen dioxide`,
    so2 = `Sulphur dioxide`,
    co  = `Carbon monoxide`,
    pm10 = `PM10 particulate matter (Hourly measured)`,
    pm25 = `PM2.5 particulate matter (Hourly measured)`,
    
    wind_dir   = `Modelled Wind Direction`,
    wind_speed = `Modelled Wind Speed`,
    temp       = `Modelled Temperature`
  )


PortTalbotAQ <- PortTalbotAQ %>%
  mutate(
    date = as.Date(date),

    time = format(as.POSIXct(time, format = "%H:%M:%S"), "%H:%M:%S"),
    datetime = as.POSIXct(paste(date, time), format = "%Y-%m-%d %H:%M:%S", tz = "UTC")
  )

cat("\n---------------- PortTalbotAQ data frame ----------------\n")
print(PortTalbotAQ)
cat("\nRows:", nrow(PortTalbotAQ), "  Columns:", ncol(PortTalbotAQ), "\n")


if (interactive()) {
  View(PortTalbotAQ)
}


make_numeric <- function(x) {
  x <- as.character(x)
  x <- trimws(x)
  x[x == "No data"] <- NA
  x[x == ""] <- NA
  as.numeric(x)
}

pollutants <- c("o3", "no", "no2", "nox", "so2", "co", "pm10", "pm25")
invalid_percent <- c()

for (p in pollutants) {
  column_values <- as.character(PortTalbotAQ[[p]])
  column_values <- trimws(column_values)
  missing_flags <- is.na(column_values) | column_values == "No data" | column_values == ""
  invalid_percent[p] <- mean(missing_flags) * 100
}

invalidData <- data.frame(
  pollutant = names(invalid_percent),
  percent_no_measurement = as.numeric(invalid_percent),
  row.names = NULL
)

cat("\n---------------- invalidData (% No data) ----------------\n")
print(invalidData)

if (interactive()) {
  View(invalidData)
}



displayPollutantGraph <- function(pollutant_name, start_month, start_year, end_month, end_year) {
  
  if (!(pollutant_name %in% names(PortTalbotAQ))) {
    stop("That pollutant column name is not in the data frame.")
  }
  

  start_date <- as.Date(paste0(start_year, "-", sprintf("%02d", start_month), "-01"))
  end_date <- as.Date(paste0(end_year, "-", sprintf("%02d", end_month), "-01"))
  end_date <- ceiling_date(end_date, "month") - days(1)
  

  df <- PortTalbotAQ %>%
    filter(date >= start_date, date <= end_date) %>%
    filter(day(date) == 1) %>%
    filter(time == "12:00:00") %>%
    mutate(value = make_numeric(.data[[pollutant_name]])) %>%
    filter(!is.na(value)) %>%
    arrange(date)
  

  ggplot(df, aes(x = date, y = value)) +
    geom_line() +
    geom_point() +
    labs(
      title = paste("Pollutant levels at noon on the 1st:", pollutant_name),
      x = "Date",
      y = pollutant_name
    )
}

cat("\n--- Task 3 graph: NOx (as NO2) from July 2015 to Nov 2024 ---\n")
print(displayPollutantGraph("nox", 7, 2015, 11, 2024))



comparePollutants <- function(pollutant_vector) {
  

  for (p in pollutant_vector) {
    if (!(p %in% names(PortTalbotAQ))) {
      stop(paste("Pollutant not found:", p))
    }
  }
  

  df <- PortTalbotAQ %>%
    filter(date >= as.Date("2015-01-01"), date <= as.Date("2024-12-31")) %>%
    filter(day(date) == 1, time == "12:00:00") %>%
    select(date, all_of(pollutant_vector))
  

  df_long <- pivot_longer(df, cols = all_of(pollutant_vector),
                          names_to = "pollutant", values_to = "raw")
  
  df_long <- df_long %>%
    mutate(value = make_numeric(raw)) %>%
    filter(!is.na(value))
  
  ggplot(df_long, aes(x = date, y = value, linetype = pollutant)) +
    geom_line() +
    labs(
      title = "Pollutant comparison (noon on 1st of month, 2015–2024)",
      x = "Date",
      y = "Level",
      linetype = "Pollutant"
    )
}

cat("\n--- Task 4 graph: compare 4 pollutants ---\n")
print(comparePollutants(c("o3", "no2", "pm10", "pm25")))



O3NOCorrelation <- function() {
  
  df <- PortTalbotAQ %>%
    transmute(
      o3_num = make_numeric(o3),
      no2_num = make_numeric(no2)
    ) %>%
    filter(!is.na(o3_num), !is.na(no2_num))
  
  r_value <- cor(df$o3_num, df$no2_num)
  
  cat("\n--- Task 5: Correlation between O3 and NO2 ---\n")
  cat("Pearson correlation =", r_value, "\n")
  
  ggplot(df, aes(x = no2_num, y = o3_num)) +
    geom_point(alpha = 0.4) +
    labs(
      title = "O3 vs NO2 scatterplot",
      subtitle = paste("Correlation =", round(r_value, 4)),
      x = "NO2",
      y = "O3"
    )
}

print(O3NOCorrelation())



df_t <- PortTalbotAQ %>%
  transmute(
    o3_num = make_numeric(o3),
    ws_num = make_numeric(wind_speed)
  ) %>%
  filter(!is.na(o3_num), !is.na(ws_num)) %>%
  mutate(wind_group = ifelse(ws_num < 5, "<5", ">=5"))

t_test_result <- t.test(o3_num ~ wind_group, data = df_t)

cat("\n--- Task 6: t-test result ---\n")
print(t_test_result)


if (t_test_result$p.value < 0.05) {
  cat("Conclusion: p < 0.05 so I reject H0 (there is a significant difference).\n")
} else {
  cat("Conclusion: p >= 0.05 so I do not reject H0 (no significant difference found).\n")
}




df_aov <- PortTalbotAQ %>%
  filter(time == "17:00:00") %>%
  filter(year(date) %in% c(2016, 2018, 2020)) %>%
  filter(month(date) %in% c(11, 12)) %>%
  transmute(
    year = factor(year(date)),
    no2_num = make_numeric(no2)
  ) %>%
  filter(!is.na(no2_num))

anova_model <- aov(no2_num ~ year, data = df_aov)

cat("\n--- Task 7: ANOVA result ---\n")
print(summary(anova_model))


p_value_anova <- summary(anova_model)[[1]][["Pr(>F)"]][1]

if (!is.na(p_value_anova) && p_value_anova < 0.05) {
  cat("ANOVA p < 0.05 so I reject H0. Now I will do Tukey HSD to see which years differ.\n")
  print(TukeyHSD(anova_model))
} else {
  cat("ANOVA p >= 0.05 so I do not reject H0. I will not do post-hoc tests.\n")
}