#Dev Patel- 2434793

# ILO1 – I am handling data correctly (import, rename, clean,type conversion, subsetting, NA handling).
# ILO2 – I am communicating results with clear visualisations (scatter, heatmap, boxplot, axes, labels, legends).
# ILO3 – I am writing R code using functions, loops, indexing, apply, simple logic.
# ILO4 – I am doing correct statistical reasoning and tests (correlation, t-test, ANOVA, Tukey, decisions)

# Q1 – IMPORT + RENAME + SHOW (ILO1)
# I am loading the CSV with headers.
# I am renaming long headers to short easy names.
# I am showing the first rows to check that it is working or not.

# I am reading the file (keeping exact header text so I can rename myself)
PortTalbotAQ <- read.csv("Port_Talbot_Hourly_AQ.csv", header = TRUE, stringsAsFactors = FALSE, check.names = FALSE) # ILO1: import with headers

# I am printing original names and a preview, just to check
names(PortTalbotAQ)  # ILO1: check column names
head(PortTalbotAQ) # ILO1: quick preview

# ILO1: I am renaming columns to short names I can type easily later
names(PortTalbotAQ)[names(PortTalbotAQ)=="Date"] <- "date"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Time"] <- "time"

names(PortTalbotAQ)[names(PortTalbotAQ)=="Ozone"] <- "o3"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Nitric oxide"] <- "no"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Nitrogen dioxide"] <- "no2"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Nitrogen oxides as nitrogen dioxide"] <- "no_no2"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Sulphur dioxide"] <- "so2"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Carbon monoxide"] <- "co"
names(PortTalbotAQ)[names(PortTalbotAQ)=="PM10 particulate matter (Hourly measured)"] <- "pm10"
names(PortTalbotAQ)[names(PortTalbotAQ)=="PM2.5 particulate matter (Hourly measured)"] <- "pm25"

# I am also renaming modelled met fields that I will need later (Q6 needs wind speed)
names(PortTalbotAQ)[names(PortTalbotAQ)=="Modelled Wind Speed"] <- "wind_speed"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Modelled Wind Direction"] <- "wind_direction"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Modelled Temperature"] <- "temp"

# ILO1: I am printing the cleaned names and a small preview
names(PortTalbotAQ)
head(PortTalbotAQ)

# Q2 – invalidData: % missing per pollutant (ILO1, ILO2)
# I am counting how many readings are "No data" or empty or NA for each pollutant.
# I am building a small data frame called invalidData with those percentages.

# I am listing only pollutant columns for this task
pollutants <- c("o3","no","no2","no_no2","so2","co","pm10","pm25")

# I am creating empty boxes to store names and percentages
pollutant_name <- NULL
empty_percentages <- NULL

# I am looping through columns and counting % "No data"
for (p in pollutants) {          # ILO1: basic loop over columns
  vals <- as.character(PortTalbotAQ[[p]])   # ILO1: I convert to text to check "No data"
  total_count <- length(vals)
  total_empty_data <- 0

  for (i in 1:total_count) {
    if (is.na(vals[i]) || vals[i]=="" || vals[i]=="No data") { # ILO1: identifying invalid entries
      total_empty_data <- total_empty_data + 1
    }
  }
  percentage_empty <- (total_empty_data / total_count) * 100  # ILO2: summarising missing data into % measure
  pollutant_name <- c(pollutant_name, p)
  empty_percentages      <- c(empty_percentages, percentage_empty)
}

# ILO2: I am making the invalidData table now
invalidData <- data.frame(pollutant = pollutant_name, percentages_of_no_data = empty_percentages)

# I am showing invalidData 
invalidData

# I am also cleaning the pollutant and met columns now for later stats:
# I change "No data"/"" to NA and convert to numeric.
clean_columns <- c(pollutants, "wind_speed", "wind_direction", "temp")
for (p in clean_columns) {
  if (!is.null(PortTalbotAQ[[p]])) {
    x <- as.character(PortTalbotAQ[[p]]) #ILO1: cleaning invalid textual entries
    for (i in 1:length(x)) {
      if (x[i]=="" || x[i]=="No data") { x[i] <- NA }
    }
    PortTalbotAQ[[p]] <- as.numeric(x)  # ILO1: wind_dir is in degrees; numeric is fine
  }
}

