bmore_arduino_aqi.Rmd

---
title: "bmore_arduino_aqi"
author: "Leona Neftaliem"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(message = FALSE, warning = FALSE, error = FALSE)
```

#Libraries
```{r}
library(dplyr)
library(lubridate)
library(tidyverse)
library(dplyr)
library(zoo)
library(data.table)
library(ggmap)
library(ggplot2)
library(viridis)
```


#Clean
```{r}
# Set the directory where your files are located
data_dir <- "~/Dropbox/Bmore_AQI/Data/"

# Get a list of all files in the directory that start with "AQLSP"
file_list <- list.files(path = data_dir, pattern = "^AQLSP", full.names = TRUE)

# Read all files and combine them into one data table
aqi_data <- rbindlist(lapply(file_list, fread))

# Rename columns to make them more readable
setnames(aqi_data, c("timestamp", "arduino_id", "program", "bme_temp_c", "bme_pressure_hpa", 
                     "bme_humidity_percent", "bme_voc_ppm", "bme_altitude_m", "k30_co2_ppm", 
                     "ozone_ppb", "sensor1_pm1_ugm3", "sensor1_pm2.5_ugm3", "sensor1_pm10_ugm3", 
                     "sensor2_pm1_ugm3", "sensor2_pm2.5_ugm3", "sensor2_pm10_ugm3"))

# Step 2: Add New Columns with Default Values
aqi_data[, `:=`(
  long = NA_real_,
  lat = NA_real_,
  neighborhood = NA_character_,
  vulnerability_score = NA_real_,
  vulnerability_level = NA_character_,
  land_use_code = NA_character_
)]

# Reorder the columns for better organization
setcolorder(aqi_data, c("timestamp", "long", "lat", "neighborhood", 
                        "vulnerability_score", "vulnerability_level", 
                        "land_use_code", "arduino_id", "program", 
                        "bme_temp_c", "bme_pressure_hpa", "bme_humidity_percent", 
                        "bme_voc_ppm", "bme_altitude_m", "k30_co2_ppm", "ozone_ppb", 
                        "sensor1_pm1_ugm3", "sensor1_pm2.5_ugm3", "sensor1_pm10_ugm3",
                        "sensor2_pm1_ugm3", "sensor2_pm2.5_ugm3", "sensor2_pm10_ugm3"))

# Step 3: Data Cleaning - Replace invalid values (-7999) with NA
aqi_data[aqi_data == -7999] <- NA


# Step 4: Convert Timestamps to Standard Format
convert_timestamp <- function(ts) {
  ts_parsed <- ymd_hms(ts, tz = "UTC")
  if (is.na(ts_parsed)) {
    ts_parsed <- mdy_hm(ts, tz = "UTC")
  }
  return(ts_parsed)
}

# Ensure timestamp is character before conversion
aqi_data[, timestamp := as.character(timestamp)]

# Apply the conversion function
aqi_data[, timestamp := sapply(timestamp, convert_timestamp)]

aqi_data <- aqi_data[!is.na(timestamp)]

# Ensure timestamp is in POSIXct format
aqi_data[, timestamp := as.POSIXct(timestamp)]

aqi_data <- aqi_data[!duplicated(aqi_data$timestamp), ]

# Floor to the nearest hour
aqi_data[, timestamp_hour := floor_date(timestamp, unit = "hour")]

# Check the cleaned data
head(aqi_data)
aqi_data$ozone_ppb[is.nan(aqi_data$ozone_ppb)]<-NA

#Ozone comes in ppb-- converting now to ppm
aqi_data$ozone_ppm <- aqi_data$ozone_ppb /1000

# Step 6: Compute Hourly Averages for PM2.5 and PM10 Data
hourly_avg <- aqi_data[, .(
  sensor1_pm2.5_ugm3_avg = mean(sensor1_pm2.5_ugm3, na.rm = TRUE),
  sensor2_pm2.5_ugm3_avg = mean(sensor2_pm2.5_ugm3, na.rm = TRUE),
  sensor1_pm10_ugm3_avg  = mean(sensor1_pm10_ugm3, na.rm = TRUE),
  sensor2_pm10_ugm3_avg  = mean(sensor2_pm10_ugm3, na.rm = TRUE),
  bme_temp_c_avg  = mean(bme_temp_c, na.rm = TRUE),
  bme_humidity_percent_avg  = mean(bme_humidity_percent, na.rm = TRUE),
  ozone_ppm_avg  = mean(ozone_ppm, na.rm = TRUE)
), by = timestamp_hour]

