# port_talbot_aq_cleaned.R
# Clean, analyze and plot Port Talbot Hourly AQ dataset
# - Robust datetime parsing from separate Date + Time columns (dd-mm-yyyy + HH:MM:SS)
# - Strip units / "V ugm-3" / "No data" and convert pollutant & met columns to numeric
# - Compute % invalid (missing) per pollutant
# - Reusable plotting/analysis functions:
#     displayPollutantGraph(), comparePollutants(), O3NOCorrelation()
# - Hypothesis tests:
#     t-test: ozone by wind speed (<5 vs >=5 m/s)
#     ANOVA: NO2 at 17:00 in Nov-Dec for years 2016, 2018, 2020 (+ Tukey HSD)
#
# Author: Shafik Nagori (2437695) 
#Run in R console. If package installation is blocked,
# install packages manually before running.

# ---------------------------
# 0) Install / load packages
# ---------------------------
pkgs <- c("dplyr", "tidyr", "lubridate", "ggplot2", "stringr", "forcats", "readr")
new_pkgs <- pkgs[!(pkgs %in% installed.packages()[, "Package"])]
if (length(new_pkgs)) {
  message("Installing missing packages: ", paste(new_pkgs, collapse = ", "))
  install.packages(new_pkgs, dependencies = TRUE)
}
suppressPackageStartupMessages({
  library(dplyr)
  library(tidyr)
  library(lubridate)
  library(ggplot2)
  library(stringr)
  library(forcats)
  library(readr)
  library(reshape2)
})

# ---------------------------
# 1) Read CSV (safe defaults)
# ---------------------------
data_file <- "Port_Talbot_Hourly_AQ.csv"
if (!file.exists(data_file)) stop("CSV not found at: ", data_file)

raw <- read.csv(data_file,
                header = TRUE,
                stringsAsFactors = FALSE,
                check.names = FALSE,
                na.strings = c("No data", "no data", "", "NA"))

# Quick preview (first few rows)
message("Preview of raw data (first 5 rows):")
print(utils::head(raw, 5))

# ---------------------------
# 2) Clean column names
# ---------------------------
# Convert names to lowercase snake_case to make coding easier.
clean_names <- function(nm) {
  nm %>%
    str_trim() %>%
    str_to_lower() %>%
    # replace punctuation with underscore
    str_replace_all("[^a-z0-9]+", "_") %>%
    # collapse multiple underscores
    str_replace_all("_+", "_") %>%
    str_replace_all("^_|_$", "")
}
names(raw) <- clean_names(names(raw))

# Show cleaned column names
message("Cleaned column names:")
print(names(raw))

# ---------------------------
# 3) Build datetime column
# ---------------------------
# Dataset uses separate date and time columns (example: Date = "01-01-2020", Time = "01:00:00")
# After cleaning names we expect "date" and "time" or similar.
date_col <- intersect(c("date", "day"), names(raw))[1]
time_col <- intersect(c("time", "hour"), names(raw))[1]

if (is.na(date_col) || is.na(time_col)) {
  stop("Could not find separate Date and Time columns. Found: ", paste(names(raw), collapse = ", "))
}

# Combine and parse: data sample shows dd-mm-yyyy and HH:MM:SS
raw <- raw %>%
  mutate(
    # paste preserves original strings; dmy_hms handles "01-01-2020 01:00:00"
    datetime = dmy_hms(paste(.data[[date_col]], .data[[time_col]]), tz = "UTC")
  )

if (all(is.na(raw$datetime))) {
  stop("Failed to parse datetime — check Date and Time formats in the CSV. Expected dd-mm-yyyy and HH:MM:SS.")
}

# Optional: reorder so datetime is first
raw <- raw %>% select(datetime, everything())

# ---------------------------
# 4) Clean numeric columns
# ---------------------------
# Many pollutant and met columns are strings like "32.39686 V ugm-3" or "No data" or
# "25.1 V ugm-3 (TEOM FDMS)". We extract the first numeric token for each column we care about.
# Define the column name candidates (based on your dataset header preview)
candidates <- c(
  "ozone",
  "nitric_oxide",             # sometimes "nitric.oxide" -> cleaned name nitric_oxide
  "nitrogen_dioxide",
  "nitrogen_oxides_as_nitrogen_dioxide",
  "sulphur_dioxide",
  "carbon_monoxide",
  "pm10_particulate_matter_hourly_measured",
  "pm2_5_particulate_matter_hourly_measured",
  "modelled_wind_direction",
  "modelled_wind_speed",
  "modelled_temperature"
)