# I am creating date/time pieces for filtering (used in Q3, Q4, Q7)
# My date format is "DD-MM-YYYY" and time looks like "H:MM:SS" or "HH:MM:SS".
# I am making year/month/day with substr because dates are fixed-width
PortTalbotAQ$year  <- as.numeric(substr(PortTalbotAQ$date, 7, 10))
PortTalbotAQ$month <- as.numeric(substr(PortTalbotAQ$date, 4, 5))
PortTalbotAQ$day   <- as.numeric(substr(PortTalbotAQ$date, 1, 2))

# I am making hour using a loop so it also works for 1-digit hours like "1:00:00"
time_vector <- as.character(PortTalbotAQ$time)
hour_value <- numeric(length(time_vector))

for (i in 1:length(time_vector)) {
  if (is.na(time_vector[i]) || time_vector[i] == "") {
    hour_value[i] <- NA
  } else {
    parts <- strsplit(time_vector[i], ":")[[1]]   # e.g. c("1","00","00")
    # I am taking the first part before ":" and converting to number
    hour_value[i] <- as.numeric(parts[1])
  }
}
PortTalbotAQ$hour <- hour_value

# I am doing quick checks
table(is.na(PortTalbotAQ$year))
table(is.na(PortTalbotAQ$month))
table(is.na(PortTalbotAQ$day))
table(is.na(PortTalbotAQ$hour))

names(PortTalbotAQ)[names(PortTalbotAQ)=="Ozone"] <- "o3"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Nitric oxide"] <- "no"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Nitrogen dioxide"] <- "no2"
names(PortTalbotAQ)[names(PortTalbotAQ)=="PM2.5 particulate matter (Hourly measured)"] <- "pm25"
names(PortTalbotAQ)[names(PortTalbotAQ)=="PM10 particulate matter (Hourly measured)"] <- "pm10"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Sulphur dioxide"] <- "so2"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Carbon monoxide"] <- "co"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Modelled Wind Speed"] <- "wind_speed"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Modelled Wind Direction"] <- "wind_direction"
names(PortTalbotAQ)[names(PortTalbotAQ)=="Air Temperature"] <- "air_temp"

# Q3 – displayPollutantGraph() (ILO2, ILO3)
# I am plotting the pollutant at NOON on the 1st of every month in a time window.
# I am using a COLOURED SCATTER plot.
# nitrogen oxides as nitrogen dioxide (no_no2), Jul 2015 to Nov 2024.
displayPollutantGraph <- function(pollutantName, startMonth, startYear, endMonth, endYear) {
  
  # I am keeping only monthly noon rows (1st day at 12:00)
  monthly_noon_rows <- PortTalbotAQ[PortTalbotAQ$day==1 & PortTalbotAQ$hour==12, ] #ILO2: filtering time/location subset
  
  # I am marking rows in the required date range
  selected_range <- rep(FALSE, nrow(monthly_noon_rows))
  for (i in 1:nrow(monthly_noon_rows)) {
    this_year <- monthly_noon_rows$year[i]; 
    this_month <- monthly_noon_rows$month[i]
    after_start <- (this_year > startYear) || (this_year==startYear && this_month>=startMonth)
    before_end  <- (this_year < endYear)   || (this_year==endYear   && this_month<=endMonth)
    if (after_start && before_end) { 
      selected_range[i] <- TRUE 
    }
  }
  filtered_data <- monthly_noon_rows[selected_range, ]
  
  # I am getting the pollutant values and removing NA
  y_values <- filtered_data[[pollutantName]]  # ILO2: selecting a specific pollutant series
  valid_index   <- which(!is.na(y_values))
  y_values <- y_values[valid_index]
  
  # I am creating x-axis positions and labels
  x_values   <- 1:length(y_values) # ILO3: preparing visual x-axis index
  x_labels <- paste(filtered_data$year[valid_index], sprintf("%02d", filtered_data$month[valid_index]), "01 12:00", sep="-")
  
  # Background color panel
  par(bg = "white")                   # Window background
  
  # ILO3: visualising pollutant trend (scatter)
  plot(x_values, y_values,
       pch = 19,                      # solid round points
       col = "steelblue",             # point color 
       xlab = "Month (1st day at 12:00)",
       ylab = pollutantName,
       main = paste("Monthly noon levels for", pollutantName),
       cex = 1.0)
  
  # ILO3: improving visual readability
  # I added lighter axis text, rotated
  axis(1, at=x_values, labels=x_labels, las=2, cex.axis=0.55, col.axis="gray40")
  # I added horizontal light grid lines to help comparison
  abline(h = pretty(y_values), col="gray85", lty=3)
}

# I am calling function to display graph
displayPollutantGraph("no_no2", 7, 2015, 11, 2024)

