report.Rmd

---
title: "South Gate PurpleAir Data Report"
author: "~~~~~"
output: pdf_document
params:
  answer: NA
  breaks: NA
  d1: NA
  d2: NA
  d3: NA
  d4: NA
  dts1: NA
  dts2: NA
  dts3: NA
  dts4: NA
  daily: NA
  down: NA
  file2: NA
  hour: NA
  highlow: NA
  matching: NA
  over: NA
  overEPA: NA
  PAfull: NA
  PAhourly: NA
  sensor: NA
  sensors: NA
  summary: NA
  under: NA
---


```{r setup, include=FALSE, echo=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(shiny)
library(ggplot2)
library(dplyr)
library(plotly)
library(tidyverse)
library(zoo)
library(gridExtra)
library(PurpleAirCEHAT)
library(lubridate)

date1 <- params$d1
dates1 <- c(params$dts1)
dates2 <- c(params$dts2)
dates3 <- params$dts3
dates4 <- params$dts4
answer <- params$answer

matchingDays <- params$matching


PAfull <- params$PAfull
PAhourly <- params$PAhourly
PAhi_lo <- params$highlow
downSensors <- params$down
avgSG <- params$summary
dailySG <- params$daily
sensor <- params$sensor

selectedSensors <- params$sensors
readings_underCT <- params$under
readings_overCT <- params$over


sg.city <- PurpleAirCEHAT::southgate()


EPAcols <- c("Good"="#00e400", "Moderate"="#ffff00","Unhealthy for Sensitive Groups" = "#ff7e00", "Unhealthy" = "#ff0000", "Very Unhealthy" = "#8f3f97", "Hazardous" = "#7e0023")

wessy_pal <- c("high"="#C93312","low"="#899DA4")


Krigehourly <- PAhourly %>% dplyr::filter(PAhourly$timestamp == as_datetime(date1)+ lubridate::hours(params$hour))

autoDF <- data.frame(PurpleAirCEHAT::krigePA(Krigehourly, as_datetime(date1)+ lubridate::hours(params$hour)))

```


## Overview 


The following sections show the plots in each component of our interactive data analysis application.We wanted to ensure that we provided users with the most basic summary calculations for air quality data, including means, medians, maximums, minimums, and diurnal patterns as well as results that are directly comparable to other AB617 communities with active AQMD PM2.5 sensors. We also wanted to allow users to monitor the functionality of their sensor network, through summative information about particular sensors and interpolation calculations.


#### Data Cleaning
The first step of processing data downloaded from the website is formatting it for data analysis. We refer to this process as 'cleaning' data. The format of the raw sensor data changed after PurpleAir deprecated its API. So, we use two different data cleaning processes, one for the old data (data from before March 30, 2021) and one for the new data (all data from after March 30,2021).

The function that cleans the old data assumes that the comma-separated values (csv) file has two PM2.5 reading channels, as well as timestamp, humidity, latitude, and longitude columns. We correct the timestamp (so that it is adjusted to the America/Los Angeles time zone) and get rid of any values that exceed the standard EPA thresholds. Then, we use a correction factor to adjust the data so that it meets EPA standards.

The function that cleans the new data assumes the CSV file has a single PM2.5 column, as well as columns for timestamp, sensor names, latitude, longitude, humidity, and mean percent difference. To clean the new data, we adjust the timestamp to the America/Los Angeles time zone, and make sure to remove any negative PM2.5 readings in the data.

The PurpleAir sensors used in generating this report are listed below:
```{r sensors, echo=FALSE}
cat(selectedSensors, sep = "\n")
```



## Summary Statistics 


For the summary statistics calculations, we refer to the general hourly PM2.5 for South Gate PurpleAir sensors. The general hourly readings are the aggregated PM2.5 values by hour, which is an average of each sensor's hourly PM2.5 readings.

