############# Air pollution ############# ## Student: Adam Smith ################## ## Student No: 2317631 ################## ########################################### ############################################### ############# REQUIRED PACKAGES ############# ##THIS SCRIPT WILL REQUIRE THE INSTALLATION OF THE PACKAGE READXL ##THIS SCRIPT WILL REQUIRE THE INSTALLATION OF TIDYVERSE PACKAGE WHICH INCLUDES GGPLOT AND DPLYR library(readxl) ## TO READ IN XLSX FILES library(tidyverse) ## FOR CLEANING AND MANIPULATING DATA AND PRODUCING GRAPHS VIA GGPLOT ################################################## ############# Function definitions ############# ## Is_invalid - function to identify the number of "no data" entries in data provided. ## DisplayPollutantGraph - function which can be called on a specific pollutant to generate a graph between ## July 2015 and November 2024. ## comparePollutants - function which will compare pollutants - requires a vector with the pollutant names to be ## created separately to reference, and when called will need the start year. ## O3NOCorrelation - function which will create a scatterplot of the correlation between Nitric Oxide and Ozone ## Will need dataframe when called. ####################################### ############# FUNCTIONS ############# ## Function to count invalids - function counts 'No data' as na and creates a vector and stores the count is_invalid <- function(x) { is.na(x) | x == 'No data' } ## Function to create graph for pollutant displayPollutantGraph <- function(pollutant, start_year, end_year, start_month, end_month) { ## Create full start and end dates using inputs start_date <- paste0(start_year, start_month, '01', sep = '-') end_date <- paste0(end_year, end_month, '01', sep = '-') ## generate function - pull pollutant, start point & end point filtered_data <- PortTalbotAQ |> filter(Date >= start_date & Date <= end_date) |> filter(day(Date) == 1 & Time == '12:00:00') |> select(Date, all_of(pollutant)) ## Filter data to be between the years specified and at noon for the 1st of the month ## Plot graph using this data ggplot(filtered_data, aes(x = as.Date(Date), y = as.numeric(.data[[pollutant]]), group = 1)) + geom_line(colour = 'blue', linewidth = 1) + ## line colour geom_point(colour = 'darkblue') + ## point colour labs(title = paste0("Variation in ", pollutant, " levels"), ## Title subtitle = paste0("Measured at noon on the 1st of each month"), ## subtitle x = "Date", ## x axis title y = paste0(pollutant, " Concentration") + ## y axis title scale_x_date(date_breaks = "1 year", date_labels = "%Y") + ## reduce data labels on x axis theme( plot.title = element_text(hjust = 0.5), ## centre plot.subtitle = element_text(hjust = 0.5), ## centre axis.text.x = element_text(angle = 45, hjust = 1) ## angle text ) ) } ## Function to create graph comparePollutants comparePollutants <- function(pollutants, start_year) { ## Define the timescale - last 10 years end_year <- start_year + 10 ## Filter data to be correct time period - the 1st of each month filtered_data <- PortTalbotAQ |> filter(Date >= start_year & Date <= end_year) |> filter(day(Date) == 1 & Time == "12:00:00") |> select(Date, all_of(pollutants)) |> ## Pivot data into longer format for comparison - this is to allow ggplot as wont populate otherwise pivot_longer( cols = all_of(pollutants), names_to = 'Pollutant', values_to = 'Concentration' ) ## plot the graph from the data ggplot( filtered_data, aes( x = as.Date(Date), y = as.numeric(Concentration), colour = Pollutant ) ) + geom_line(linewidth = 1) + geom_point() + labs( title = paste0( "Comparison of Pollutant Levels (", ## Title start_year, "-", end_year, ")", sep = "" ), subtitle = paste0("Measured at noon on the 1st of each month"), ## Subtitle x = "Date", ## x axis y = "Concentration" ## y axis ) + scale_x_date(date_breaks = "1 year", date_labels = "%Y") + ## reduce data labels on x axis theme( plot.title = element_text(hjust = 0.5), ## centre plot.subtitle = element_text(hjust = 0.5), ## centre axis.text.x = element_text(angle = 45, hjust = 1), ## angle text legend_position = "bottom" ## legend below plot line ) + facet_wrap(vars(Pollutant)) ## split graphs into respective groups } ## Function to create graph O3NOCorrelation O3NOCorrelation <- function(x) { ## Ensure the columns are numeric to avoid errors filtered_data <- x ozone <- as.numeric(filtered_data$O3) nitrogen_dioxide <- as.numeric(filtered_data$NO2) ## Calculate and print correlation coefficient cor_coeff <- cor(ozone, nitrogen_dioxide, use = "complete.obs") print(cor_coeff) ## Create scatterplot ggplot(PortTalbotAQ, aes(x = nitrogen_dioxide, y = ozone)) + geom_point(alpha = 0.4, colour = 'blue') + ## set colour geom_smooth(method = 'lm', colour = 'red', linetype = 'dashed') + ## set colour labs(title = "Correlation between Ozone (O3) and Nitrogen Dioxide (NO2)", ## title subtitle = paste0("Pearson Correlation Coefficient:", round(cor_coeff, 4)), ## subtitle x = "Nitrogen Dioxide (NO2) Concentration", ## X axis title y = "Ozone (O3) Concentration") + ## y axis title theme_minimal() ## apply minimal theme - keeps graph tidy and looking professional } ######################################### ############# MAIN SCRIPT ############# ## set working directory for where files are stored on my / this device #setwd("C:/Users/23628/OneDrive/Desktop/R directory") #amend direction of slashes as throws an error if copied from location ## Load data file into dataframe #PortTalbotHourlyAQ<-read.csv("PortTalbotHourlyAQ.csv", header=TRUE, stringsAsFactors = TRUE) ## this brings through the whole data file allowing you to read the additional data, but doesn't ## bring across the headers, and throws the columns out. ##### Task 1 ##### ## Load just data file without the first 10 lines. PortTalbotAQ <- read_xlsx ("PortTalbotHourlyAQ.xlsx", skip = 10) ## changed to use readxl as read.csv was importing additional columns of NA's ## Skips the first 10 rows of the data file - setting the columns correctly - Doesn't allow you to read the first ## lines so both have been created. ## strings as factors is not available in readxl ## header = true is not available in readxl ## Warnings will appear where "status" columns will have to change and be numbered. These can be ignored. ##rename columns to shortened forms. Names are pretty succinct, so have used chemical formulae. names(PortTalbotAQ) <- c("Date", "Time", "O3", "Status", "NO", "Status.1", "NO2", "Status.2", "NO_as_NO2", "status.3", "SO2", "Status.4", "CO", "Status.5", "PM10_Prtclt_mtr_hrly", "Status.6", "PM2.5_Prtclt_mtr_hrly", "Status.7", "Wind_drctn", "Status.8", "Wind_Spd", "Status.9", "Temp", "Status.10") ## In the main file the column "status" are not numbered but R requires them be numbered to differentiate them. ##### Task 2 ##### ## calculate the percentage of invalid data entries - this takes the average and then multiplies by 100 to provide ## the percentage and then enters the data into a dataframe called invalidData. invalidData <- data.frame( Column = names(PortTalbotAQ), Percent_Invalid = sapply(PortTalbotAQ, function(x) mean(is_invalid(x)) * 100) ) ## Creates dataframe invalidData and pulls across the headers from the DF(PortTalbotAQ) - calculates a column ## called Percent_Invalid which uses the function to pull invalid data (No Data) takes the mean value and x 100 ## to calculate the average. ## This brings across all headers from the dataframe - Assignment specifically asks for just pollutants. invalidData <- invalidData[c(3,5,7,9,11,13), ] ##Rebuild the invalidData dataframe with only the pollutants as requested. ## Display the dataframe. print(invalidData) ##### Task 3 ##### ## call function to produce graph, specify the pollutant needed and dates required (July 2015 - Nov 2024) displayPollutantGraph('NO_as_NO2', 2014, 2024, 6, 11) ##### Task 4 ##### ## create a vector of the 4 pollutants to compare polls <- c('NO', 'NO2', 'SO2', 'CO') ## call function to compare with the vector of pollutants and start year (will graph 9 years from this point) comparePollutants(polls, 2014) ##### Task 5 ##### ## call function to produce graph comparing O3 and NO2 O3NOCorrelation(PortTalbotAQ) ##### Task 6 ##### ## T-test to look at significance between wind speed and ozone levels. ## standard T Test hypothesis is that the means of both groups will be equal, or not have a statistical significance between them - my experimental ## hypothesis would be the reverse - that there IS a statistical difference between the two groups. ## Create a classification column for wind and remove any NA's in the dataset o3_test <- PortTalbotAQ |> mutate(wind_speed_group = as.factor(ifelse(`Wind_Spd` < 5, "low speed", "high speed"))) ## create o3 test with additional column - split Wind+Spd into 2 groups <5 or >5 name them "low speed" and "high speed" ## Run T-test using wind_speed_group and O3 columns result <- t.test(as.numeric(O3) ~ wind_speed_group, data = o3_test, na.rm = TRUE) result ## the P value is equal to .000000000000000022 which is significantly below the threshold of 0.05 - as such the null hypothesis is rejected - there IS a ## statistically significant difference between groups. ## Task 7 ##### ## ANOVA test looking at nitrogen dioxide levels over three distinct periods of time ## Create vector of just the required time and dates to allow for easier filtering ## Filter data to correct years, months in the year and readings taken at 17:00 only anova_test <- PortTalbotAQ |> mutate(Date = as.Date(Date), Year = year(Date), Month = month(Date)) |> filter(hour(hms(Time)) == 17, Month %in% c(11, 12), Year %in% c(2016, 2018, 2020)) |> ## change Year column into a factor for ANOVA test mutate(Year = as.factor(Year)) ## Perform ANOVA test result_anova <- aov(as.numeric(NO2) ~ Year, data = anova_test) ## Print ANOVA test results summary(result_anova) ## Significance value of 7.2e-07 is equal to 0.000000072 - as such the null hypothesis is rejected - there IS a significant difference between the years. ## ANOVA test confirms a difference exists but it doesn't identify where - to identify where the difference lies additional tests (TukeyHSD) are needed. ## Perform TukeyHSD test on the results from the Anova test. TukeyHSD(result_anova) ## The 'p adj' column displays significance values. ## There is a significant difference between 2016-2018 (0.0000202) and 2016-2020 (0.0000057) ## There is no significant difference between 2018-2020 (0.9553023)