# Convert the timestamp to daily timestamps
aqi_data[, timestamp_day := floor_date(timestamp, unit = "day")]


# Compute Daily Averages for PM2.5, PM10, Temperature, Humidity, and Ozone Data
daily_avg <- aqi_data[, .(
  sensor1_pm2.5_ugm3_avg = mean(sensor1_pm2.5_ugm3, na.rm = TRUE),
  sensor2_pm2.5_ugm3_avg = mean(sensor2_pm2.5_ugm3, na.rm = TRUE),
  sensor1_pm10_ugm3_avg  = mean(sensor1_pm10_ugm3, na.rm = TRUE),
  sensor2_pm10_ugm3_avg  = mean(sensor2_pm10_ugm3, na.rm = TRUE),
  bme_temp_c_avg  = mean(bme_temp_c, na.rm = TRUE),
  bme_humidity_percent_avg  = mean(bme_humidity_percent, na.rm = TRUE),
  ozone_ppm_avg  = mean(ozone_ppm, na.rm = TRUE)
), by = timestamp_day]

# View the first few rows of the daily average data
head(daily_avg)


```

#PM Hourly Graphics
```{r}
#PM2.5
# Convert hourly_avg to data.frame for ggplot
hourly_avg <- as.data.frame(hourly_avg)
hourly_avg <- hourly_avg %>%
  mutate(day = floor_date(timestamp_hour, unit = "day"))



# Step 6: Plot the Hourly Averaged Data for PM2.5 with vibrant colors
pm25 <- ggplot(hourly_avg, aes(x = timestamp_hour)) +
  geom_line(aes(y = sensor1_pm2.5_ugm3_avg, color = "Sensor 1 PM2.5"), size = 1) +  # First y axis
  geom_line(aes(y = sensor2_pm2.5_ugm3_avg, color = "Sensor 2 PM2.5"), size = 1, linetype = "dashed") +  # Second y axis
  scale_y_continuous(
    name = "Sensor 1 PM2.5 (µg/m³)", 
    sec.axis = sec_axis(~ ., name = "Sensor 2 PM2.5 (µg/m³)")
  ) +
  scale_color_manual(values = c("Sensor 1 PM2.5" = "#FF5733", "Sensor 2 PM2.5" = "#007AFF")) + 
  labs(x = "Timestamp (Hourly Average)", 
       color = "PM2.5 Sensors", 
       title = "Hourly Averaged PM2.5 Concentrations",
       subtitle = "Comparison of Sensor 1 and Sensor 2") +
  theme_minimal() +
  theme(
    axis.title.y = element_text(color = "#FF5733"),
    axis.title.y.right = element_text(color = "#007AFF"),
    legend.position = "top",
    plot.title = element_text(size = 18, face = "bold"),
    plot.subtitle = element_text(size = 14),
    axis.text = element_text(size = 12)
  )

# Step 7: Display the Plot
print(pm25)

#PM10
# Plot the Hourly Averaged Data for PM10 with vibrant colors
pm10 <- ggplot(hourly_avg, aes(x = timestamp_hour)) +
  geom_line(aes(y = sensor1_pm10_ugm3_avg, color = "Sensor 1 PM10"), size = 1) +  # First y axis
  geom_line(aes(y = sensor2_pm10_ugm3_avg, color = "Sensor 2 PM10"), size = 1, linetype = "dashed") +  # Second y axis
  scale_y_continuous(
    name = "Sensor 1 PM10 (µg/m³)", 
    sec.axis = sec_axis(~ ., name = "Sensor 2 PM10 (µg/m³)")
  ) +
  scale_color_manual(values = c("Sensor 1 PM10" = "#FF5733", "Sensor 2 PM10" = "#007AFF")) +
  labs(x = "Timestamp (Hourly Average)", 
       color = "PM10 Sensors", 
       title = "Hourly Averaged PM10 Concentrations",
       subtitle = "Comparison of Sensor 1 and Sensor 2") +
  theme_minimal() +
  theme(
    axis.title.y = element_text(color = "#FF5733"),
    axis.title.y.right = element_text(color = "#007AFF"),
    legend.position = "top",
    plot.title = element_text(size = 18, face = "bold"),
    plot.subtitle = element_text(size = 14),
    axis.text = element_text(size = 12)
  )