```{r past3, echo=FALSE}

Average <- mean(avgSG$average_PM2.5)

historical <- ggplot(avgSG[avgSG$timestamp <= as_datetime(params$d2) & avgSG$timestamp >= as_datetime(params$d2)-lubridate::days(3),], aes(x=as_datetime(timestamp),y=average_PM2.5, group=category))+
    geom_hline(aes(yintercept = Average), color="blue", linetype="dashed")+
    geom_col(aes(fill=category),col=1, lwd=0.5)+
    labs(x = "Day", y = "PM2.5 (micrograms/m3)") +
    ggtitle("Hourly Average PM2.5")+
    theme_minimal()+
    scale_fill_discrete(type=EPAcols)+
    theme(axis.text.x = element_text(angle = 75, hjust=1))


print(historical)
```


Similarly, we calculate the daily PM2.5 averages for South Gate PurpleAir sensors which are aggregated PM2.5 every 24 hours.

```{r pastMonth, echo=FALSE}
Average <- mean(dailySG$average_PM2.5)

historical <- ggplot(dailySG[dailySG$day <= as.Date(params$d3) & dailySG$day >= as.Date(params$d3)-lubridate::days(30),], aes(x=day,y=average_PM2.5, group=category))+
    geom_hline(aes(yintercept = Average), color="blue", linetype="dashed")+
    geom_col(aes(fill=category),col=1, lwd=0.5)+
    labs(x = "Day", y = "PM2.5 (micrograms/m3)") +
    ggtitle("Daily Average PM2.5")+
    theme_minimal()+
    scale_fill_discrete(type=EPAcols)+
    theme(axis.text.x = element_text(angle = 75, hjust=1))


print(historical)

```


Below we use a lollipop chart to report the highs and lows of general PM2.5 values.

```{r highlow, echo=FALSE}

HighAverage <- mean(dailySG$max)

LowAverage <- mean(dailySG$min)

hilo <- ggplot(data=dailySG[dailySG$day <= as.Date(params$d4) & dailySG$day >= as.Date(params$d4)-lubridate::days(30),],   aes(x=day, group = day)) +
        geom_segment(aes(x=day, xend=day, y=min,yend=max),lwd=1)+
        geom_point(aes(y=max, col="high"), size=5) +
        geom_point(aes(y=min, col="low"), size=5)+
        geom_hline(aes(yintercept = HighAverage), color="#C93312", linetype="dashed")+
        geom_hline(aes(yintercept = LowAverage), color="#899DA4", linetype="dashed")+
        labs(x = 'Day', y = 'PM2.5 (micrograms/m3)') +
        scale_colour_manual(name="Type", values=wessy_pal, guide = guide_legend(override.aes=aes(fill=NA)) ) +
        ggtitle("Daily Highs and Lows") +
        theme_minimal()

print(hilo)

```



```{r overEPA, echo=FALSE}

# numDays <- length(unique(lubridate::mday(summarySG$timestamp)))
# 
# # finding the days over EPA threshold
# overEPA <- params$overEPA
# 
# 
# if(length(overEPA$timestamp) > 0) {
#    ourData <- PurpleAirCEHAT::overEPA_hist(overEPA,numDays)
#    epahist <- ggplot(ourData, aes(x=day,y=freq)) +
#    geom_bar(position="dodge", stat="identity") +
#    ggtitle("Days over EPA threshold in South Gate")}
# 
# 
# print(epahist)

```


The following two graphs display 8 hour averages for two distinct, specified time periods.

