# Save all console output into a text file
sink("Analysis_Output.txt", split=TRUE)

cat("PORT TALBOT AIR QUALITY ANALYSIS\n")
cat("--------------------------------\n\n")

# Load the packages required for data manipulation, visualisation and statistical analysis.
pkgs <- c("ggplot2","dplyr","lubridate","corrplot")

for(p in pkgs){
  if(!require(p, character.only=TRUE)){
    install.packages(p,repos="https://cloud.r-project.org")
    library(p,character.only=TRUE)
  }
}

# Import the air quality dataset into a data frame while treating "No data" as missing values.
AirQuality <- read.csv(
  "Port_Talbot_Hourly_AQ.csv",
  stringsAsFactors = FALSE,
  na.strings = "No data",
  check.names = TRUE
)

# Display column names to verify that the dataset has been imported correctly.
cat("COLUMN NAMES\n")
print(names(AirQuality))

# Check the data types to identify variables that require conversion.
cat("\nDATA TYPES\n")
str(AirQuality)

# Convert pollutant and weather measurements into numeric format for analysis.
numCols <- c(
"Ozone",
"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"
)

numCols <- intersect(numCols,names(AirQuality))

AirQuality[numCols] <- lapply(AirQuality[numCols],as.numeric)

# Convert the date and time columns into R date-time objects for time-based analysis.
AirQuality$Date <- dmy(AirQuality$Date)
AirQuality$DateTime <- ymd_hms(paste(AirQuality$Date,AirQuality$Time))

# Extract the month from each observation so seasons can be assigned.
AirQuality$Month <- month(AirQuality$Date)

# Create a season variable to compare pollutant concentrations throughout the year.
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",
TRUE ~ "Autumn"
)

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

# Display the first few observations to inspect the imported data.
cat("\nHEAD\n")
print(head(AirQuality))

# Display the number of rows and columns in the dataset.
cat("\nDIMENSIONS\n")
print(dim(AirQuality))

# Produce descriptive statistics to summarise pollutant and weather variables.
cat("\nSUMMARY\n")
print(summary(AirQuality[numCols]))

# Create a folder to store all generated visualisations.
dir.create("Plots",showWarnings=FALSE)

# Calculate monthly average pollutant values to reduce short-term fluctuations and highlight long-term trends.
MonthlyAQ <- AirQuality %>%
mutate(MonthDate=floor_date(Date,"month")) %>%
group_by(MonthDate) %>%
summarise(
Ozone=mean(Ozone,na.rm=TRUE),
NO2=mean(Nitrogen.dioxide,na.rm=TRUE),
WindSpeed=mean(Modelled.Wind.Speed,na.rm=TRUE),
Temperature=mean(Modelled.Temperature,na.rm=TRUE),
.groups="drop"
)

# A line chart is used because it clearly shows changes in ozone concentrations over time.
p1 <- ggplot(MonthlyAQ,aes(MonthDate,Ozone))+
geom_line(colour="steelblue",linewidth=.5)+
theme_minimal()+
labs(title="Monthly Average Ozone",x="Month",y="Ozone (ug/m3)")

ggsave("Plots/Ozone_TimeSeries.png",p1,width=12,height=6,dpi=300)

# A second time series is produced to compare trends in nitrogen dioxide over the same period.
p2 <- ggplot(MonthlyAQ,aes(MonthDate,NO2))+
geom_line(colour="firebrick",linewidth=.5)+
theme_minimal()+
labs(title="Monthly Average Nitrogen Dioxide",x="Month",y="NO2 (ug/m3)")

ggsave("Plots/NO2_TimeSeries.png",p2,width=12,height=6,dpi=300)

# A scatter plot is used to investigate whether ozone concentration is related to wind speed.
# A fitted regression line is included to highlight the overall linear trend.
p3 <- ggplot(AirQuality,
aes(Modelled.Wind.Speed,Ozone))+
geom_point(alpha=.08,size=.3)+
geom_smooth(method="lm",colour="red",se=FALSE)+
theme_minimal()+
labs(title="Ozone vs Wind Speed")

ggsave("Plots/Ozone_WindSpeed.png",p3,width=10,height=6,dpi=300)