# Display the PM10 plot
print(pm10)

#PM2.5 and Temp (C)

p_temp_pm25_heatmap <- ggplot(hourly_avg, aes(x = bme_temp_c_avg, y = sensor1_pm2.5_ugm3_avg)) +
  geom_bin2d(bins = 30, aes(fill = ..count..)) +
  scale_fill_viridis_c(option = "plasma") +
  labs(x = "Average Temperature (°C)", 
       y = "Average PM2.5 (µg/m³)",
       fill = "Count",
       title = "Heatmap of PM2.5 Concentrations") +
  theme_minimal(base_size = 14)

p_temp_pm25_heatmap

#-------
# Load necessary library for correlation calculation
library(ggpubr)

# Calculate correlation between PM2.5 and temperature
cor_temp_pm25_sensor1 <- cor(hourly_avg$bme_temp_c_avg, hourly_avg$sensor1_pm2.5_ugm3_avg, use = "complete.obs")
cor_temp_pm25_sensor2 <- cor(hourly_avg$bme_temp_c_avg, hourly_avg$sensor2_pm2.5_ugm3_avg, use = "complete.obs")

# Calculate correlation between PM2.5 and humidity
cor_humidity_pm25_sensor1 <- cor(hourly_avg$bme_humidity_percent_avg, hourly_avg$sensor1_pm2.5_ugm3_avg, use = "complete.obs")
cor_humidity_pm25_sensor2 <- cor(hourly_avg$bme_humidity_percent_avg, hourly_avg$sensor2_pm2.5_ugm3_avg, use = "complete.obs")

# Scatter plot with LOESS smoothing for temperature vs PM2.5 with correlation annotations
p_temp_pm25_loess <- ggplot(hourly_avg, aes(x = bme_temp_c_avg)) +
  geom_point(aes(y = sensor1_pm2.5_ugm3_avg, color = "Sensor 1 PM2.5"), size = 2, alpha = 0.6) +
  geom_point(aes(y = sensor2_pm2.5_ugm3_avg, color = "Sensor 2 PM2.5"), size = 2, alpha = 0.6) +
  geom_smooth(aes(y = sensor1_pm2.5_ugm3_avg, color = "Sensor 1 PM2.5"), method = "loess", se = FALSE) +
  geom_smooth(aes(y = sensor2_pm2.5_ugm3_avg, color = "Sensor 2 PM2.5"), method = "loess", se = FALSE, linetype = "dashed") +
  scale_color_manual(values = c("Sensor 1 PM2.5" = "#FF5733", "Sensor 2 PM2.5" = "#007AFF")) +
  labs(x = "Average Temperature (°C)", 
       y = "Average PM2.5 (µg/m³)",
       color = "PM2.5 Sensors",
       title = "PM2.5 Concentrations vs Temperature") +
  annotate("text", x = min(hourly_avg$bme_temp_c_avg, na.rm = TRUE), 
           y = max(hourly_avg$sensor1_pm2.5_ugm3_avg, na.rm = TRUE), 
           label = paste("r =", round(cor_temp_pm25_sensor1, 2)), 
           hjust = 0, color = "#FF5733", size = 5) +
  annotate("text", x = min(hourly_avg$bme_temp_c_avg, na.rm = TRUE), 
           y = max(hourly_avg$sensor2_pm2.5_ugm3_avg, na.rm = TRUE) - 5, 
           label = paste("r =", round(cor_temp_pm25_sensor2, 2)), 
           hjust = 0, color = "#007AFF", size = 5) +
  theme_minimal(base_size = 14)

# Print the temperature plot
print(p_temp_pm25_loess)