```{r avgs, echo=FALSE, fig.show="hold", out.width="50%"}
par(mar = c(4, 4, .1, .1))

Title1 <- paste("8 Hour Averages From", paste(as.character(dates1), collapse = " to "))

avgs1 <- zoo::rollmean(avgSG$average_PM2.5[lubridate::date(avgSG$timestamp) >= toString(dates1[1]) & lubridate::date(avgSG$timestamp) <= toString(dates1[2])], k=8)

        
        
Title2 <- paste("8 Hour Averages From", paste(as.character(dates2), collapse = " to "))

avgs2 <- zoo::rollmean(avgSG$average_PM2.5[lubridate::date(avgSG$timestamp) >= toString(dates2[1]) & lubridate::date(avgSG$timestamp) <= toString(dates2[2])], k=8)

paste(paste(Title1,"(left)", sep=" "), paste(Title2,"(right)", sep=" "), sep=" ")

plot(avgs1,type="l", xlab = "Time", ylab = "PM2.5 (micrograms/m3)")
plot(avgs2,type="l", xlab = "Time", ylab = "PM2.5 (micrograms/m3)")

```




The following three graphs display the diurnal patterns for average, peak, the range of PM2.5 values over a 24-hour period.

```{r diAVG, echo=FALSE}

dates <- params$dts4

diurnalR <- aggregate(cbind(average_PM2.5) ~ hour, data = avgSG, FUN= function(x) {round(mean(x),2)} )

addDiurnal <- aggregate(cbind(average_PM2.5) ~ hour,
                        data = avgSG[as.Date(avgSG$timestamp, tz = "America/Los_Angeles") >= as.Date(dates4[1], tz = "America/Los_Angeles") &
                                         as.Date(avgSG$timestamp, tz = "America/Los_Angeles") <= as.Date(dates4[2], tz = "America/Los_Angeles"),],
                        FUN= function(x) {round(mean(x),2)} )

#creating a dynamic scale for the plot
min1 <- min(diurnalR)
min2 <- min(addDiurnal)

max1 <- max(diurnalR)
max2 <- max(addDiurnal)

plot(diurnalR,type="o", lwd=1.5, main = "Average PM2.5 Values", xlab= "Hour", ylab="PM2.5 (microgams/m3)", ylim=c(min(min1,min2), max(max1,max2)))
lines(addDiurnal,type="o", lwd=1.5, col="blue")
grid()

```



```{r diMax, echo=FALSE}

diurnalR <- aggregate(cbind(max) ~ hour, data = avgSG, FUN= function(x) {round(mean(x),2)} )

addDiurnal <- aggregate(cbind(max) ~ hour,
                        data = avgSG[lubridate::date(avgSG$timestamp) >= toString(dates4[1]) &
                                         lubridate::date(avgSG$timestamp) <= toString(dates4[2]),],
                        FUN= function(x) {round(mean(x),2)} )

#creating a dynamic scale for the plot
min1 <- min(diurnalR)
min2 <- min(addDiurnal)

max1 <- max(diurnalR)
max2 <- max(addDiurnal)

plot(diurnalR,type="o", lwd=1.5, main = "Peak PM2.5 Values", xlab= "Hour", ylab="PM2.5 (micrograms/m3)", ylim=c(min(min1,min2), max(max1,max2)))
lines(addDiurnal, type="o", lwd=1.5, col="blue")
grid()

```



```{r diRange, echo=FALSE}

diurnalR <- aggregate(cbind(range) ~ hour, data = avgSG, FUN= function(x) {round(mean(x),2)} )

addDiurnal <- aggregate(cbind(range) ~ hour,
                        data = avgSG[lubridate::date(avgSG$timestamp) >= toString(dates4[1]) &
                                         lubridate::date(avgSG$timestamp) <= toString(dates4[2]),],
                        FUN= function(x) {round(mean(x),2)} )

#creating a dynamic scale for the plot
min1 <- min(diurnalR)
min2 <- min(addDiurnal)

max1 <- max(diurnalR)
max2 <- max(addDiurnal)

plot(diurnalR,type="o", lwd=1.5, main = "Range of PM2.5 Values", xlab= "Hour", ylab="PM2.5 (micrograms/m3)", ylim=c(min(min1,min2), max(max1,max2)))
lines(addDiurnal,type="o", lwd=1.5, col="blue")
grid()

```


