library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)

# initializing pollution data
pollutionData <- read.csv("/home/nlc62/csds133/nlc62/24sp-csds133-nlc62/5-FinalProject/pollution_2000_2023.csv")
pollutionData$State <- factor(pollutionData$State)
pollutionData$County <- factor(pollutionData$County)
#pollutionData$City <- factor(pollutionData$City)
pollutionData$Date <- as.POSIXct(pollutionData$Date, format = "%Y-%m-%d")
pollutionData$Year <- as.numeric(format(pollutionData$Date, "%Y"))
pollutionData <- pollutionData %>% filter(Year >= 2016) %>% filter(Year < 2023)

# initializing weather data
weatherData <- read.csv("/home/nlc62/csds133/nlc62/24sp-csds133-nlc62/5-FinalProject/WeatherEvents_Jan2016-Dec2022.csv")
weatherData$Severity <- factor(weatherData$Severity)
weatherData$Type <- factor(weatherData$Type)
weatherData$AirportCode <- factor(weatherData$AirportCode)
weatherData$TimeZone <- factor(weatherData$TimeZone)
weatherData$City <- factor(weatherData$City)
weatherData$County <- factor(weatherData$County)
weatherData$State <- factor(weatherData$State)
weatherData$ZipCode <- factor(weatherData$ZipCode)

# exploring the data
#levels(weatherData$Severity) # severity levels are "Light", "Moderate", "Heavy", "Severe", "Other", and "UNK" (unknown)
#levels(weatherData$Type) # types of weather are "Cold", "Fog", "Hail", "Precipitation", "Rain", "Snow", and "Storm"
#weatherData %>% filter(Severity == "Other") %>% filter(Type == "Cold" |Type == "Fog" |Type == "Precipitation" | Type =="Rain" | Type =="Snow" | Type =="Storm")
# "Other" severity is only used for hail
#weatherData %>% filter(Severity=="UNK")
#weatherData %>% filter(Severity == "UNK") %>% filter(Type == "Cold" |Type == "Fog" |Type == "Hail" | Type =="Rain" | Type =="Snow" | Type =="Storm")
# "UNK" severity is only used for precipitation

# check what states are included in each data set
#levels(weatherData$State)
#levels(pollutionData$State)
# filtering by state but not necessary since filtering by cities
#weatherData <- weatherData %>% filter(State!="MT") %>% filter(State!="NE") %>% filter(State!="WV")
#pollutionData <- pollutionData %>% filter(State!="Alaska") %>% filter(State!="District of Columbia") %>% filter(State!="Hawaii")

# filtering so that only cities that exist between both data sets are used
pollutionData$StateAbb <- state.abb[match(pollutionData$State, state.name)]
pollutionData$CityAndState <- paste(pollutionData$City, pollutionData$StateAbb)
weatherData$CityAndState <- paste(weatherData$City, weatherData$State)
pollutionCities <- unique(pollutionData$CityAndState)
weatherCities <- unique(weatherData$CityAndState)
overlapCitiesWeatherData <- weatherData %>% filter(weatherData$CityAndState %in% pollutionCities)
overlapCitiesPollutionData <- pollutionData %>% filter(pollutionData$CityAndState %in% weatherCities)
# ordering by city
overlapCitiesPollutionData <- overlapCitiesPollutionData[order(overlapCitiesPollutionData$City),]
overlapCitiesWeatherData <- overlapCitiesWeatherData[order(overlapCitiesWeatherData$City),]

overlapCitiesPollutionData$City <- factor(overlapCitiesPollutionData$City)

# ensure no repeats (don't want to have certain datapoints overrepresented)
overlapCitiesPollutionData <- distinct(overlapCitiesPollutionData)

# checking if cities were filtered correctly
uniqueCities1 <- unique(overlapCitiesPollutionData$City)
uniqueCities2 <- unique(overlapCitiesWeatherData$City)
print(uniqueCities1 %in% uniqueCities2)
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE
# adding columns to separate the start/end dates and times, making data type date
overlapCitiesWeatherData <- overlapCitiesWeatherData %>% separate(StartTime.UTC.,c("startDate", "startTime"), sep = " ", remove = FALSE)
overlapCitiesWeatherData <- overlapCitiesWeatherData %>% separate(EndTime.UTC.,c("endDate", "endTime"), sep = " ", remove = FALSE)
overlapCitiesWeatherData$startDate <- as.POSIXct(overlapCitiesWeatherData$startDate, format="%Y-%m-%d")
overlapCitiesWeatherData$endDate <- as.POSIXct(overlapCitiesWeatherData$endDate, format="%Y-%m-%d")

# calculating amount of time event lasted
overlapCitiesWeatherData$timeEventLasted <- difftime(overlapCitiesWeatherData$EndTime.UTC., overlapCitiesWeatherData$StartTime.UTC., units = "mins")
LosAngelesPollutionData <- overlapCitiesPollutionData %>% filter(City=='Los Angeles') %>% filter(State=="California")
LosAngelesWeatherData <- overlapCitiesWeatherData %>% filter(City == 'Los Angeles') %>% filter(State=="CA")
LosAngelesWeatherData$startDate <- as.POSIXct(LosAngelesWeatherData$startDate, format="%Y-%m-%d")
LosAngelesWeatherData$endDate <- as.POSIXct(LosAngelesWeatherData$endDate, format="%Y-%m-%d")
OriginalLAWeatherData <-overlapCitiesWeatherData %>% filter(City == 'Los Angeles')

# Reset row names
row.names(LosAngelesWeatherData) <- NULL

# loop to separate weather events that span two days, correct time of event
k=nrow(LosAngelesWeatherData)
LosAngelesWeatherData$StartTime.UTC. <- as.POSIXct(LosAngelesWeatherData$StartTime.UTC., format = "%Y-%m-%d %H:%M:%S")


# HERE I need help figuring out how to make this work. I'm trying to make it such that if a weather event spans between two days (ex: 23:40 first day to 03:27 the second day, I want to split these up into their own weather events so that we can look at each day specifically, and change the time of the event to reflect this.) I think it also would be good for events with precipitation to try to split the amount of precipitation by the same proportion as the time is split between the two days.
for(i in 1:k) {
  if(as.Date(LosAngelesWeatherData$EndTime.UTC.[i]) > as.Date(LosAngelesWeatherData$StartTime.UTC.[i])) {
    LosAngelesWeatherData <- rbind(LosAngelesWeatherData, LosAngelesWeatherData[i,])
    row.names(LosAngelesWeatherData) <- NULL
    
    date_only <- as.Date(LosAngelesWeatherData$StartTime.UTC.[i])
    next_day <- date_only + 1
    midnight_next_day <- as.POSIXct(paste(next_day, "00:00:00"), format = "%Y-%m-%d %H:%M:%S")
    time_diff <- difftime(midnight_next_day, LosAngelesWeatherData$StartTime.UTC.[i], units = "mins")
    
    LosAngelesWeatherData$timeEventLasted[nrow(LosAngelesWeatherData)] <- LosAngelesWeatherData$timeEventLasted[i]-time_diff
    LosAngelesWeatherData$timeEventLasted[i] <- time_diff
    
    LosAngelesWeatherData$StartTime.UTC.[nrow(LosAngelesWeatherData)] <- as.POSIXct(paste(LosAngelesWeatherData$endDate[i], "00:00:01"))#, format = "%Y-%m-%d %H:%M:%S")
    LosAngelesWeatherData$startDate[nrow(LosAngelesWeatherData)] <- as.POSIXct(LosAngelesWeatherData$endDate[i], format = "%Y-%m-%d")
    LosAngelesWeatherData$startTime[nrow(LosAngelesWeatherData)] <- "00:00:00"
    LosAngelesWeatherData$EndTime.UTC.[i] <- (paste((LosAngelesWeatherData$startDate[i]),"23:59:59"))#, format = "%Y-%m-%d %H:%M:%S")
    LosAngelesWeatherData$endDate[i] <- as.POSIXct(LosAngelesWeatherData$startDate[i], format = "%Y-%m-%d")
    LosAngelesWeatherData$endTime[i] <- "23:59:59"
  } 
}