# Because names changed by cleaning, try to find the actual columns present using contains
find_existing <- function(patterns, nm_vec) {
  found <- character(0)
  for (p in patterns) {
    matches <- nm_vec[str_detect(nm_vec, fixed(p))]
    if (length(matches) > 0) found <- c(found, matches[1])
  }
  unique(found)
}

existing_cols <- find_existing(candidates, names(raw))
message("Detected columns to clean (candidates matched):")
print(existing_cols)

# Function to extract the first number (handles integers, decimals, negative)
extract_num <- function(x) {
  # Replace common non-breaking spaces or weird whitespace
  x2 <- ifelse(is.na(x), NA_character_, str_replace_all(as.character(x), "\\u00A0", " "))
  # Use readr::parse_number which is robust (handles commas)
  suppressWarnings(parse_number(x2, locale = locale(decimal_mark = ".")))
}

# Create numeric versions with suffix "_num"
for (col in existing_cols) {
  newcol <- paste0(col, "_num")
  raw[[newcol]] <- extract_num(raw[[col]])
}

# Show a quick numeric preview
message("Preview numeric columns (first 6 rows):")
num_preview <- raw %>% select(datetime, ends_with("_num")) %>% head(6)
print(num_preview)

# ---------------------------
# 5) Calculate invalid (missing) percentages
# ---------------------------
# Use the numeric versions to compute percent missing.
numeric_cols <- names(raw)[str_ends(names(raw), "_num")]
invalidData <- data.frame(
  pollutant = numeric_cols,
  percent_missing = sapply(raw[numeric_cols], function(x) mean(is.na(x)) * 100),
  row.names = NULL,
  stringsAsFactors = FALSE
)
message("Percent missing (using numeric columns):")
print(invalidData)

# ---------------------------
# 6) Helper: sequence of first-of-month at noon
# ---------------------------
first_of_month_noon_seq <- function(start_month, start_year, end_month, end_year) {
  start <- make_date(start_year, start_month, 1)
  end   <- make_date(end_year, end_month, 1)
  seq(from = as.POSIXct(start, tz = "UTC") + hours(12),
      to   = as.POSIXct(end, tz = "UTC") + hours(12),
      by = "1 month")
}

# ---------------------------
# 7) Reusable plotting/analysis functions
# ---------------------------

# 7a) displayPollutantGraph:
#    Time-series for pollutant at noon on 1st of each month between start and end.
#    'pollutant' should be the original cleaned column name (e.g., "ozone_num" or "no2_num").
displayPollutantGraph <- function(df, pollutant_num_col, start_month, start_year, end_month, end_year) {
  if (!(pollutant_num_col %in% names(df))) stop("Pollutant numeric column not found: ", pollutant_num_col)
  times <- first_of_month_noon_seq(start_month, start_year, end_month, end_year)
  plot_df <- df %>% filter(datetime %in% times) %>% select(datetime, value = all_of(pollutant_num_col))
  if (nrow(plot_df) == 0) stop("No observations matched the first-of-month-noon filter for that range.")
  p <- ggplot(plot_df, aes(x = datetime, y = value)) +
    geom_line() +
    geom_point() +
    labs(title = paste("Monthly (1st @ noon) levels:", pollutant_num_col),
         x = "Month (1st @ 12:00)", y = "Value") +
    theme_minimal()
  print(p)
  invisible(p)
}

# 7b) comparePollutants:
#    Overlay multiple pollutant numeric columns (vector of names) at first-of-month noon points.
comparePollutants <- function(df, pollutant_num_cols, start_month = NULL, start_year = NULL, end_month = NULL, end_year = NULL) {
  pollutant_num_cols <- intersect(pollutant_num_cols, names(df))
  if (length(pollutant_num_cols) < 1) stop("No pollutant numeric columns found to compare.")
  if (!is.null(start_month)) {
    times <- first_of_month_noon_seq(start_month, start_year, end_month, end_year)
    df_sub <- df %>% filter(datetime %in% times)
  } else {
    df_sub <- df %>% filter(day(datetime) == 1, hour(datetime) == 12, minute(datetime) %in% c(0, NA))
  }
  plot_df <- df_sub %>%
    select(datetime, all_of(pollutant_num_cols)) %>%
    pivot_longer(cols = all_of(pollutant_num_cols), names_to = "pollutant", values_to = "value")
  if (nrow(plot_df) == 0) stop("No rows available for chosen filter.")
  p <- ggplot(plot_df, aes(x = datetime, y = value, color = pollutant)) +
    geom_line() + geom_point() +
    labs(title = "Comparison of pollutants (1st of month @ noon)",
         x = "Month", y = "Value") +
    theme_minimal()
  print(p)
  invisible(p)
}