# Select pollutant variables for correlation analysis.
corrVars <- c(
"Ozone",
"Nitric.oxide",
"Nitrogen.dioxide",
"PM10.particulate.matter..Hourly.measured.",
"PM2.5.particulate.matter..Hourly.measured."
)

corrVars <- intersect(corrVars,names(AirQuality))

# Remove missing observations before calculating correlations.
corrData <- na.omit(AirQuality[corrVars])

# Calculate the Pearson correlation matrix to measure linear relationships between pollutants.
corrMat <- cor(corrData)

# Display the correlation matrix for interpretation.
cat("\nCORRELATION MATRIX\n")
print(round(corrMat,3))

# A heatmap is used because it provides a clear visual summary of correlation strength and direction.
png("Plots/Correlation_Heatmap.png",1400,1400,res=220)

corrplot(corrMat,
method="color",
type="upper",
addCoef.col="black",
tl.col="black")

dev.off()

# Create a scatter plot to examine the relationship between ozone and nitrogen dioxide.
# The regression line helps identify any overall association.
p4 <- ggplot(corrData,
aes(Ozone,Nitrogen.dioxide))+
geom_point(alpha=.08,size=.3)+
geom_smooth(method="lm",colour="red",se=FALSE)+
theme_minimal()+
labs(title="Ozone vs Nitrogen Dioxide")

ggsave("Plots/Scatter_Correlation.png",p4,width=10,height=6,dpi=300)

# Divide observations into high and low wind speed groups using the median wind speed.
# Using the median creates two groups of similar size for comparison.
medianWind <- median(AirQuality$Modelled.Wind.Speed,na.rm=TRUE)

AirQuality$WindGroup <- ifelse(
AirQuality$Modelled.Wind.Speed>=medianWind,
"High","Low")

cat("\n====================\n")
cat("T TEST\n")
cat("====================\n")
cat("H0: Mean ozone is equal for high and low wind speed.\n")
cat("H1: Mean ozone differs.\n\n")

# Perform a Welch two-sample t-test because ozone concentrations are compared between two independent groups.
tt <- t.test(Ozone~WindGroup,data=AirQuality)

print(tt)

# Reject or fail to reject the null hypothesis based on the p-value.
if(tt$p.value<0.05){
cat("Conclusion: Reject H0.\n")
}else{
cat("Conclusion: Fail to Reject H0.\n")
}

cat("\n====================\n")
cat("ANOVA\n")
cat("====================\n")
cat("H0: Mean ozone is equal across seasons.\n")
cat("H1: At least one season differs.\n\n")

# Perform a one-way ANOVA because ozone concentrations are compared across four independent seasonal groups.
fit <- aov(Ozone~Season,data=AirQuality)

print(summary(fit))

# Conduct Tukey's post hoc test when the ANOVA result is significant to identify which seasons differ.
if(summary(fit)[[1]][1,"Pr(>F)"]<0.05){
cat("\nPost Hoc Test\n")
print(TukeyHSD(fit))
}

# A boxplot is used because it clearly compares the distribution, spread and median ozone concentration across seasons.
p5 <- ggplot(AirQuality,
aes(Season,Ozone,fill=Season))+
geom_boxplot(outlier.alpha=.1)+
theme_minimal()+
theme(legend.position="none")+
labs(title="Seasonal Ozone Distribution")

ggsave("Plots/Ozone_Season_Boxplot.png",p5,width=8,height=6,dpi=300)

# Brief interpretation of the main findings.
cat("\nINTERPRETATION\n")
cat("- Monthly ozone concentrations showed clear seasonal fluctuations but no strong long-term upward or downward trend.\n")
cat("- Monthly nitrogen dioxide concentrations generally declined over the monitoring period with noticeable seasonal peaks.\n")
cat("- Ozone showed a weak positive relationship with wind speed, indicating higher ozone concentrations at higher wind speeds.\n")
cat("- Ozone had a moderate negative correlation with nitrogen dioxide (r = -0.49), while PM10 and PM2.5 showed a strong positive correlation (r = 0.61).\n")
cat("- The t-test found significantly higher ozone concentrations under high wind speed conditions (p < 0.001).\n")
cat("- ANOVA showed significant seasonal differences in ozone concentrations, with Spring generally recording the highest levels (p < 0.001).\n")


sink()