# planning how to calculate times
# Extract the date portion from the datetime column
#date_only <- as.Date(LosAngelesWeatherData$StartTime.UTC.[27])
# Add one day to the date portion to get the date of the next day
#next_day <- date_only + 1
# Combine the next day's date with midnight time to get midnight of the next day
#midnight_next_day <- as.POSIXct(paste(next_day, "00:00:00"), format = "%Y-%m-%d %H:%M:%S")
# Calculate the time difference from midnight of the next day for each row
#time_diff <- difftime(midnight_next_day, LosAngelesWeatherData$StartTime.UTC.[27], units = "mins")
#print(time_diff)
#print(midnight_next_day)
#print(LosAngelesWeatherData$StartTime.UTC.[27])
#print(difftime(next_day,LosAngelesWeatherData$StartTime.UTC.[27], units = "mins"))
#print(difftime(as.POSIXct(paste(as.Date(LosAngelesWeatherData$StartTime.UTC.[27]+1), "00:00:00")),LosAngelesWeatherData$StartTime.UTC.[27], units = "mins"))
#print(difftime(as.POSIXct(LosAngelesWeatherData$endTime[5855], format = "%H:%M:%D"), as.POSIXct(LosAngelesWeatherData$startTime[5855],format = "%H:%M:%D"), units = "hours"))
#print(difftime(LosAngelesWeatherData$endTime[5855], LosAngelesWeatherData$startTime[5855], units = "minutes"))
#LosAngelesWeatherData$timeEventLasted <- difftime(LosAngelesWeatherData$EndTime.UTC., LosAngelesWeatherData$StartTime.UTC., units = "mins")

# graph for Los Angeles

# loop to put true/false if weather event occurred on a given day in LA

LosAngelesWeatherData <- LosAngelesWeatherData[order(LosAngelesWeatherData$startDate),]
row.names(LosAngelesWeatherData) <- NULL

LosAngelesPollutionData$weatherEvents <- NA
LosAngelesPollutionData <- LosAngelesPollutionData[!duplicated(LosAngelesPollutionData$Date), ]
c
## function (...)  .Primitive("c")
j <- 1
i <- 1
for(i in 1:nrow(LosAngelesPollutionData)) {
  while((as.numeric(difftime(LosAngelesPollutionData$Date[i],LosAngelesWeatherData$startDate[j], units="days"))) >0 && j <= nrow(LosAngelesWeatherData)) {
    j <- j+1
  }
  if(as.Date(LosAngelesPollutionData$Date[i])==as.Date(LosAngelesWeatherData$startDate[j]) && j <= nrow(LosAngelesWeatherData)) {
    LosAngelesPollutionData$weatherEvents[i] <- TRUE
  } 
  else {
    LosAngelesPollutionData$weatherEvents[i] <- FALSE
  }
}

# making separate data sets for days with vs without weather events
LosAngelesDataNoWeatherEvents <- LosAngelesPollutionData %>% filter(weatherEvents==FALSE)
LosAngelesDataNoWeatherEvents <- LosAngelesDataNoWeatherEvents[!duplicated(LosAngelesDataNoWeatherEvents$Date), ]

LosAngelesDataWeatherEvents <- LosAngelesPollutionData %>% filter(weatherEvents==TRUE)
LosAngelesDataWeatherEvents <- LosAngelesDataWeatherEvents[!duplicated(LosAngelesDataWeatherEvents$Date), ]
#LosAngelesDataWeatherEvents$Date <- as.Date(LosAngelesDataWeatherEvents$Date, format="%Y-%m-%d")
#LosAngelesWeatherData$startDate <- as.Date(LosAngelesWeatherData$startDate, format="%Y-%m-%d")
# violin plot and histogram: no weather events, SO2 AQI
p <- ggplot(LosAngelesDataNoWeatherEvents, aes(x=City, y=SO2.AQI)) + 
  geom_violin()
p