# Q4 – comparePollutants() (ILO2, ILO3)
# I am comparing several pollutants together at monthly noon readings (all 10 yrs).
# I am using a HEATMAP (base R image())
# Rows are months across the 10 years; columns are chosen pollutants.
# I am min–max scaling each column to 0..1 so colors are comparable.

comparePollutants <- function(pollutantVector) {
  
  monthly_noon_rows <- PortTalbotAQ[PortTalbotAQ$day==1 & PortTalbotAQ$hour==12, ]
  
  total_months <- nrow(monthly_noon_rows)
  total_pol  <- length(pollutantVector)
  scaled_matrix <- matrix(NA_real_, nrow=total_months, ncol=total_pol) #ILO2: building structured comparison matrix
  
  for (poll_index in 1:total_pol) {
    pollutant_values <- monthly_noon_rows[[pollutantVector[poll_index]]]
    #ILO2: scaling columns for fair comparison
    if (all(is.na(pollutant_values))) {
      scaled_matrix[,poll_index] <- NA_real_
    } else {
      min_value <- min(pollutant_values, na.rm=TRUE)
      max_value <- max(pollutant_values, na.rm=TRUE)
      if (max_value > min_value) {
        scaled_matrix[,poll_index] <- (pollutant_values - min_value) / (max_value - min_value)
      } else {
        scaled_matrix[,poll_index] <- 0
      }
    }
  }
  
  # ILO3: I am using Nice smooth color palette so it will be easy to do comparison
  color_scale <- colorRampPalette(c("navy", "cyan", "yellow", "orange", "red"))(100)
  
  # ILO3: I am drawing heatmap to visually compare pollutants
  image(t(apply(scaled_matrix, 2, rev)),
        axes=FALSE,
        col=color_scale ,
        main="Heatmap (monthly noon) for selected pollutants\n0 = low | 1 = high (scaled per pollutant)",
        xlab="Pollutants",
        ylab="Time (old → recent)")
  
  #ILO3: enhancing interpretability of visual output
  axis(1, at=seq(0,1,length.out=total_pol), labels=pollutantVector, las=2, cex.axis=0.8)
  
  reversed_years <- rev(monthly_noon_rows$year)
  year_ticks <- unique(which(!duplicated(reversed_years)))
  if (length(year_ticks) > 1) {
    axis(2,
         at=(year_ticks-1)/(length(reversed_years)-1),
         labels=reversed_years[year_ticks],
         las=2, cex.axis=0.7)
  }
  
  # I am adding color legend for clarity
  par(xpd=TRUE)
  legend("topright",
         legend=c("LOW","HIGH"),
         fill=c(color_scale[1], color_scale[length(color_scale)]),
         border=NA, bty="n", cex=0.9)
  par(xpd=FALSE)
}

# I am calling function to display graph (I can choose any 4 pollutants for comparison)
comparePollutants(c("so2","no2","pm25","o3"))

# Q5 – O3NOCorrelation() (ILO1, ILO2, ILO3, ILO4)
# I am calculating Pearson correlation between O3 and NO2 for the full period.
# I am plotting a SCATTERPLOT (O3 vs NO2) to show the relationship.

O3NOCorrelation <- function() {
  o3_values  <- PortTalbotAQ$o3
  no2_values <- PortTalbotAQ$no2
  
  valid_rows <- complete.cases(o3_values, no2_values) # ILO1: removing invalid rows before analysis
  no2_clean <- no2_values[valid_rows]
  o3_clean <- o3_values[valid_rows]
  
  corr_val <- cor(o3_clean, no2_clean) # ILO4: calculating Pearson correlation measure
  print(paste("Correlation between O3 and NO2 =", round(corr_val,3)))
  
  # I am using Pearson correlation because it is a parametric test.
  # I am applying it to continuous numeric variables (O3 and NO2).
  # I am assuming linear relation + roughly normal distribution, so Pearson is suitable.

  # ILO3: I am using semi-transparent scatter plot visualising relationship
  plot(no2_clean, o3_clean,
       pch = 16,
       col = rgb(231, 76, 60, maxColorValue = 255, alpha = 110),  # tomato with transparency
       xlab = "NO2",
       ylab = "O3",
       main = paste("Scatter of O3 vs NO2 (corr =", round(corr_val, 3), ")"))
  
  # ILO3: I added a simple smoothed trend line (LOWESS) to show pattern clearly
  lines(lowess(no2_clean, o3_clean, f = 2/3), lwd = 2, col = "steelblue")
  # I added light gridlines
  abline(h = pretty(o3_clean), v = pretty(no2_clean), col = "gray90", lty = 3)
  
  return(corr_val)
}