# Scatter plot with LOESS smoothing for humidity vs PM2.5 with correlation annotations
p_humidity_pm25_loess <- ggplot(hourly_avg, aes(x = bme_humidity_percent_avg)) +
  geom_point(aes(y = sensor1_pm2.5_ugm3_avg, color = "Sensor 1 PM2.5"), size = 2, alpha = 0.6) +
  geom_point(aes(y = sensor2_pm2.5_ugm3_avg, color = "Sensor 2 PM2.5"), size = 2, alpha = 0.6) +
  geom_smooth(aes(y = sensor1_pm2.5_ugm3_avg, color = "Sensor 1 PM2.5"), method = "loess", se = FALSE) +
  geom_smooth(aes(y = sensor2_pm2.5_ugm3_avg, color = "Sensor 2 PM2.5"), method = "loess", se = FALSE, linetype = "dashed") +
  scale_color_manual(values = c("Sensor 1 PM2.5" = "#FF5733", "Sensor 2 PM2.5" = "#007AFF")) +
  labs(x = "Average Humidity (%)", 
       y = "Average PM2.5 (µg/m³)",
       color = "PM2.5 Sensors",
       title = "PM2.5 Concentrations vs Humidity") +
  annotate("text", x = min(hourly_avg$bme_humidity_percent_avg, na.rm = TRUE), 
           y = max(hourly_avg$sensor1_pm2.5_ugm3_avg, na.rm = TRUE), 
           label = paste("r =", round(cor_humidity_pm25_sensor1, 2)), 
           hjust = 0, color = "#FF5733", size = 5) +
  annotate("text", x = min(hourly_avg$bme_humidity_percent_avg, na.rm = TRUE), 
           y = max(hourly_avg$sensor2_pm2.5_ugm3_avg, na.rm = TRUE) - 5, 
           label = paste("r =", round(cor_humidity_pm25_sensor2, 2)), 
           hjust = 0, color = "#007AFF", size = 5) +
  theme_minimal(base_size = 14)

# Print the humidity plot
print(p_humidity_pm25_loess)

#----- Ozone
# Create a LOESS plot
ggplot(hourly_avg, aes(x = timestamp_hour, y = ozone_ppm_avg)) +
  geom_point(color = "blue", alpha = 0.5) +  # Scatter plot of the data
  geom_smooth(method = "loess", color = "red", se = TRUE) +  # LOESS smoothing line with confidence interval
  labs(
    title = "Hourly Averaged Ozone Levels",
    x = "Timestamp",
    y = "Ozone (ppm)"
  ) +
  theme_minimal(base_size = 15) +  # Clean theme with larger text
  #scale_x_datetime(labels = date_format("%Y-%m-%d %H:%M"), date_breaks = "1 day") +  # Better date format for x-axis
  theme(axis.text.x = element_text(angle = 45, hjust = 1))  # Rotate x-axis labels for better readability

#Ultimately: make a census tract level of green to red
##Color codes (green to red): 0-7 (green), 7-8, 89, 9-10, 10-11, 11-12, 12-30 (red) with PM2.5 level
```

#PM Daily Avg Graphics
```{r}

#PM2.5
# Convert daily_avg to data.frame for ggplot
daily_avg <- as.data.frame(daily_avg)
daily_avg <- daily_avg %>%
  mutate(day = floor_date(timestamp_day, unit = "day"))


# Convert daily_avg to data.frame for ggplot
daily_avg <- as.data.frame(daily_avg)
daily_avg <- daily_avg %>%
  mutate(day = floor_date(timestamp_day, unit = "day"))

# Plot the Daily Averaged Data for PM2.5 with vibrant colors and reference line
pm25 <- ggplot(daily_avg, aes(x = timestamp_day)) +
  geom_line(aes(y = sensor1_pm2.5_ugm3_avg, color = "Sensor 1 PM2.5"), size = 1) +  # First y axis
  geom_line(aes(y = sensor2_pm2.5_ugm3_avg, color = "Sensor 2 PM2.5"), size = 1, linetype = "dashed") +  # Second y axis
  geom_hline(yintercept = 35, linetype = "dotted", color = "red") +  # Reference line
  scale_y_continuous(
    name = "Sensor 1 PM2.5 (µg/m³)", 
    sec.axis = sec_axis(~ ., name = "Sensor 2 PM2.5 (µg/m³)")
  ) +
  scale_color_manual(values = c("Sensor 1 PM2.5" = "#FF5733", "Sensor 2 PM2.5" = "#007AFF")) + 
  labs(x = "Timestamp (Daily Average)", 
       color = "PM2.5 Sensors", 
       title = "Daily Averaged PM2.5 Concentrations",
       subtitle = "Comparison of Sensor 1 and Sensor 2") +
  theme_minimal() +
  theme(
    axis.title.y = element_text(color = "#FF5733"),
    axis.title.y.right = element_text(color = "#007AFF"),
    legend.position = "top",
    plot.title = element_text(size = 18, face = "bold"),
    plot.subtitle = element_text(size = 14),
    axis.text = element_text(size = 12)
  )