This chart displays the density of peak PM2.5 values over a specified 24-hour period.
```{r density, echo=FALSE}

suppressWarnings({   
dens <- ggplot(PAhi_lo[PAhi_lo$type == "high" & lubridate::date(PAhi_lo$timestamp) >= toString(dates3[1]) & lubridate::date(PAhi_lo$timestamp) <= toString(dates3[2]),], aes(x=hour, group = timeofday))+
  geom_histogram(aes(y=after_stat(count/nrow(PAhi_lo[PAhi_lo$type == "high",])), color=timeofday, fill=timeofday), alpha=0.7, stat="count", bins=4, lwd=1)+
  geom_density(alpha = 0.2, fill = "grey")+
  labs(x = "Hour", y = "Density") +
  ggtitle("Density of Peak PM2.5 Values Over 24-hour Period")+
  scale_color_brewer(palette="Set1")+
  scale_fill_brewer(palette="Set1")+
  theme_minimal()

print(dens)
})

```

## Sensor Summaries

The Sensor Summaries plots allow users to inspect the consistency of each sensor, identify which sensors typically report higher/lower values of PM2.5, and identify which sensors typically report outliers. These plots are particularly helpful in identifying any geographically specific air quality patterns that may affect a sensor's typical functionality, and also where air quality in South Gate is better (or worse). 



#### Network-Based Observations

The plots below provide a holistic inspection of the sensor network, showing where sensors report higher (or lower) values, which sensors typically report outliers, and finally, how often each of the sensors go down.


```{r catsSensors, echo=FALSE}

suppressWarnings({
plot <- ggplot(PAhourly, aes(x=names, group=category)) +
            geom_histogram(aes(y=after_stat(count), fill=category), position = position_stack(reverse = TRUE), stat="count", lwd=0.5, col="black")+
            labs(x = "Sensors", y = "Count") +
            ggtitle("Category of Readings by Sensor")+
            theme_minimal()+
            scale_fill_discrete(type=EPAcols)+
            theme(axis.text.x = element_text(angle = 75, hjust=1))

print(plot)
})
```



```{r downDays, echo=FALSE}

down <- ggplot(downSensors, aes(y=numDownDays))+
    geom_col(aes(x= names, y=numDownDays), lwd=1)+
    labs(x = "Sensors", y = "Number of Days") +
    ggtitle("Number of 'Down' Days")+
    coord_cartesian(ylim= c(0,max(downSensors[,2])))+
    scale_color_brewer(palette="RdYlBu")+
    scale_fill_brewer(palette="RdYlBu")+
    theme_minimal()+
    theme(axis.text.x = element_text(angle = 45, hjust=1))

print(down)

```



The following charts and spatial plots illustrate which sensors typically report higher or lower values. They report on the number of readings over/under the median that each sensor records, and the magnitude of this over/under in terms of percent difference. The median here is calculated at each hour for all the sensors in the network, and the readings are all taken to exclusively be over/under the median. The count of readings over the median for each sensor are further striated into categories of incremental percent difference that inform about how much greater, or lesser, than the median the observed readings are. This chart is useful for both identifying which sensors are typically reporting higher values, and of course, which sensors are consistently reporting outliers.



```{r overMedian, echo=FALSE}


suppressWarnings({
pct_diffCols <- c("0-15%"="#FEEDDE", "15%-50%"="#FDBE85","50%-100%" = "#FD8D3C", "200% or more"="#E6550D","total readings"="#80808050")

plot <- ggplot(readings_overCT, aes(x=names)) +
    geom_col(aes(y=total_readings, fill= "total readings"), col=1) +
    geom_col(aes(y=`0-15%`, fill="0-15%"), col=1) +
    geom_col(aes(y=`15%-50%`, fill="15%-50%"), col=1) +
    geom_col(aes(y=`50%-100%`, fill="50%-100%"), col=1) +
    geom_col(aes(y=`above_200%`, fill="200% or more"), col=1) +
    labs(x = "Sensors", y = "Count") +
    ggtitle("Readings Over Median by Sensor")+
    theme_minimal()+
    scale_colour_manual(name="Values",values=pct_diffCols, guide = guide_legend(override.aes=aes(fill=NA)) ) +
    scale_fill_manual(name="Values",values=pct_diffCols) +
    theme(axis.text.x = element_text(angle = 75, hjust=1))

print(plot)
})

```



