required_packages <- c("ggplot2", "dplyr", "tidyr", "lubridate")
missing_packages <- required_packages[!sapply(required_packages, requireNamespace, quietly = TRUE)]
if (length(missing_packages) > 0) {
  stop(paste("Install:", paste(missing_packages, collapse = ", ")))
}

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

file_path <- "AirQualityUCI.csv"

AirQuality <- read.csv(
  file_path,
  header = TRUE,
  stringsAsFactors = FALSE,
  na.strings = c("NA", "", "-200", "-200.0")
)

AirQuality <- AirQuality[, !grepl("^Unnamed", names(AirQuality)), drop = FALSE]
AirQuality <- AirQuality[, colSums(!is.na(AirQuality)) > 0, drop = FALSE]

names(AirQuality) <- gsub("[^A-Za-z0-9]+", "_", names(AirQuality))
names(AirQuality) <- gsub("^_|_$", "", names(AirQuality))

for (col_name in names(AirQuality)) {
  if (!(col_name %in% c("Date", "Time"))) {
    AirQuality[[col_name]] <- suppressWarnings(as.numeric(AirQuality[[col_name]]))
  }
}

AirQuality$DateTime <- dmy_hms(paste(AirQuality$Date, AirQuality$Time), quiet = TRUE)
AirQuality$Month <- month(AirQuality$DateTime)
AirQuality$Hour <- hour(AirQuality$DateTime)

AirQuality$Season <- case_when(
  AirQuality$Month %in% c(12, 1, 2) ~ "Winter",
  AirQuality$Month %in% c(3, 4, 5) ~ "Spring",
  AirQuality$Month %in% c(6, 7, 8) ~ "Summer",
  AirQuality$Month %in% c(9, 10, 11) ~ "Autumn",
  TRUE ~ NA_character_
)

AirQuality$Season <- factor(AirQuality$Season, levels = c("Winter", "Spring", "Summer", "Autumn"))

co_col <- "CO_GT"
no2_col <- "NO2_GT"
nox_col <- "NOx_GT"
temp_col <- "T"
rh_col <- "RH"
ah_col <- "AH"
benz_col <- "C6H6_GT"
nmhc_col <- "NMHC_GT"
pt08_1 <- "PT08_S1_CO"
pt08_2 <- "PT08_S2_NMHC"
pt08_3 <- "PT08_S3_NOx"
pt08_4 <- "PT08_S4_NO2"
pt08_5 <- "PT08_S5_O3"

plot_list <- list()
plot_names <- c()

add_plot <- function(p, name) {
  plot_list[[length(plot_list) + 1]] <<- p
  plot_names[length(plot_names) + 1] <<- name
}

add_plot(
  ggplot(AirQuality, aes(x = .data[[co_col]])) +
    geom_histogram(bins = 30) +
    labs(title = "Histogram of CO(GT)", x = "CO(GT)", y = "Frequency") +
    theme_minimal(),
  "01_histogram_co"
)

add_plot(
  ggplot(AirQuality, aes(x = .data[[no2_col]])) +
    geom_histogram(bins = 30) +
    labs(title = "Histogram of NO2(GT)", x = "NO2(GT)", y = "Frequency") +
    theme_minimal(),
  "02_histogram_no2"
)

add_plot(
  ggplot(AirQuality, aes(x = .data[[nox_col]])) +
    geom_histogram(bins = 30) +
    labs(title = "Histogram of NOx(GT)", x = "NOx(GT)", y = "Frequency") +
    theme_minimal(),
  "03_histogram_nox"
)

add_plot(
  ggplot(AirQuality, aes(x = .data[[temp_col]])) +
    geom_histogram(bins = 30) +
    labs(title = "Histogram of Temperature", x = "T", y = "Frequency") +
    theme_minimal(),
  "04_histogram_temp"
)

add_plot(
  ggplot(AirQuality, aes(x = .data[[rh_col]])) +
    geom_histogram(bins = 30) +
    labs(title = "Histogram of Relative Humidity", x = "RH", y = "Frequency") +
    theme_minimal(),
  "05_histogram_rh"
)

add_plot(
  ggplot(AirQuality, aes(x = .data[[co_col]])) +
    geom_density() +
    labs(title = "Density Plot of CO(GT)", x = "CO(GT)", y = "Density") +
    theme_minimal(),
  "06_density_co"
)

add_plot(
  ggplot(AirQuality, aes(x = .data[[no2_col]])) +
    geom_density() +
    labs(title = "Density Plot of NO2(GT)", x = "NO2(GT)", y = "Density") +
    theme_minimal(),
  "07_density_no2"
)

add_plot(
  ggplot(AirQuality, aes(x = DateTime, y = .data[[co_col]])) +
    geom_line() +
    labs(title = "Time Series of CO(GT)", x = "DateTime", y = "CO(GT)") +
    theme_minimal(),
  "08_timeseries_co"
)