# Print the plot
print(pm25)

#PM10
# Plot the Daily Averaged Data for PM10 with vibrant colors
pm10 <- ggplot(daily_avg, aes(x = timestamp_day)) +
  geom_line(aes(y = sensor1_pm10_ugm3_avg, color = "Sensor 1 PM10"), size = 1) +  # First y axis
  geom_line(aes(y = sensor2_pm10_ugm3_avg, color = "Sensor 2 PM10"), size = 1, linetype = "dashed") +  # Second y axis
  scale_y_continuous(
    name = "Sensor 1 PM10 (µg/m³)", 
    sec.axis = sec_axis(~ ., name = "Sensor 2 PM10 (µg/m³)")
  ) +
  scale_color_manual(values = c("Sensor 1 PM10" = "#FF5733", "Sensor 2 PM10" = "#007AFF")) +
  labs(x = "Timestamp (Daily Average)", 
       color = "PM10 Sensors", 
       title = "Daily Averaged PM10 Concentrations",
       subtitle = "Comparison of Sensor 1 and Sensor 2") +
  theme_minimal() +
  theme(
    axis.title.y = element_text(color = "#FF5733"),
    axis.title.y.right = element_text(color = "#007AFF"),
    legend.position = "top",
    plot.title = element_text(size = 18, face = "bold"),
    plot.subtitle = element_text(size = 14),
    axis.text = element_text(size = 12)
  )

# Display the PM10 plot
print(pm10)

#PM2.5 and Temp (C)

p_temp_pm25_heatmap <- ggplot(daily_avg, aes(x = bme_temp_c_avg, y = sensor1_pm2.5_ugm3_avg)) +
  geom_bin2d(bins = 30, aes(fill = ..count..)) +
  scale_fill_viridis_c(option = "plasma") +
  labs(x = "Average Temperature (°C)", 
       y = "Average PM2.5 (µg/m³)",
       fill = "Count",
       title = "Heatmap of PM2.5 Concentrations") +
  theme_minimal(base_size = 14)

p_temp_pm25_heatmap

#-------
# Load necessary library for correlation calculation
library(ggpubr)

# Calculate correlation between PM2.5 and temperature
cor_temp_pm25_sensor1 <- cor(daily_avg$bme_temp_c_avg, daily_avg$sensor1_pm2.5_ugm3_avg, use = "complete.obs")
cor_temp_pm25_sensor2 <- cor(daily_avg$bme_temp_c_avg, daily_avg$sensor2_pm2.5_ugm3_avg, use = "complete.obs")

# Calculate correlation between PM2.5 and humidity
cor_humidity_pm25_sensor1 <- cor(daily_avg$bme_humidity_percent_avg, daily_avg$sensor1_pm2.5_ugm3_avg, use = "complete.obs")
cor_humidity_pm25_sensor2 <- cor(daily_avg$bme_humidity_percent_avg, daily_avg$sensor2_pm2.5_ugm3_avg, use = "complete.obs")

# Scatter plot with LOESS smoothing for temperature vs PM2.5 with correlation annotations
p_temp_pm25_loess <- ggplot(daily_avg, aes(x = bme_temp_c_avg)) +
  geom_point(aes(y = sensor1_pm2.5_ugm3_avg, color = "Sensor 1 PM2.5"), size = 2, alpha = 0.6) +
  geom_point(aes(y = sensor2_pm2.5_ugm3_avg, color = "Sensor 2 PM2.5"), size = 2, alpha = 0.6) +
  geom_smooth(aes(y = sensor1_pm2.5_ugm3_avg, color = "Sensor 1 PM2.5"), method = "loess", se = FALSE) +
  geom_smooth(aes(y = sensor2_pm2.5_ugm3_avg, color = "Sensor 2 PM2.5"), method = "loess", se = FALSE, linetype = "dashed") +
  scale_color_manual(values = c("Sensor 1 PM2.5" = "#FF5733", "Sensor 2 PM2.5" = "#007AFF")) +
  labs(x = "Average Temperature (°C)", 
       y = "Average PM2.5 (µg/m³)",
       color = "PM2.5 Sensors",
       title = "PM2.5 Concentrations vs Temperature") +
  annotate("text", x = min(daily_avg$bme_temp_c_avg, na.rm = TRUE), 
           y = max(daily_avg$sensor1_pm2.5_ugm3_avg, na.rm = TRUE), 
           label = paste("r =", round(cor_temp_pm25_sensor1, 2)), 
           hjust = 0, color = "#FF5733", size = 5) +
  annotate("text", x = min(daily_avg$bme_temp_c_avg, na.rm = TRUE), 
           y = max(daily_avg$sensor2_pm2.5_ugm3_avg, na.rm = TRUE) - 5, 
           label = paste("r =", round(cor_temp_pm25_sensor2, 2)), 
           hjust = 0, color = "#007AFF", size = 5) +
  theme_minimal(base_size = 14)