```{r overMap, echo=FALSE}

suppressWarnings({
k <- ggplot(readings_overCT, aes(longitude, latitude)) +
            geom_path(data = sg.city, aes(long, lat, group=id), color='black')+
            geom_point(aes(size= count/total_readings, color=count/total_readings)) +
            scale_color_gradient(low='gold', high='red', name= "Normalized High Readings") +
            guides(size=FALSE) +
            xlim(-118.2325,-118.155) +
            ylim(33.91029, 33.96837)+
            ggtitle("Readings Over Median by Sensor")+
            geom_text(aes(label=names), check_overlap = F, show.legend = F, size = 3, vjust = 2)+
            theme_minimal()

print(k)
})
```



```{r underMedian, echo=FALSE}

suppressWarnings({
pct_diffCols <- c("0-15%"="#EDF8FB", "15%-50%"="#B3CDE3","50%-100%" = "#8C96C6", "200% or more"="#88419D","total readings"="#80808050")

plot <- ggplot(readings_underCT, aes(x=names)) +
    geom_col(aes(y=total_readings, fill= "total readings"), col=1) +
    geom_col(aes(y=`0-15%`, fill="0-15%"), col=1) +
    geom_col(aes(y=`15%-50%`, fill="15%-50%"), col=1) +
    geom_col(aes(y=`50%-100%`, fill="50%-100%"), col=1) +
    geom_col(aes(y=`above_200%`, fill="200% or more"), col=1) +
    labs(x = "Sensors", y = "Count") +
    ggtitle("Readings Under Median by Sensor") +
    theme_minimal()+
    scale_colour_manual(name="Values",values=pct_diffCols, guide = guide_legend(override.aes=aes(fill=NA)) ) +
    scale_fill_manual(name="Values",values=pct_diffCols) +
    theme(axis.text.x = element_text(angle = 75, hjust=1))

print(plot)
})
```




```{r underMap, echo=FALSE}
suppressWarnings({
k <- ggplot(readings_underCT, aes(longitude, latitude)) +
    geom_path(data = sg.city, aes(long, lat, group=id), color='black')+
    geom_point(aes(size= count/total_readings, color=count/total_readings)) +
    scale_color_gradient(low='violet', high='blue', name= "Normalized Low Values") +
    xlim(-118.2325,-118.155) +
    ylim(33.91029, 33.96837)+
    guides(size=FALSE) +
    ggtitle("Readings Under Median by Sensor")+
    geom_text(aes(label=names), check_overlap = F, show.legend = F, size = 3, vjust = 2)+
    theme_minimal()


print(k)
})
```






#### Single-sensor Observations

Here we have plots and observations for an individual sensor. There are visualizations for high and low values, EPA categories by percentage, and mean percent difference data (for data downloaded after March, 30, 2021).

```{r catsbysensor, echo=FALSE}

suppressWarnings({
print(ggplot(PAhourly[PAhourly$names == sensor,], aes(x=category, group=category)) +
          geom_histogram(aes(y=after_stat(count/nrow(PAhourly[PAhourly$names == sensor,])), fill=category), stat="count", lwd=0.5, col="black")+
          labs(x = "Sensors", y = "Count") +
          ggtitle(paste("Category of Readings for", sensor)) +
          theme_minimal()+
          scale_fill_discrete(type=EPAcols)+
          theme(axis.text.x = element_text(angle = 75, hjust=1)))
})

```