# 7c) O3NOCorrelation:
#    Correlation test (Pearson) and scatterplot between ozone and nitrogen dioxide numeric columns.
O3NOCorrelation <- function(df, o3_col = "ozone_num", no2_col = "nitrogen_dioxide_num") {
  if (!(o3_col %in% names(df))) stop("Ozone numeric column missing:", o3_col)
  if (!(no2_col %in% names(df))) stop("NO2 numeric column missing:", no2_col)
  cor_df <- df %>% filter(!is.na(.data[[o3_col]]), !is.na(.data[[no2_col]]))
  if (nrow(cor_df) < 3) stop("Not enough complete observations for correlation test.")
  ct <- cor.test(cor_df[[o3_col]], cor_df[[no2_col]], use = "complete.obs", method = "pearson")
  print(ct)
  p <- ggplot(cor_df, aes_string(x = no2_col, y = o3_col)) +
    geom_point(alpha = 0.3) +
    geom_smooth(method = "lm", se = TRUE) +
    labs(title = paste("Scatterplot & fit:", o3_col, "vs", no2_col),
         x = no2_col, y = o3_col) +
    theme_minimal()
  print(p)
  invisible(list(plot = p, cor_test = ct))
}

# ---------------------------
# 8) Hypothesis tests
# ---------------------------

# 8a) t-test: Ozone levels vs wind speed (<5 vs >=5 m/s)
# Identify ozone and wind-speed numeric columns
ozone_col <- ifelse("ozone_num" %in% names(raw), "ozone_num", NA)
wind_col  <- ifelse("modelled_wind_speed_num" %in% names(raw), "modelled_wind_speed_num", NA)

if (!is.na(ozone_col) && !is.na(wind_col)) {
  tt_df <- raw %>%
    select(datetime, ozone = .data[[ozone_col]], wind = .data[[wind_col]]) %>%
    filter(!is.na(ozone), !is.na(wind)) %>%
    mutate(wind_group = ifelse(wind < 5, "<5", ">=5"))
  # Ensure both groups exist
  if (length(unique(tt_df$wind_group)) == 2 && nrow(tt_df) >= 10) {
    # Check variance equality and run t-test (Welch's or pooled depending on var.test)
    var_test <- var.test(ozone ~ wind_group, data = tt_df)
    message("Variance test p-value: ", signif(var_test$p.value, 4))
    # Use Welch two-sample t-test by default (robust)
    t_result <- t.test(ozone ~ wind_group, data = tt_df)
    message("T-test (Ozone by wind_group):")
    print(t_result)
  } else {
    message("Skipping t-test: need at least two wind groups '<5' and '>=5' with sufficient observations.")
  }
} else {
  message("Skipping t-test: ozone or wind speed numeric column missing.")
}

# 8b) ANOVA: NO2 across Nov-Dec for years 2016, 2018, 2020 at 17:00
# Identify NO2 numeric column
no2_col <- ifelse("nitrogen_dioxide_num" %in% names(raw), "nitrogen_dioxide_num", NA)
if (!is.na(no2_col)) {
  anova_df <- raw %>%
    filter(year(datetime) %in% c(2016, 2018, 2020),
           month(datetime) %in% c(11, 12),
           hour(datetime) == 17,
           minute(datetime) %in% c(0, NA)) %>%
    select(datetime, no2 = .data[[no2_col]]) %>%
    mutate(year = factor(year(datetime))) %>%
    filter(!is.na(no2))
  message("ANOVA subset rows:", nrow(anova_df))
  if (nrow(anova_df) >= 3 && length(unique(anova_df$year)) >= 2) {
    fit <- aov(no2 ~ year, data = anova_df)
    message("ANOVA summary:")
    print(summary(fit))
    # Check significance and perform Tukey HSD post-hoc if appropriate
    pval <- summary(fit)[[1]][["Pr(>F)"]][1]
    if (!is.null(pval) && pval < 0.05) {
      message("ANOVA significant (p < 0.05). Running TukeyHSD post-hoc:")
      tryCatch({
        tuk <- TukeyHSD(fit)
        print(tuk)
      }, error = function(e) {
        message("TukeyHSD failed: ", e$message)
      })
    } else {
      message("ANOVA not significant (p >= 0.05) — skipping post-hoc.")
    }
  } else {
    message("Skipping ANOVA: not enough observations in requested subset.")
  }
} else {
  message("Skipping ANOVA: NO2 numeric column not found.")
}