# Print the temperature plot
print(p_temp_pm25_loess)

# Scatter plot with LOESS smoothing for humidity vs PM2.5 with correlation annotations
p_humidity_pm25_loess <- ggplot(daily_avg, aes(x = bme_humidity_percent_avg)) +
  geom_point(aes(y = sensor1_pm2.5_ugm3_avg, color = "Sensor 1 PM2.5"), size = 2, alpha = 0.6) +
  geom_point(aes(y = sensor2_pm2.5_ugm3_avg, color = "Sensor 2 PM2.5"), size = 2, alpha = 0.6) +
  geom_smooth(aes(y = sensor1_pm2.5_ugm3_avg, color = "Sensor 1 PM2.5"), method = "loess", se = FALSE) +
  geom_smooth(aes(y = sensor2_pm2.5_ugm3_avg, color = "Sensor 2 PM2.5"), method = "loess", se = FALSE, linetype = "dashed") +
  scale_color_manual(values = c("Sensor 1 PM2.5" = "#FF5733", "Sensor 2 PM2.5" = "#007AFF")) +
  labs(x = "Average Humidity (%)", 
       y = "Average PM2.5 (µg/m³)",
       color = "PM2.5 Sensors",
       title = "PM2.5 Concentrations vs Humidity") +
  annotate("text", x = min(daily_avg$bme_humidity_percent_avg, na.rm = TRUE), 
           y = max(daily_avg$sensor1_pm2.5_ugm3_avg, na.rm = TRUE), 
           label = paste("r =", round(cor_humidity_pm25_sensor1, 2)), 
           hjust = 0, color = "#FF5733", size = 5) +
  annotate("text", x = min(daily_avg$bme_humidity_percent_avg, na.rm = TRUE), 
           y = max(daily_avg$sensor2_pm2.5_ugm3_avg, na.rm = TRUE) - 5, 
           label = paste("r =", round(cor_humidity_pm25_sensor2, 2)), 
           hjust = 0, color = "#007AFF", size = 5) +
  theme_minimal(base_size = 14)

# Print the humidity plot
print(p_humidity_pm25_loess)

#----- Ozone
# Create a LOESS plot
ggplot(daily_avg, aes(x = timestamp_day, y = ozone_ppm_avg)) +
  geom_point(color = "blue", alpha = 0.5) +  # Scatter plot of the data
  geom_smooth(method = "loess", color = "red", se = TRUE) +  # LOESS smoothing line with confidence interval
  labs(
    title = "Daily Averaged Ozone Levels",
    x = "Timestamp",
    y = "Ozone (ppm)"
  ) +
  theme_minimal(base_size = 15) +  # Clean theme with larger text
  #scale_x_datetime(labels = date_format("%Y-%m-%d %H:%M"), date_breaks = "1 day") +  # Better date format for x-axis
  theme(axis.text.x = element_text(angle = 45, hjust = 1))  # Rotate x-axis labels for better readability

```

#AQI Graphs
```{r}