```{r pctDiff, echo=FALSE}
if (answer == "Y"){
   title <- paste("Mean Percent Difference for", sensor, sep=" ")
   plot(x=PAfull$timestamp[PAfull$names==sensor], y=PAfull$percent.diff[PAfull$names==sensor], xlab = "Time", ylab = "Percent Difference", main = title, ylim = c(0,100), type="l")
   } else { 
      df <- as.data.frame(x=c(),percent.diff=c())
      plot(df$x,df$percent.diff, xlim = c(0,10),ylim =c(0,100))}

```


```{r highlowSensor, echo=FALSE}

hilo <- ggplot(data=PAhi_lo[PAhi_lo$names == sensor,], aes(x=date, y=PM2.5, group = date)) +
        geom_line(lwd=1)+
        geom_point(data=PAhi_lo[PAhi_lo$type == "high"  & PAhi_lo$names == sensor, ],
                   aes(x=date, y=PM2.5, group = type, col="high"), size=5)+
        geom_point(data=PAhi_lo[PAhi_lo$type == "low" & PAhi_lo$names == sensor, ],
                   aes(x=date, y=PM2.5, group = type, col="low"), size=5)+
        labs(x = 'Day', y = 'PM2.5 (μg/m3)') +
        scale_colour_manual(name="Type",values=wessy_pal, guide = guide_legend(override.aes=aes(fill=NA)) ) +
        scale_fill_manual(name="Type",values=wessy_pal) +
        ggtitle(sensor) +
        theme_minimal()

print(hilo)


```




## Interpolation and Sensor Placement

In the following plots, one can find the three outputs of the ordinary kriging results. This includes the predictions along with the variance and corresponding standard deviation of those predictions. It is important to note that better coverage reduces variance and standard deviation. So, if the user wants to know where to place more sensors, they should pick places with a higher variance/standard deviation to improve the accuracy of estimation and make the network more effective.

```{r prediction, echo=FALSE}

names(autoDF)[3] <- "predicted"

autoPred <- ggplot() + geom_tile(autoDF, mapping = aes(x,y,fill=predicted), alpha=0.90) +
    geom_point(Krigehourly[,c('PM2.5',"longitude","latitude")], color ="black", size=2, pch=21, mapping = aes(longitude, latitude, fill=PM2.5), inherit.aes = TRUE) +
    coord_equal() +
    scale_fill_continuous(type = "viridis") +
    labs(x = "longitude", y="latitude")+
    theme_bw() +
    ggtitle("Ordinary Kriging PM2.5 Predictions")

print(autoPred)

```


```{r variance, echo=FALSE}

names(autoDF)[4] <- "variance"

autoVars <- ggplot() + geom_tile(autoDF, mapping = aes(x,y,fill=variance), alpha=0.90) +
    geom_point(Krigehourly[,c('PM2.5',"longitude","latitude")], color ="black", size=2, pch=21, mapping =aes(longitude, latitude, fill=PM2.5), inherit.aes = TRUE) +
    coord_equal() +
    scale_fill_continuous(type = "viridis") +
    labs(x = "longitude", y="latitude")+
    theme_bw() +
    ggtitle("Ordinary Kriging PM2.5 Variance")

print(autoVars)

```


```{r stdev, echo=FALSE}

names(autoDF)[5] <- "standard deviation"

autoStDev <- ggplot() + geom_tile(autoDF, mapping = aes(x,y,fill=`standard deviation`), alpha=0.90) +
    geom_point(Krigehourly[,c('PM2.5',"longitude","latitude")], color ="black", size=2, pch=21 , mapping =aes(longitude, latitude, fill=PM2.5), inherit.aes = TRUE) +
    coord_equal() +
    scale_fill_continuous(type = "viridis") +
    labs(x = "longitude", y="latitude")+
    theme_bw() +
    ggtitle("Ordinary Kriging PM2.5 Standard Deviation")

print(autoStDev)

```