# I am calling function to display graph
O3NOCorrelation()

# Q6 – t-test: O3 by wind speed (<5 vs >=5) (ILO1, ILO2)
# I am testing if mean O3 is different when wind_speed<5 and when wind_speed>=5.
# I am stating H0/H1 and printing the test result and my decision.

# I am writing HYPOTHESES as comments:
# H0 (null): mean O3 is the SAME for low wind (<5 m/s) and high wind (>=5 m/s)
# H1 (alt):  mean O3 is DIFFERENT between those two wind groups

# I am building the two groups now (removing NA)
lowWindGroup  <- PortTalbotAQ$o3[ PortTalbotAQ$wind_speed < 5 ] # ILO1: forming group 1 from dataset
highWindGroup <- PortTalbotAQ$o3[ PortTalbotAQ$wind_speed >= 5 ] #ILO1: forming group 2 from dataset
lowWindGroup  <- lowWindGroup[ !is.na(lowWindGroup) ]
highWindGroup <- highWindGroup[ !is.na(highWindGroup) ]

# ILO4 I am running the 2-sample t-test (two-sided)
wind_ttest_result <- t.test(lowWindGroup, highWindGroup, alternative="two.sided") 
print(wind_ttest_result) # ILO4: interpreting hypothesis decision result

# I am using a t-test, which is a parametric method based on comparing means.
# I am applying it to two numeric samples (O3 values) from two different wind groups.
# I am assuming each sample comes from an approximately normal distribution and variances are not extremely different, so the t-test is suitable here.

# I am writing the decision rule as comments:
# If p-value < 0.05 → I REJECT H0 (there is a significant difference).
# If p-value >= 0.05 → I FAIL TO REJECT H0 (no significant difference found).

# Q7 – ANOVA (NO2 at 17:00 in Nov–Dec for 2016, 2018, 2020) (ILO1, ILO2, ILO3, ILO4)
# I am filtering: hour==17 AND month in {11,12} AND year in {2016,2018,2020}.
# I am running one-way ANOVA: no2 ~ year_group
# If significant, I am doing TukeyHSD to find which pairs differ.
# I am also drawing a BOXPLOT (base R) to show distributions.

# I am writing HYPOTHESES as comments:
# H0 (null): mean NO2 (17:00, Nov–Dec) is the SAME for 2016, 2018, 2020.
# H1 (alt):  at least ONE of the years has a DIFFERENT mean NO2.

# ILO1: filtering specific subset for test
anova_data <- PortTalbotAQ[
  PortTalbotAQ$hour==17 &
    (PortTalbotAQ$month==11 | PortTalbotAQ$month==12) &
    (PortTalbotAQ$year==2016 | PortTalbotAQ$year==2018 | PortTalbotAQ$year==2020),
]

# I am keeping rows with real NO2 values
anova_data <- anova_data[ !is.na(anova_data$no2), ]

# ILO2: I am making the factor (group) variable for comparison
anova_data$year_factor <- as.factor(anova_data$year) 

# I am using ANOVA, which is a parametric test comparing group means.
# I am assuming each group comes from a roughly normal distribution.
# I am assuming equal (or similar) variances between the three groups.
# Because NO2 is numeric and groups are independent by year, ANOVA is appropriate here.

# ILO4: I am performing one-way ANOVA test
anova_result <- aov(no2 ~ year_factor, data=anova_data)
summary(anova_result)

# ILO4: I am doing post-hoc (Tukey) to see WHICH pairs are different (if ANOVA is significant)
tukey_results <- TukeyHSD(anova_result)
tukey_results

# ILO3: I am drawing a BOXPLOT to visualise the three groups for better comarison
boxplot(no2 ~ year_factor,
        data = anova_data,
        main  = "NO2 at 17:00 (Nov–Dec): 2016 vs 2018 vs 2020",
        xlab  = "Year",
        ylab  = "NO2",
        col   = c("#FFCC80", "#90CAF9", "#A5D6A7"),  # soft orange/blue/green
        border = "gray40",
        notch = TRUE)                                # small visual for medians
grid(nx = NA, ny = NULL, col = "gray90", lty = 3)   # horizontal grid only

# I am writing decision comments:
# If ANOVA p-value < 0.05 → I REJECT H0 (there IS a difference).
# Then I read TukeyHSD to say which year-pairs differ (e.g., 2016–2018).
# If ANOVA p-value >= 0.05 → I FAIL TO REJECT H0 (no evidence of difference).