# Function to calculate AQI from PM2.5
calculate_aqi_pm25 <- function(pm25) {
  if (pm25 <= 9.0) {
    return(50 * pm25 / 9)  # Good
  } else if (pm25 <= 35.4) {
    return(50 + (100 - 50) * (pm25 - 9.1) / (35.4 - 9.1))  # Moderate
  } else if (pm25 <= 55.4) {
    return(100 + (150 - 100) * (pm25 - 35.5) / (55.4 - 35.5))  # Unhealthy for Sensitive Groups
  } else if (pm25 <= 125.4) {
    return(150 + (200 - 150) * ((pm25 - 55.5)^3) / ((125.4 - 55.5)^3))  # Unhealthy
  } else if (pm25 <= 225.4) {
    return(200 + (300 - 200) * ((pm25 - 125.5)^3) / ((225.4 - 125.5)^3))  # Very Unhealthy
  } else {
    return(301)  # Hazardous
  }
}

# Function to calculate AQI from PM10
calculate_aqi_pm10 <- function(pm10) {
  if (pm10 <= 54) {
    return(50 * pm10 / 54)  # Good
  } else if (pm10 <= 154) {
    return(50 + (100 - 50) * (pm10 - 55) / (154 - 55))  # Moderate
  } else if (pm10 <= 254) {
    return(100 + (150 - 100) * (pm10 - 155) / (254 - 155))  # Unhealthy for Sensitive Groups
  } else if (pm10 <= 354) {
    return(150 + (200 - 150) * (pm10 - 255) / (354 - 255))  # Unhealthy
  } else if (pm10 <= 424) {
    return(200 + (300 - 200) * (pm10 - 355) / (424 - 355))  # Very Unhealthy
  } else {
    return(301)  # Hazardous
  }
}


# Function to calculate AQI from Ozone 
calculate_aqi_o3 <- function(o3) {
  if (is.na(o3)) {
    return(NA)  # Return NA if o3 is NA
  } else if (o3 <= 0.054) {
    return(50 * o3 / 0.054)  # Good
  } else if (o3 <= 0.070) {
    return(50 + (100 - 50) * (o3 - 0.055) / (0.070 - 0.055))  # Moderate
  } else if (o3 <= 0.085) {
    return(100 + (150 - 100) * (o3 - 0.071) / (0.085 - 0.071))  # Unhealthy for Sensitive Groups
  } else if (o3 <= 0.105) {
    return(150 + (200 - 150) * (o3 - 0.086) / (0.105 - 0.086))  # Unhealthy
  } else if (o3 <= 0.200) {
    return(200 + (300 - 200) * (o3 - 0.106) / (0.200 - 0.106))  # Very Unhealthy
  } else {
    return(301)  # Hazardous
  }
}

#Calculate the 8-hour rolling average for ozone

# Apply the AQI calculation to the dataset
daily_avg$aqi_sensor1_pm25 <- sapply(daily_avg$sensor1_pm2.5_ugm3_avg, calculate_aqi_pm25)
daily_avg$aqi_sensor2_pm25 <- sapply(daily_avg$sensor2_pm2.5_ugm3_avg, calculate_aqi_pm25)
daily_avg$aqi_sensor1_pm10 <- sapply(daily_avg$sensor1_pm10_ugm3_avg, calculate_aqi_pm10)
daily_avg$aqi_sensor2_pm10 <- sapply(daily_avg$sensor2_pm10_ugm3_avg, calculate_aqi_pm10)
#daily_avg$aqi_o3 <- sapply(daily_avg$ozone_8hr_avg, calculate_aqi_o3)

# Define the corrected color palette
aqi_colors <- c(
  "Good" = "#66c2a5",  # Green
  "Moderate" = "#ffd972",  # Yellow
  "Unhealthy for Sensitive Groups" = "#fc8d62",  # Orange
  "Unhealthy" = "#d90429",  # Red
  "Very Unhealthy" = "#a01a7d"  # Purple
)

# Convert the category column to a factor with the desired order
aqi_categories <- data.frame(
  category = c("Good", "Moderate", "Unhealthy for Sensitive Groups", "Unhealthy", "Very Unhealthy"),
  ymin = c(0, 50, 100, 150, 200),  # Corresponding lower bounds
  ymax = c(50, 100, 150, 200, 300)  # Corresponding upper bounds
)