hist(LosAngelesDataNoWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in LA on Days with No Recorded Weather")

# violin plot and histogram: no weather events, O3 AQI
p <- ggplot(LosAngelesDataNoWeatherEvents, aes(x=City, y=O3.AQI)) + 
  geom_violin()
p

hist(LosAngelesDataNoWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in LA on Days with No Recorded Weather")

# violin plot and histogram: recorded weather events, SO2 AQI
p <- ggplot(LosAngelesDataWeatherEvents, aes(x=City, y=SO2.AQI)) + 
  geom_violin()
p

hist(LosAngelesDataWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in LA on Days with Recorded Weather")

# violin plot and histogram: recorded weather events, O3 AQI
p <- ggplot(LosAngelesDataWeatherEvents, aes(x=City, y=O3.AQI)) + 
  geom_violin()
p

hist(LosAngelesDataWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in LA on Days with Recorded Weather")

# cleaning data for the loop (order, duplicates for pollution)
row.names(LosAngelesDataWeatherEvents) <- NULL
LosAngelesDataWeatherEvents <- LosAngelesDataWeatherEvents[!duplicated(LosAngelesDataWeatherEvents$Date), ]
LosAngelesDataWeatherEvents <- LosAngelesDataWeatherEvents[order(LosAngelesDataWeatherEvents$Date),]
LosAngelesWeatherData <- LosAngelesWeatherData[order(LosAngelesWeatherData$startDate),]

# calculate total time of weather events for the day
LosAngelesDataWeatherEvents$timeLasted <- NA
j <- 1
i <- 1
for(i in 1:nrow(LosAngelesDataWeatherEvents)) {
  p <- 0
  while((as.numeric(difftime(as.Date(LosAngelesDataWeatherEvents$Date[i]),as.Date(LosAngelesWeatherData$startDate[j]), units="days")))>0 && j <= nrow(LosAngelesWeatherData)) {
    j <- j+1
  }
  while((as.numeric(difftime(as.Date(LosAngelesDataWeatherEvents$Date[i]),as.Date(LosAngelesWeatherData$startDate[j]), units="days")))==0 && j <= nrow(LosAngelesWeatherData)) {
    p <- p + LosAngelesWeatherData$timeEventLasted[j]
    j <- j+1
  }
  LosAngelesDataWeatherEvents$timeLasted[i] <- p
  p <- 0
}
# plotting boxplots for O3 AQI for all cities
boxplot(O3.AQI~City,
        data=overlapCitiesPollutionData,
        main="O3 AQI for Different Cities",
        xlab="City",
        ylab="O3 AQI",
        col="pink",
        border="purple", las = 2
)

# boxplots for SO2 AQI for all cities shared between weather and pollution datasets
boxplot(SO2.AQI~City,
        data=overlapCitiesPollutionData,
        main="SO2 AQI for Different Cities",
        xlab="City",
        ylab="SO2 AQI",
        col="pink",
        border="purple", las = 2
)

# boxplots for NO2 AQI for all cities shared between weather and pollution datasets
boxplot(NO2.AQI~City,
        data=overlapCitiesPollutionData,
        main="NO2 AQI for Different Cities",
        xlab="City",
        ylab="NO2 AQI",
        col="pink",
        border="purple", las = 2
)

# boxplots for CO AQI for all cities shared between weather and pollution datasets
boxplot(CO.AQI~City,
        data=overlapCitiesPollutionData,
        main="CO AQI for Different Cities",
        xlab="City",
        ylab="CO AQI",
        col="pink",
        border="purple", las = 2
)

ggplot(overlapCitiesPollutionData, aes(x = O3.AQI, y = SO2.AQI, color=City)) + 
  geom_point() +
  labs(x = "O3 AQI", y = "SO2 AQI",
       title = "Relationship between SO2 AQI and O3 AQI",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0))

ggplot(overlapCitiesPollutionData, aes(x = O3.AQI, y = CO.AQI, color=City)) + 
  geom_point() +
  labs(x = "O3 AQI", y = "CO AQI",
       title = "Relationship between CO AQI and O3 AQI",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0))

ggplot(overlapCitiesPollutionData, aes(x = NO2.AQI, y = CO.AQI, color=City)) + 
  geom_point() +
  labs(x = "NO2 AQI", y = "CO AQI",
       title = "Relationship between NO2 AQI and O3 AQI",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0))

# Plot heatmap
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
pollutionMeasures <- overlapCitiesPollutionData[, c("O3.AQI", "SO2.AQI", "NO2.AQI", "CO.AQI")]
correlation_matrix <- cor(pollutionMeasures)
correlation_long <- melt(correlation_matrix)
ggplot(data = correlation_long, aes(x = Var1, y = Var2, fill = value)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(-1,1), space = "Lab",
                       name="Pearson\nCorrelation") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 1, 
                                   size = 12, hjust = 1)) +
  labs(x = "Variables", y = "Variables")

# graph time of weather events vs SO2 AQI
summary(LosAngelesDataWeatherEvents)
##        X               Date                          Address         
##  Min.   :438394   Min.   :2016-01-01 00:00:00.00   Length:1283       
##  1st Qu.:470028   1st Qu.:2017-11-15 12:00:00.00   Class :character  
##  Median :530951   Median :2019-08-02 00:00:00.00   Mode  :character  
##  Mean   :533139   Mean   :2019-07-20 18:44:31.24                     
##  3rd Qu.:593778   3rd Qu.:2021-04-24 12:00:00.00                     
##  Max.   :626077   Max.   :2022-12-31 00:00:00.00                     
##                                                                      
##         State                  County              City         O3.Mean       
##  California:1283   Los Angeles    :1283   Los Angeles:1283   Min.   :0.00075  
##  Alabama   :   0   Ada            :   0   Albuquerque:   0   1st Qu.:0.01961  
##  Alaska    :   0   Adair          :   0   Baton Rouge:   0   Median :0.02829  
##  Arizona   :   0   Adams          :   0   Birmingham :   0   Mean   :0.02722  
##  Arkansas  :   0   Alameda        :   0   Bismarck   :   0   3rd Qu.:0.03471  
##  Colorado  :   0   Alexandria City:   0   Cheyenne   :   0   Max.   :0.06577  
##  (Other)   :   0   (Other)        :   0   (Other)    :   0                    
##  O3.1st.Max.Value  O3.1st.Max.Hour     O3.AQI          CO.Mean      
##  Min.   :0.00200   Min.   : 7.00   Min.   :  2.00   Min.   :0.1458  
##  1st Qu.:0.03300   1st Qu.:10.00   1st Qu.: 31.00   1st Qu.:0.2896  
##  Median :0.04200   Median :10.00   Median : 39.00   Median :0.3917  
##  Mean   :0.04198   Mean   :10.48   Mean   : 41.74   Mean   :0.4596  
##  3rd Qu.:0.05100   3rd Qu.:11.00   3rd Qu.: 47.00   3rd Qu.:0.5675  
##  Max.   :0.11800   Max.   :23.00   Max.   :214.00   Max.   :1.3958  
##                                                                     
##  CO.1st.Max.Value CO.1st.Max.Hour      CO.AQI          SO2.Mean      
##  Min.   :0.2000   Min.   : 0.000   Min.   : 2.000   Min.   :-0.4667  
##  1st Qu.:0.4000   1st Qu.: 1.000   1st Qu.: 5.000   1st Qu.: 0.1167  
##  Median :0.5000   Median : 6.000   Median : 6.000   Median : 0.2625  
##  Mean   :0.6266   Mean   : 6.712   Mean   : 7.185   Mean   : 0.3028  
##  3rd Qu.:0.8000   3rd Qu.: 9.000   3rd Qu.: 9.000   3rd Qu.: 0.4521  
##  Max.   :1.8000   Max.   :23.000   Max.   :20.000   Max.   : 1.4917  
##                                                                      
##  SO2.1st.Max.Value SO2.1st.Max.Hour    SO2.AQI           NO2.Mean     
##  Min.   :-0.3000   Min.   : 0.00    Min.   : 0.0000   Min.   : 2.729  
##  1st Qu.: 0.3000   1st Qu.: 7.00    1st Qu.: 0.0000   1st Qu.:11.825  
##  Median : 0.6000   Median :11.00    Median : 0.0000   Median :16.279  
##  Mean   : 0.7158   Mean   :10.37    Mean   : 0.3297   Mean   :18.486  
##  3rd Qu.: 0.9000   3rd Qu.:13.00    3rd Qu.: 0.0000   3rd Qu.:24.036  
##  Max.   :17.9000   Max.   :23.00    Max.   :24.0000   Max.   :55.604  
##                                                                       
##  NO2.1st.Max.Value NO2.1st.Max.Hour    NO2.AQI           Year     
##  Min.   : 5.10     Min.   : 0.00    Min.   : 5.00   Min.   :2016  
##  1st Qu.:22.95     1st Qu.: 6.00    1st Qu.:21.50   1st Qu.:2017  
##  Median :30.10     Median : 9.00    Median :28.00   Median :2019  
##  Mean   :31.86     Mean   :11.22    Mean   :29.67   Mean   :2019  
##  3rd Qu.:39.80     3rd Qu.:20.00    3rd Qu.:37.00   3rd Qu.:2021  
##  Max.   :80.00     Max.   :23.00    Max.   :79.00   Max.   :2022  
##                                                                   
##    StateAbb         CityAndState       weatherEvents    timeLasted    
##  Length:1283        Length:1283        Mode:logical   Min.   :   3.0  
##  Class :character   Class :character   TRUE:1283      1st Qu.: 113.5  
##  Mode  :character   Mode  :character                  Median : 264.0  
##                                                       Mean   : 432.1  
##                                                       3rd Qu.: 625.5  
##                                                       Max.   :2981.0  
## 
plot(LosAngelesDataWeatherEvents$timeLasted, LosAngelesDataWeatherEvents$SO2.AQI)