add_plot(
  ggplot(AirQuality, aes(x = DateTime, y = .data[[no2_col]])) +
    geom_line() +
    labs(title = "Time Series of NO2(GT)", x = "DateTime", y = "NO2(GT)") +
    theme_minimal(),
  "09_timeseries_no2"
)

long_ts <- AirQuality %>%
  select(DateTime, all_of(c(co_col, no2_col))) %>%
  pivot_longer(cols = -DateTime, names_to = "Pollutant", values_to = "Value")

add_plot(
  ggplot(long_ts, aes(x = DateTime, y = Value)) +
    geom_line() +
    facet_wrap(~Pollutant, scales = "free_y", ncol = 1) +
    labs(title = "Time Trends of Two Pollutants", x = "DateTime", y = "Value") +
    theme_minimal(),
  "10_time_trends_two_pollutants"
)

add_plot(
  ggplot(AirQuality, aes(x = .data[[temp_col]], y = .data[[co_col]])) +
    geom_point(alpha = 0.5) +
    geom_smooth(method = "lm", se = FALSE) +
    labs(title = "CO(GT) vs Temperature", x = "T", y = "CO(GT)") +
    theme_minimal(),
  "11_scatter_co_temp"
)

add_plot(
  ggplot(AirQuality, aes(x = .data[[rh_col]], y = .data[[no2_col]])) +
    geom_point(alpha = 0.5) +
    geom_smooth(method = "lm", se = FALSE) +
    labs(title = "NO2(GT) vs RH", x = "RH", y = "NO2(GT)") +
    theme_minimal(),
  "12_scatter_no2_rh"
)

add_plot(
  ggplot(AirQuality, aes(x = .data[[temp_col]], y = .data[[nox_col]])) +
    geom_point(alpha = 0.5) +
    geom_smooth(method = "lm", se = FALSE) +
    labs(title = "NOx(GT) vs Temperature", x = "T", y = "NOx(GT)") +
    theme_minimal(),
  "13_scatter_nox_temp"
)

add_plot(
  ggplot(AirQuality, aes(x = Season, y = .data[[co_col]])) +
    geom_boxplot() +
    labs(title = "Seasonal Boxplot of CO(GT)", x = "Season", y = "CO(GT)") +
    theme_minimal(),
  "14_boxplot_co_season"
)

add_plot(
  ggplot(AirQuality, aes(x = Season, y = .data[[no2_col]])) +
    geom_boxplot() +
    labs(title = "Seasonal Boxplot of NO2(GT)", x = "Season", y = "NO2(GT)") +
    theme_minimal(),
  "15_boxplot_no2_season"
)

hourly_co <- AirQuality %>%
  group_by(Hour) %>%
  summarise(Value = mean(.data[[co_col]], na.rm = TRUE), .groups = "drop")

add_plot(
  ggplot(hourly_co, aes(x = Hour, y = Value)) +
    geom_line() +
    geom_point() +
    labs(title = "Average Hourly CO(GT)", x = "Hour", y = "CO(GT)") +
    theme_minimal(),
  "16_hourly_co"
)

hourly_no2 <- AirQuality %>%
  group_by(Hour) %>%
  summarise(Value = mean(.data[[no2_col]], na.rm = TRUE), .groups = "drop")

add_plot(
  ggplot(hourly_no2, aes(x = Hour, y = Value)) +
    geom_line() +
    geom_point() +
    labs(title = "Average Hourly NO2(GT)", x = "Hour", y = "NO2(GT)") +
    theme_minimal(),
  "17_hourly_no2"
)

corr_vars <- c(co_col, no2_col, nox_col, temp_col, rh_col, ah_col, benz_col)
corr_data <- AirQuality %>% select(all_of(corr_vars)) %>% na.omit()
corr_matrix <- cor(corr_data, method = "pearson")
print(round(corr_matrix, 3))

corr_df <- as.data.frame(as.table(corr_matrix))
names(corr_df) <- c("Var1", "Var2", "Correlation")

add_plot(
  ggplot(corr_df, aes(x = Var1, y = Var2, fill = Correlation)) +
    geom_tile() +
    geom_text(aes(label = round(Correlation, 2))) +
    scale_fill_gradient2(midpoint = 0) +
    labs(title = "Correlation Heatmap", x = "", y = "") +
    theme_minimal(),
  "18_correlation_heatmap"
)

png("19_scatterplot_matrix.png", width = 1400, height = 1200, res = 150)
pairs(AirQuality[, c(co_col, no2_col, nox_col, temp_col, rh_col)], main = "Scatterplot Matrix")
dev.off()

dir.create("plots_output", showWarnings = FALSE)

for (i in seq_along(plot_list)) {
  ggsave(
    filename = file.path("plots_output", paste0(plot_names[i], ".png")),
    plot = plot_list[[i]],
    width = 8,
    height = 5,
    dpi = 300
  )
}

for (i in seq_along(plot_list)) {
  print(plot_list[[i]])
}

cat("Total plots created:", length(plot_list), "\n")
cat("All ggplots saved in folder: plots_output\n")
cat("Scatterplot matrix saved as: 19_scatterplot_matrix.png\n")
print(names(AirQuality))