# Create the plot
pm_aqi_plot <- ggplot() +
  # Add shaded regions for AQI categories
  geom_rect(data = aqi_categories, aes(xmin = min(daily_avg$timestamp_day), xmax = max(daily_avg$timestamp_day), ymin = ymin, ymax = ymax, fill = category), alpha = 0.55) +
  # Add line plots for AQI with refined colors
  geom_line(data = daily_avg, aes(x = timestamp_day, y = aqi_sensor1_pm25, color = "PM2.5 Sensor 1 AQI"), size = 1.2) +
  geom_line(data = daily_avg, aes(x = timestamp_day, y = aqi_sensor2_pm25, color = "PM2.5 Sensor 2 AQI"), size = 1.2, linetype = "dotted") +
  geom_line(data = daily_avg, aes(x = timestamp_day, y = aqi_sensor1_pm10, color = "PM10 Sensor 1 AQI"), size = 1.2) +
  geom_line(data = daily_avg, aes(x = timestamp_day, y = aqi_sensor2_pm10, color = "PM10 Sensor 2 AQI"), size = 1.2, linetype = "dotted") +
  #geom_line(data = daily_avg, aes(x = timestamp_day, y = aqi_o3, color = "Ozone 8-hour AQI"), size = 1.2) +
  
  # Define colors for AQI categories and sensor lines
  scale_fill_manual(
    name = "AQI Categories",
    values = aqi_colors
  ) +
  scale_color_manual(values = c(
  #"PM2.5 Sensor 1 AQI" = viridis(6)[1], 
  "PM2.5 Sensor 1 AQI" = "#4F2C8C",
  "PM2.5 Sensor 2 AQI" = viridis(6)[2], 
  "PM10 Sensor 1 AQI" = "#00B8B5",
  #"PM10 Sensor 1 AQI" = viridis(6)[3], 
  "PM10 Sensor 2 AQI" = "#00B8B5"
  #"Ozone 8-hour AQI" = "#4AE114"  # A44200
    )) +
  
  # Enhance labeling and title presentation
  labs(x = "Timestamp", 
       y = "AQI", 
       color = "Sensors") +
  
  # Minimal theme with slight adjustments
  theme_minimal() +
  theme(
    legend.position = "right",
    plot.title = element_text(size = 20, face = "bold"),
    plot.subtitle = element_text(size = 16),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12),
    panel.grid.minor = element_blank()  # Remove minor gridlines for clarity
  ) +
  guides(
    fill = guide_legend(override.aes = list(alpha = 0.55))  # Ensure transparency is shown in the legend
  )

# Display the plot
print(pm_aqi_plot)

#--------------
# Calculate the 8-hour rolling average for ozone based on hourly data
hourly_avg$ozone_8hr_avg <- rollmean(hourly_avg$ozone_ppm_avg, 8, align = "right", fill = NA)

#8-hour ozone
hourly_avg$aqi_o3 <- sapply(hourly_avg$ozone_8hr_avg, calculate_aqi_o3)

# Create the plot
pm_o3_plot <- ggplot() +
  # Add shaded regions for AQI categories
  geom_rect(data = aqi_categories, aes(xmin = min(hourly_avg$timestamp_hour), xmax = max(hourly_avg$timestamp_hour), ymin = ymin, ymax = ymax, fill = category), alpha = 0.55) +
  # Add line plots for AQI with refined colors
  geom_line(data = hourly_avg, aes(x = timestamp_hour, y = aqi_o3, color = "Ozone 8-hour AQI"), size = 1.2) +
  
  # Define colors for AQI categories and sensor lines
  scale_fill_manual(
    name = "AQI Categories",
    values = aqi_colors
  ) +
  scale_color_manual(values = c(
  "Ozone 8-hour AQI" = "#006400"  # A44200
    )) +
  
  # Enhance labeling and title presentation
  labs(x = "Timestamp", 
       y = "AQI", 
       color = "Sensors") +
  
  # Minimal theme with slight adjustments
  theme_minimal() +
  theme(
    legend.position = "right",
    plot.title = element_text(size = 20, face = "bold"),
    plot.subtitle = element_text(size = 16),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12),
    panel.grid.minor = element_blank()  # Remove minor gridlines for clarity
  ) +
  guides(
    fill = guide_legend(override.aes = list(alpha = 0.55))  # Ensure transparency is shown in the legend
  )
pm_o3_plot
 #How to interpret: highest AQI of the three pollutants would be the index for the day/hour.


```

#Next

  - make a census tract map with aqi 
  - when neighborhoods come in, do autocorrelation analysis + code
  - when all data comes in across all neighborhoods, do spatial econ regression