LosAngelesDataNoWeatherEvents$timeLasted <- 0


# scatter plot of SO2 AQI values vs time of weather events
mergedData <- rbind(LosAngelesDataNoWeatherEvents, LosAngelesDataWeatherEvents)
ggplot(mergedData, aes(x = timeLasted, y = SO2.AQI, color=weatherEvents)) + 
    geom_point() +
    scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
    labs(x = "Time of Weather Events in a Day", y = "SO2 AQI",
         title = "Daily SO2 AQI in Los Angeles on Days with vs without Weather Events",
         caption = "Data: Kaggle") +
    theme(plot.title = element_text(size = 14, face = "bold.italic"),
          plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
              values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(mergedData, aes(x = Date, y = SO2.AQI, color=weatherEvents)) + 
  geom_point() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(x = "Date", y = "SO2 AQI",
       title = "SO2 AQI from 2016-2022 in Los Angeles",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
                     values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# boxplot of SO2 AQI for LA 
boxplot(SO2.AQI~weatherEvents,
        data=LosAngelesPollutionData,
        main="SO2 AQI for LA on Days With vs Without Weather Events",
        xlab="Weather Events Occured",
        ylab="SO2 AQI",
        col="purple",
        border="blue", las = 2
)

# boxplot of O3 AQI for LA
boxplot(O3.AQI~weatherEvents,
        data=LosAngelesPollutionData,
        main="O3 AQI for LA on Days With vs Without Weather Events",
        xlab="Weather Events Occured",
        ylab="O3 AQI",
        col="purple",
        border="blue", las = 2
)

# adding fog as variable to los angeles pollution dataset
LosAngelesFogData <- LosAngelesWeatherData %>% filter(Type=="Fog")
fogDays <- unique(LosAngelesFogData$startDate)
LAFogDaysPollution <- LosAngelesPollutionData %>% filter(Date %in% fogDays)
LosAngelesPollutionData$Fog <- ifelse(LosAngelesPollutionData$Date %in% fogDays, "Foggy", "No fog")

# graphing pollution levels on foggy vs no fog days
ggplot(LosAngelesPollutionData, aes(x = Fog, y = SO2.AQI, color=Fog)) + 
  geom_point() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(x = "Date", y = "SO2 AQI",
       title = "SO2 AQI on Foggy and Non-Foggy Days in Los Angeles",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(breaks = c("Foggy", "Not foggy"),
                     values=c("cyan2", "darkslategray"))

ggplot(LosAngelesPollutionData, aes(x = Fog, y = O3.AQI, color=Fog)) + 
  geom_point() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(x = "Date", y = "O3 AQI",
       title = "O3 AQI on Foggy and Non-Foggy Days in Los Angeles",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(breaks = c("Foggy", "Not foggy"),
                     values=c("cyan2", "darkslategray"))

library(ggstatsplot)
## You can cite this package as:
##      Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
##      Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.5
## ✔ lubridate 1.9.3     ✔ stringr   1.5.1
## ✔ purrr     1.0.2     ✔ tibble    3.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# violin plot with statistical analysis of fog vs no fog for LA
plt <- ggbetweenstats(
  data = LosAngelesPollutionData,
  x = Fog,
  y = SO2.AQI
) + labs(
  x = "Presence of Fog",
  y = "SO2 AQI",
  title = "SO2 AQI Compared Between Foggy and Non-Foggy in Los Angeles"
)
plt

plt2 <- ggbetweenstats(
  data = LosAngelesPollutionData,
  x = Fog,
  y = O3.AQI
) + labs(
  x = "Presence of Fog",
  y = "O3 AQI",
  title = "O3 AQI Compared Between Foggy and Non-Foggy in Los Angeles"
)
plt2

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(LosAngelesPollutionData, aes(x = Date, y = SO2.AQI, color=Fog)) + 
  geom_point() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(x = "Date", y = "SO2 AQI",
       title = "SO2 AQI from 2016-2022 in Los Angeles",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(name = "Legend", breaks = c("Foggy", "No fog"),
                     values=c("cyan2", "darkslategray"), labels = c("Foggy", "No fog"))

# scatter plot of O3 AQI values over time with distinction between days with vs without weather events
ggplot(LosAngelesPollutionData, aes(x = Date, y = O3.AQI, color=Fog)) + 
  geom_point() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(x = "Date", y = "SO2 AQI",
       title = "O3 AQI from 2016-2022 in Los Angeles",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(name = "Legend", breaks = c("Foggy", "No fog"),
                     values=c("cyan2", "darkslategray"), labels = c("Foggy", "No fog"))

DenverPollutionData <- overlapCitiesPollutionData %>% filter(City=='Denver')%>% filter(State=='Colorado')
DenverWeatherData <- overlapCitiesWeatherData %>% filter(City == 'Denver') %>% filter(State=='CO')
DenverWeatherData$startDate <- as.POSIXct(DenverWeatherData$startDate, format="%Y-%m-%d")
DenverWeatherData$endDate <- as.POSIXct(DenverWeatherData$endDate, format="%Y-%m-%d")

# Reset row names
row.names(DenverWeatherData) <- NULL

# loop to separate weather events that span two days, correct time of event
k=nrow(DenverWeatherData)
DenverWeatherData$StartTime.UTC. <- as.POSIXct(DenverWeatherData$StartTime.UTC., format = "%Y-%m-%d %H:%M:%S")

for(i in 1:k) {
  if(as.Date(DenverWeatherData$EndTime.UTC.[i]) > as.Date(DenverWeatherData$StartTime.UTC.[i])) {
    DenverWeatherData <- rbind(DenverWeatherData, DenverWeatherData[i,])
    row.names(DenverWeatherData) <- NULL
    
    date_only <- as.Date(DenverWeatherData$StartTime.UTC.[i])
    next_day <- date_only + 1
    midnight_next_day <- as.POSIXct(paste(next_day, "00:00:00"), format = "%Y-%m-%d %H:%M:%S")
    time_diff <- difftime(midnight_next_day, DenverWeatherData$StartTime.UTC.[i], units = "mins")
    
    DenverWeatherData$timeEventLasted[nrow(DenverWeatherData)] <- DenverWeatherData$timeEventLasted[i]-time_diff
    DenverWeatherData$timeEventLasted[i] <- time_diff
    
    DenverWeatherData$StartTime.UTC.[nrow(DenverWeatherData)] <- as.POSIXct(paste(DenverWeatherData$endDate[i], "00:00:01"))#, format = "%Y-%m-%d %H:%M:%S")
    DenverWeatherData$startDate[nrow(DenverWeatherData)] <- as.POSIXct(DenverWeatherData$endDate[i], format = "%Y-%m-%d")
    DenverWeatherData$startTime[nrow(DenverWeatherData)] <- "00:00:00"
    DenverWeatherData$EndTime.UTC.[i] <- (paste((DenverWeatherData$startDate[i]),"23:59:59"))#, format = "%Y-%m-%d %H:%M:%S")
    DenverWeatherData$endDate[i] <- as.POSIXct(DenverWeatherData$startDate[i], format = "%Y-%m-%d")
    DenverWeatherData$endTime[i] <- "23:59:59"
  } 
}


# loop to put true/false if weather event occurred on a given day in LA

DenverWeatherData <- DenverWeatherData[order(DenverWeatherData$startDate),]
row.names(DenverWeatherData) <- NULL

DenverPollutionData$weatherEvents <- NA
DenverPollutionData <- DenverPollutionData[!duplicated(DenverPollutionData$Date), ]


j <- 1
i <- 1
for(i in 1:nrow(DenverPollutionData)) {
  while((as.numeric(difftime(DenverPollutionData$Date[i],DenverWeatherData$startDate[j], units="days"))) >0 && j <= nrow(DenverWeatherData)) {
    j <- j+1
  }
  if(as.Date(DenverPollutionData$Date[i])==as.Date(DenverWeatherData$startDate[j]) && j <= nrow(DenverWeatherData)) {
    DenverPollutionData$weatherEvents[i] <- TRUE
  } 
  else {
    DenverPollutionData$weatherEvents[i] <- FALSE
  }
}

# making separate data sets for days with vs without weather events
DenverDataNoWeatherEvents <- DenverPollutionData %>% filter(weatherEvents==FALSE)
DenverDataNoWeatherEvents <- DenverDataNoWeatherEvents[!duplicated(DenverDataNoWeatherEvents$Date), ]

DenverDataWeatherEvents <- DenverPollutionData %>% filter(weatherEvents==TRUE)
DenverDataWeatherEvents <- DenverDataWeatherEvents[!duplicated(DenverDataWeatherEvents$Date), ]
# cleaning data for the loop (order, duplicates for pollution)
row.names(DenverDataWeatherEvents) <- NULL
DenverDataWeatherEvents <- DenverDataWeatherEvents[!duplicated(DenverDataWeatherEvents$Date), ]
DenverDataWeatherEvents <- DenverDataWeatherEvents[order(DenverDataWeatherEvents$Date),]
DenverWeatherData <- DenverWeatherData[order(DenverWeatherData$startDate),]

# calculate total time of weather events for the day
DenverDataWeatherEvents$timeLasted <- NA
j <- 1
i <- 1
for(i in 1:nrow(DenverDataWeatherEvents)) {
  p <- 0
  while((as.numeric(difftime(as.Date(DenverDataWeatherEvents$Date[i]),as.Date(DenverWeatherData$startDate[j]), units="days")))>0 && j <= nrow(DenverWeatherData)) {
    j <- j+1
  }
  while((as.numeric(difftime(as.Date(DenverDataWeatherEvents$Date[i]),as.Date(DenverWeatherData$startDate[j]), units="days")))==0 && j <= nrow(DenverWeatherData)) {
    p <- p + DenverWeatherData$timeEventLasted[j]
    j <- j+1
  }
  DenverDataWeatherEvents$timeLasted[i] <- p
  p <- 0
}
# violin plot and histogram: no weather events, SO2 AQI
p <- ggplot(DenverDataNoWeatherEvents, aes(x=City, y=SO2.AQI)) + 
  geom_violin()
p

hist(DenverDataNoWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in Denver on Days with No Recorded Weather")

# violin plot and histogram: no weather events, O3 AQI
p <- ggplot(DenverDataNoWeatherEvents, aes(x=City, y=O3.AQI)) + 
  geom_violin()
p

hist(DenverDataNoWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in Denver on Days with No Recorded Weather")

# violin plot and histogram: recorded weather events, SO2 AQI
p <- ggplot(DenverDataWeatherEvents, aes(x=City, y=SO2.AQI)) + 
  geom_violin()
p

hist(DenverDataWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in Denver on Days with Recorded Weather")

# violin plot and histogram: recorded weather events, O3 AQI
p <- ggplot(DenverDataWeatherEvents, aes(x=City, y=O3.AQI)) + 
  geom_violin()
p

hist(DenverDataWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in Denver on Days with Recorded Weather")

ClevelandPollutionData <- overlapCitiesPollutionData %>% filter(City=='Cleveland') %>% filter(State=='Ohio')
ClevelandWeatherData <- overlapCitiesWeatherData %>% filter(City == 'Cleveland') %>% filter(State=="OH")
ClevelandWeatherData$startDate <- as.POSIXct(ClevelandWeatherData$startDate, format="%Y-%m-%d")
ClevelandWeatherData$endDate <- as.POSIXct(ClevelandWeatherData$endDate, format="%Y-%m-%d")

# Reset row names
row.names(ClevelandWeatherData) <- NULL

# loop to separate weather events that span two days, correct time of event
k=nrow(ClevelandWeatherData)
ClevelandWeatherData$StartTime.UTC. <- as.POSIXct(ClevelandWeatherData$StartTime.UTC., format = "%Y-%m-%d %H:%M:%S")

for(i in 1:k) {
  if(as.Date(ClevelandWeatherData$EndTime.UTC.[i]) > as.Date(ClevelandWeatherData$StartTime.UTC.[i])) {
    ClevelandWeatherData <- rbind(ClevelandWeatherData, ClevelandWeatherData[i,])
    row.names(ClevelandWeatherData) <- NULL
    
    date_only <- as.Date(ClevelandWeatherData$StartTime.UTC.[i])
    next_day <- date_only + 1
    midnight_next_day <- as.POSIXct(paste(next_day, "00:00:00"), format = "%Y-%m-%d %H:%M:%S")
    time_diff <- difftime(midnight_next_day, ClevelandWeatherData$StartTime.UTC.[i], units = "mins")
    
    ClevelandWeatherData$timeEventLasted[nrow(ClevelandWeatherData)] <- ClevelandWeatherData$timeEventLasted[i]-time_diff
    ClevelandWeatherData$timeEventLasted[i] <- time_diff
    
    ClevelandWeatherData$StartTime.UTC.[nrow(ClevelandWeatherData)] <- as.POSIXct(paste(ClevelandWeatherData$endDate[i], "00:00:00"))#, format = "%Y-%m-%d %H:%M:%S") # need to be 00 or 01??
    ClevelandWeatherData$startDate[nrow(ClevelandWeatherData)] <- as.POSIXct(ClevelandWeatherData$endDate[i], format = "%Y-%m-%d")
    ClevelandWeatherData$startTime[nrow(ClevelandWeatherData)] <- "00:00:00"
    ClevelandWeatherData$EndTime.UTC.[i] <- (paste((ClevelandWeatherData$startDate[i]),"23:59:59"))#, format = "%Y-%m-%d %H:%M:%S")
    ClevelandWeatherData$endDate[i] <- as.POSIXct(ClevelandWeatherData$startDate[i], format = "%Y-%m-%d")
    ClevelandWeatherData$endTime[i] <- "23:59:59"
  } 
}

# loop to put true/false if weather event occurred on a given day in Cleveland

ClevelandWeatherData <- ClevelandWeatherData[order(ClevelandWeatherData$startDate),]
row.names(ClevelandWeatherData) <- NULL

ClevelandPollutionData$weatherEvents <- NA
ClevelandPollutionData <- ClevelandPollutionData[!duplicated(ClevelandPollutionData$Date), ]
c
## function (...)  .Primitive("c")
j <- 1
i <- 1
for(i in 1:nrow(ClevelandPollutionData)) {
  while((as.numeric(difftime(ClevelandPollutionData$Date[i],ClevelandWeatherData$startDate[j], units="days"))) >0 && j <= nrow(ClevelandWeatherData)) {
    j <- j+1
  }
  if(as.Date(ClevelandPollutionData$Date[i])==as.Date(ClevelandWeatherData$startDate[j]) && j <= nrow(ClevelandWeatherData)) {
    ClevelandPollutionData$weatherEvents[i] <- TRUE
  } 
  else {
    ClevelandPollutionData$weatherEvents[i] <- FALSE
  }
}

# making separate data sets for days with vs without weather events
ClevelandDataNoWeatherEvents <- ClevelandPollutionData %>% filter(weatherEvents==FALSE)
ClevelandDataNoWeatherEvents <- ClevelandDataNoWeatherEvents[!duplicated(ClevelandDataNoWeatherEvents$Date), ]

ClevelandDataWeatherEvents <- ClevelandPollutionData %>% filter(weatherEvents==TRUE)
ClevelandDataWeatherEvents <- ClevelandDataWeatherEvents[!duplicated(ClevelandDataWeatherEvents$Date), ]
# cleaning data for the loop (order, duplicates for pollution)
row.names(ClevelandDataWeatherEvents) <- NULL
ClevelandDataWeatherEvents <- ClevelandDataWeatherEvents[!duplicated(ClevelandDataWeatherEvents$Date), ]
ClevelandDataWeatherEvents <- ClevelandDataWeatherEvents[order(ClevelandDataWeatherEvents$Date),]
ClevelandWeatherData <- ClevelandWeatherData[order(ClevelandWeatherData$startDate),]

# calculate total time of weather events for the day
ClevelandDataWeatherEvents$timeLasted <- NA
j <- 1
i <- 1
for(i in 1:nrow(ClevelandDataWeatherEvents)) {
  p <- 0
  while((as.numeric(difftime(as.Date(ClevelandDataWeatherEvents$Date[i]),as.Date(ClevelandWeatherData$startDate[j]), units="days")))>0 && j <= nrow(ClevelandWeatherData)) {
    j <- j+1
  }
  while((as.numeric(difftime(as.Date(ClevelandDataWeatherEvents$Date[i]),as.Date(ClevelandWeatherData$startDate[j]), units="days")))==0 && j <= nrow(ClevelandWeatherData)) {
    p <- p + ClevelandWeatherData$timeEventLasted[j]
    j <- j+1
  }
  ClevelandDataWeatherEvents$timeLasted[i] <- p
  p <- 0
}
# violin plot and histogram: no weather events, SO2 AQI
p <- ggplot(ClevelandDataNoWeatherEvents, aes(x=City, y=SO2.AQI)) + 
  geom_violin()
p

hist(ClevelandDataNoWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in Cleveland on Days with No Recorded Weather")

# violin plot and histogram: no weather events, O3 AQI
p <- ggplot(ClevelandDataNoWeatherEvents, aes(x=City, y=O3.AQI)) + 
  geom_violin()
p

hist(ClevelandDataNoWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in Cleveland on Days with No Recorded Weather")

# violin plot and histogram: recorded weather events, SO2 AQI
p <- ggplot(ClevelandDataWeatherEvents, aes(x=City, y=SO2.AQI)) + 
  geom_violin()
p

hist(ClevelandDataWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in Cleveland on Days with Recorded Weather")

# violin plot and histogram: recorded weather events, O3 AQI
p <- ggplot(ClevelandDataWeatherEvents, aes(x=City, y=O3.AQI)) + 
  geom_violin()
p

hist(ClevelandDataWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in Cleveland on Days with Recorded Weather")

# graph time of weather events vs SO2 AQI
plot(ClevelandDataWeatherEvents$timeLasted, ClevelandDataWeatherEvents$SO2.AQI)

ClevelandDataNoWeatherEvents$timeLasted <- 0

# scatter plot of SO2 AQI values vs time of weather events
mergedClevelandData <- rbind(ClevelandDataNoWeatherEvents, ClevelandDataWeatherEvents)
ggplot(mergedClevelandData, aes(x = timeLasted, y = SO2.AQI, color=weatherEvents)) + 
    geom_point() +
    scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
    labs(x = "Time of Weather Events in a Day", y = "SO2 AQI",
         title = "Daily SO2 AQI in Cleveland on Days with vs without Weather Events",
         caption = "Data: Kaggle") +
    theme(plot.title = element_text(size = 14, face = "bold.italic"),
          plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
              values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(mergedClevelandData, aes(x = Date, y = SO2.AQI, color=weatherEvents)) + 
  geom_point() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(x = "Date", y = "SO2 AQI",
       title = "SO2 AQI from 2016-2022 in Cleveland",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
                     values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# boxplot of SO2 AQI for Cleveland 
boxplot(SO2.AQI~weatherEvents,
        data=ClevelandPollutionData,
        main="SO2 AQI for Cleveland on Days With vs Without Weather Events",
        xlab="Weather Events Occured",
        ylab="SO2 AQI",
        col="purple",
        border="blue", las = 2
)

# boxplot of O3 AQI for Cleveland
boxplot(O3.AQI~weatherEvents,
        data=ClevelandPollutionData,
        main="O3 AQI for Cleveland on Days With vs Without Weather Events",
        xlab="Weather Events Occured",
        ylab="O3 AQI",
        col="purple",
        border="blue", las = 2
)

# adding fog as variable to Cleveland pollution dataset
ClevelandFogData <- ClevelandWeatherData %>% filter(Type=="Fog")
fogDays <- unique(ClevelandFogData$startDate)
ClevelandFogDaysPollution <- ClevelandPollutionData %>% filter(Date %in% fogDays)
ClevelandPollutionData$Fog <- ifelse(ClevelandPollutionData$Date %in% fogDays, "Foggy", "No fog")

# graphing pollution levels on foggy vs no fog days
ggplot(ClevelandPollutionData, aes(x = Fog, y = SO2.AQI, color=Fog)) + 
  geom_point() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(x = "Date", y = "SO2 AQI",
       title = "SO2 AQI on Foggy and Non-Foggy Days in Cleveland",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(breaks = c("Foggy", "No fog"),
                     values=c("cyan2", "darkslategray"))

library(ggstatsplot)
library(tidyverse)

# violin plot with statistical analysis of fog vs no fog for Cleveland
plt <- ggbetweenstats(
  data = ClevelandPollutionData,
  x = Fog,
  y = SO2.AQI
) + labs(
  x = "Presence of Fog",
  y = "SO2 AQI",
  title = "SO2 AQI Compared Between Foggy and Non-Foggy in Cleveland"
)
plt

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(ClevelandPollutionData, aes(x = Date, y = SO2.AQI, color=Fog)) + 
  geom_point() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(x = "Date", y = "SO2 AQI",
       title = "SO2 AQI from 2016-2022 in Cleveland",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(name = "Legend", breaks = c("Foggy", "No fog"),
                     values=c("cyan2", "darkslategray"), labels = c("Foggy", "No fog"))

summary(ClevelandDataWeatherEvents$timeLasted)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     7.0   180.0   580.5   864.5  1287.0  4971.0
# graph time of weather events vs SO2 AQI
plot(DenverDataWeatherEvents$timeLasted, DenverDataWeatherEvents$SO2.AQI)

DenverDataNoWeatherEvents$timeLasted <- 0


# scatter plot of SO2 AQI values vs time of weather events
mergedDenverData <- rbind(DenverDataNoWeatherEvents, DenverDataWeatherEvents)
ggplot(mergedDenverData, aes(x = timeLasted, y = SO2.AQI, color=weatherEvents)) + 
    geom_point() +
    scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
    labs(x = "Time of Weather Events in a Day", y = "SO2 AQI",
         title = "Daily SO2 AQI in Denver on Days with vs without Weather Events",
         caption = "Data: Kaggle") +
    theme(plot.title = element_text(size = 14, face = "bold.italic"),
          plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
              values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(mergedDenverData, aes(x = Date, y = SO2.AQI, color=weatherEvents)) + 
  geom_point() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(x = "Date", y = "SO2 AQI",
       title = "SO2 AQI from 2016-2022 in Denver",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
                     values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# boxplot of SO2 AQI for Denver 
boxplot(SO2.AQI~weatherEvents,
        data=DenverPollutionData,
        main="SO2 AQI for Denver on Days With vs Without Weather Events",
        xlab="Weather Events Occured",
        ylab="SO2 AQI",
        col="purple",
        border="blue", las = 2
)

# boxplot of O3 AQI for Denver
boxplot(O3.AQI~weatherEvents,
        data=DenverPollutionData,
        main="O3 AQI for Denver on Days With vs Without Weather Events",
        xlab="Weather Events Occured",
        ylab="O3 AQI",
        col="purple",
        border="blue", las = 2
)

# adding fog as variable to Denver pollution dataset
DenverFogData <- DenverWeatherData %>% filter(Type=="Fog")
fogDays <- unique(DenverFogData$startDate)
DenverFogDaysPollution <- DenverPollutionData %>% filter(Date %in% fogDays)
DenverPollutionData$Fog <- ifelse(DenverPollutionData$Date %in% fogDays, "Foggy", "No fog")

# graphing pollution levels on foggy vs no fog days
ggplot(DenverPollutionData, aes(x = Fog, y = SO2.AQI, color=Fog)) + 
  geom_point() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(x = "Date", y = "SO2 AQI",
       title = "SO2 AQI on Foggy and Non-Foggy Days in Denver",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(breaks = c("Foggy", "No fog"),
                     values=c("cyan2", "darkslategray"))

library(ggstatsplot)
library(tidyverse)

# violin plot with statistical analysis of fog vs no fog for LA
plt <- ggbetweenstats(
  data = DenverPollutionData,
  x = Fog,
  y = SO2.AQI
) + labs(
  x = "Presence of Fog",
  y = "SO2 AQI",
  title = "SO2 AQI Compared Between Foggy and Non-Foggy in Denver"
)
plt

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(DenverPollutionData, aes(x = Date, y = SO2.AQI, color=Fog)) + 
  geom_point() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(x = "Date", y = "SO2 AQI",
       title = "SO2 AQI from 2016-2022 in Denver",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(name = "Legend", breaks = c("Foggy", "No fog"),
                     values=c("cyan2", "darkslategray"), labels = c("Foggy", "No fog"))

pltDen <- ggbetweenstats(
  data = DenverPollutionData,
  x = weatherEvents,
  y = SO2.AQI
) + labs(
  x = "Presence of Weather Events",
  y = "SO2 AQI",
  title = "SO2 AQI Between Days With and Without Weather in Denver"
)
plt

pltLA <- ggbetweenstats(
  data = LosAngelesPollutionData,
  x = weatherEvents,
  y = SO2.AQI
) + labs(
  x = "Presence of Weather Events",
  y = "SO2 AQI",
  title = "SO2 AQI Between Days With and Without Weather in LA"
)
pltLA

pltCle <- ggbetweenstats(
  data = ClevelandPollutionData,
  x = weatherEvents,
  y = SO2.AQI
) + labs(
  x = "Presence of Weather Events",
  y = "SO2 AQI",
  title = "SO2 AQI Between Days With and Without Weather in Cleveland"
)
pltCle

CleLADenPollution <- rbind(ClevelandPollutionData, LosAngelesPollutionData, DenverPollutionData)

GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, 
                           draw_group = function(self, data, ..., draw_quantiles = NULL) {
  data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
  grp <- data[1, "group"]
  newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
  newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
  newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])

  if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
    stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
      1))
    quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
    aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
    aesthetics$alpha <- rep(1, nrow(quantiles))
    both <- cbind(quantiles, aesthetics)
    quantile_grob <- GeomPath$draw_panel(both, ...)
    ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
  }
  else {
    ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
  }
})

geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., 
                              draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, 
                              show.legend = NA, inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, 
        position = position, show.legend = show.legend, inherit.aes = inherit.aes, 
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}

ggplot(CleLADenPollution, aes(City, SO2.AQI, fill = weatherEvents)) + geom_split_violin()

ggplot(data = CleLADenPollution, aes (x = weatherEvents, y = SO2.AQI, fill=weatherEvents)) + geom_violin(position = "dodge") + facet_wrap(~City)

t.test(DenverDataNoWeatherEvents$SO2.AQI, DenverDataWeatherEvents$SO2.AQI, alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  DenverDataNoWeatherEvents$SO2.AQI and DenverDataWeatherEvents$SO2.AQI
## t = 6.0243, df = 1864, p-value = 2.042e-09
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.5138864 1.0099953
## sample estimates:
## mean of x mean of y 
##  2.119004  1.357063
t.test(LosAngelesDataNoWeatherEvents$SO2.AQI, LosAngelesDataWeatherEvents$SO2.AQI, alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  LosAngelesDataNoWeatherEvents$SO2.AQI and LosAngelesDataWeatherEvents$SO2.AQI
## t = 1.3113, df = 2390.9, p-value = 0.1899
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.02453863  0.12360448
## sample estimates:
## mean of x mean of y 
##  0.379229  0.329696
t.test(ClevelandDataNoWeatherEvents$SO2.AQI, ClevelandDataWeatherEvents$SO2.AQI, alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  ClevelandDataNoWeatherEvents$SO2.AQI and ClevelandDataWeatherEvents$SO2.AQI
## t = 8.1233, df = 1074.5, p-value = 1.237e-15
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  2.731518 4.471355
## sample estimates:
## mean of x mean of y 
##  6.885204  3.283768
DenverFogPollution <- DenverPollutionData %>% filter(Fog == "Foggy")
DenverNoFogPollution <- DenverPollutionData %>% filter(Fog == "No fog")
t.test(DenverFogPollution$SO2.AQI, DenverNoFogPollution$SO2.AQI, alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  DenverFogPollution$SO2.AQI and DenverNoFogPollution$SO2.AQI
## t = -2.7745, df = 560.04, p-value = 0.005713
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.6451957 -0.1103281
## sample estimates:
## mean of x mean of y 
##  1.520124  1.897886
LANoFogPollution <- LosAngelesPollutionData %>% filter(Fog == "No fog")
t.test(LAFogDaysPollution$SO2.AQI, LANoFogPollution$SO2.AQI, alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  LAFogDaysPollution$SO2.AQI and LANoFogPollution$SO2.AQI
## t = 1.4167, df = 1568.5, p-value = 0.1568
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.02289013  0.14194291
## sample estimates:
## mean of x mean of y 
## 0.3900293 0.3305029
ClevelandNoFog <- ClevelandPollutionData %>% filter(Fog=="No fog")
t.test(ClevelandFogDaysPollution$SO2.AQI, ClevelandNoFog$SO2.AQI, alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  ClevelandFogDaysPollution$SO2.AQI and ClevelandNoFog$SO2.AQI
## t = -3.4819, df = 1399.5, p-value = 0.0005131
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.0662502 -0.5770584
## sample estimates:
## mean of x mean of y 
##  3.478571  4.800226
plot1 <- ggplot(data = LosAngelesPollutionData, aes(x=Fog, y=SO2.AQI)) + geom_boxplot(aes(fill=Fog)) +   labs(x = "", y = "SO2 AQI",
       title = "Los Angeles") + theme(legend.position="none") + ylim(0,70)
plot2 <- ggplot(data = DenverPollutionData, aes(x=Fog, y=SO2.AQI)) + geom_boxplot(aes(fill=Fog)) +   labs(x = "", y = "SO2 AQI",
       title = "Denver") + theme(legend.position="none")
plot3 <- ggplot(data = ClevelandPollutionData, aes(x=Fog, y=SO2.AQI)) + geom_boxplot(aes(fill=Fog)) +   labs(x = "", y = "SO2 AQI",
       title = "Cleveland") + theme(legend.position="none")
require(gridExtra)
## Loading required package: gridExtra
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
grid.arrange(plot1, plot2, plot3, ncol=3)

LAFogDaysPollution$Fog <- "Foggy"
DenverFogDaysPollution$Fog <- "Foggy"
ClevelandFogDaysPollution$Fog <- "Foggy"
mergedFogData <- rbind(LAFogDaysPollution, DenverFogDaysPollution, ClevelandFogDaysPollution)
mergedNoFogData <- rbind(LANoFogPollution, DenverNoFogPollution, ClevelandNoFog)

mergedFogAndNoFogData <- rbind(mergedFogData, mergedNoFogData)
t.test(mergedFogData$SO2.AQI, mergedNoFogData$SO2.AQI, alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  mergedFogData$SO2.AQI and mergedNoFogData$SO2.AQI
## t = -5.5553, df = 4081.8, p-value = 2.948e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.040657 -0.497734
## sample estimates:
## mean of x mean of y 
##  1.625122  2.394318
ggplot(mergedFogAndNoFogData, aes(x = Fog, y = SO2.AQI, color=Fog)) + 
  geom_violin() +
  scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
  labs(y = "SO2 AQI",
       title = "SO2 AQI on Foggy and Non-Foggy Days in Los Angeles and Denver",
       caption = "Data: Kaggle") +
  theme(plot.title = element_text(size = 14, face = "bold.italic"),
        plot.caption = element_text(hjust = 0)) + 
  scale_color_manual(breaks = c("Foggy", "No fog"),
                     values=c("cyan2", "darkslategray"))

plt <- ggbetweenstats(
  data = mergedFogAndNoFogData,
  x = Fog,
  y = SO2.AQI
) + labs(
  x = "Presence of Fog",
  y = "SO2 AQI",
  title = "SO2 AQI Compared Between Foggy and Non-Foggy in Los Angeles, Denver, and Cleveland"
)
plt

ggplot(pollutionData, aes(x = Date, y = SO2.AQI, color=City)) + 
  geom_point() + theme(legend.position="none")

ggplot(overlapCitiesPollutionData, aes(x = Date, y = SO2.AQI)) +
  geom_line(position = "dodge", stat = "identity", aes(col = City)) +
  facet_wrap(~City) +
  theme(legend.position = "none") +
  labs(title = "SO2 AQI vs Date", x = "Date", y = "SO2 AQI", caption = "source: Kaggle")
## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`

ggplot(pollutionData, aes(x = Date, y = O3.AQI, color=City)) + 
  geom_point() + theme(legend.position="none")

ggplot(overlapCitiesPollutionData, aes(x = Date, y = O3.AQI)) +
  geom_line(position = "dodge", stat = "identity", aes(col = City)) +
  facet_wrap(~City) +
  theme(legend.position = "none") +
  labs(title = "O3 AQI vs Date", x = "Date", y = "O3 AQI", caption = "source: Kaggle")
## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`

ggplot(pollutionData, aes(x = Date, y = CO.AQI, color=City)) + 
  geom_point() + theme(legend.position="none")

ggplot(overlapCitiesPollutionData, aes(x = Date, y = CO.AQI)) +
  geom_line(position = "dodge", stat = "identity", aes(col = City)) +
  facet_wrap(~City) +
  theme(legend.position = "none") +
  labs(title = "CO AQI vs Date", x = "Date", y = "CO AQI", caption = "source: Kaggle")
## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`

ggplot(pollutionData, aes(x = Date, y =  NO2.AQI, color=City)) + 
  geom_point() + theme(legend.position="none")

ggplot(overlapCitiesPollutionData, aes(x = Date, y = NO2.AQI)) +
  geom_line(position = "dodge", stat = "identity", aes(col = City)) +
  facet_wrap(~City) +
  theme(legend.position = "none") +
  labs(title = "NO2 AQI vs Date", x = "Date", y = "NO2 AQI", caption = "source: Kaggle")
## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`