---
title: "Analysis of Ozone Season NO~X~ Emissions Data for Coal-Fired EGUs in Four Mid-Atlantic States"
author:
date: "`r format(Sys.time(), '%B %d, %Y')`"
output:
  bookdown::pdf_document2: 
    number_sections: true
    fig_caption: yes
  bookdown::word_document2: 
    number_sections: true
    fig_caption: yes
  bookdown::html_document2:
    df_print: paged
---

# Introduction

The Ozone Transport Commission (OTC) is recommending that EPA require Pennsylvania to set daily nitrogen oxides (NO~X~) emissions limits during the ozone season for coal-fired electric generating units (EGUs) with already-installed selective catalytic reduction (SCR) or selective non-catalytic reduction (SNCR) controls. Further, the OTC recommends that these limits be as stringent as limits set by Delaware, Maryland, and New Jersey.

This code accompanies the Technical Support Document (TSD) titled "Analysis of Ozone Season NO~X~ Emissions Data for Coal-Fired EGUs in Four Mid-Atlantic States" for the "Notice of Public Hearing and Supplemental Information–Ozone Transport Commission Recommendation that EPA Require Daily Limits for Emissions of Nitrogen Oxides from Certain Sources in Pennsylvania". The purpose of the TSD is to explore recent NO~X~ emission rates and trends at coal-fired EGUs equipped with SCR or SNCR in Pennsylvania, Delaware, Maryland, and New Jersey. For comparison, coal-fired EGUs not equipped with SCR or SNCR are also included in the analysis. The analysis primarily focuses on ozone season data for 2009-2019. This analysis relies on EPA Clean Air Markets Division's (CAMD's) [Power Sector Emissions Data](link: https://www.epa.gov/airmarkets/power-sector-emissions-data).

This code will import and prepare the relevant data as well as generate the unit-level plots of ozone season NO~X~ emission rates, NO~X~ emissions, and operating hours that appear in Section 3 of the TSD. The data and code are provided for transparency and replication of the data analysis.

```{r setup, include = FALSE, message = FALSE}
#clear out the global environment
rm(list = ls())

#libraries
library(readxl)
library(dplyr)
library(tidyverse)
library(ggplot2)
library(grid)
library(gridExtra)
library(knitr)
library(scales)
library(bookdown)

#set working directory
#This path specification will ensure that any user can knit the notebook so long as the Rmd document is located in the same directory as the 'data' subfolder.
setwd("./")

#prevent warnings from showing up in output
knitr::opts_chunk$set(message=FALSE, warning=FALSE, error=FALSE)
```

```{r data-import-prep, include = FALSE, message = FALSE}
#set column specifications for reading in data
spec <- cols(
  .default = col_character(),
  UNIT_ID = col_double(),
  OP_DATE = col_datetime(format = ""),
  OP_HOUR = col_double(),
  OP_TIME = col_double(),
  HEAT_INPUT = col_double(),
  GLOAD = col_double(),
  NOX_MASS = col_double(),
  NOX_RATE = col_double(),
  OP_YEAR = col_double(),
  ORIS_CODE = col_double(),
  EDAT_LOAD_VALUE = col_double(),
  HI_MEASURED_FLAG = col_double(),
  HI_SUB_FLAG = col_double(),
  HI_M_S_FLAG = col_double(),
  HI_LME = col_double(),
  HI_CALC = col_double(),
  HI_OTHER = col_double(),
  NOXM_MEASURED_FLAG = col_double(),
  NOXM_SUB_FLAG = col_double(),
  NOXM_M_S_FLAG = col_double(),
  NOXM_LME = col_double(),
  NOXM_CALC = col_double(),
  NOXR_MEASURED_FLAG = col_double(),
  NOXR_SUB_FLAG = col_double(),
  NOXR_M_S_FLAG = col_double()
)

#read in data for 2009-2019
data_2009 <- read_csv("data\\2009_data.csv", col_types = spec)
data_2010 <- read_csv("data\\2010_data.csv", col_types = spec)
data_2011 <- read_csv("data\\2011_data.csv", col_types = spec)
data_2012 <- read_csv("data\\2012_data.csv", col_types = spec)
data_2013 <- read_csv("data\\2013_data.csv", col_types = spec)
data_2014 <- read_csv("data\\2014_data.csv", col_types = spec)
data_2015 <- read_csv("data\\2015_data.csv", col_types = spec)
data_2016 <- read_csv("data\\2016_data.csv", col_types = spec)
data_2017 <- read_csv("data\\2017_data.csv", col_types = spec)
data_2018 <- read_csv("data\\2018_data.csv", col_types = spec)
data_2019 <- read_csv("data\\2019_data.csv", col_types = spec)

#combine all data years into one data frame
data_raw <-
  rbind(
    data_2009,
    data_2010,
    data_2011,
    data_2012,
    data_2013,
    data_2014,
    data_2015,
    data_2016,
    data_2017,
    data_2018,
    data_2019
  )

#data edits
  #filter out partial hours (keep only hours where the unit operates for the full hour)
  #Chalk Point Unit 2 has "SACR" which in this analysis will be classified as SNCR
  #Colver, Northampton, and Panther Creek report "Ammonia Injection" (NH3) as a control, which will be classified as SNCR    in this analysis
  #Indian River Unit 4 indicates both SCR and SNCR in the data, it will be grouped into SCR in this analysis
data <- data_raw %>%
  filter(OP_TIME == 1) %>%
  mutate(SNCR_FLAG = if_else(FACILITY_NAME == "Chalk Point" & UNITID == 2, "Y", SNCR_FLAG)) %>%
  mutate(SNCR_FLAG = if_else(FACILITY_NAME == "Colver Green Energy", "Y", SNCR_FLAG)) %>%
  mutate(SNCR_FLAG = if_else(FACILITY_NAME == "Northampton Generating Plant", "Y", SNCR_FLAG)) %>%
  mutate(SNCR_FLAG = if_else(FACILITY_NAME == "Panther Creek Energy Facility", "Y", SNCR_FLAG))
data$SNCR_FLAG[data$FACILITY_NAME == "Indian River" & data$UNITID == 4] <- NA

#clean up
rm(data_2009)
rm(data_2010)
rm(data_2011)
rm(data_2012)
rm(data_2013)
rm(data_2014)
rm(data_2015)
rm(data_2016)
rm(data_2017)
rm(data_2018)
rm(data_2019)
rm(data_raw)

#define variable
rate_line <- 0.12 #PA RACT II level (although it is 30-day rolling average value) for plots

```

```{r remove-substitute-data, include = FALSE, message = FALSE}
#count substitute data
noxr_sub <- sum(data$NOXR_SUB_FLAG)
hi_sub <- sum(data$HI_SUB_FLAG & data$NOXR_SUB_FLAG == 0)
noxm_sub <- 
  sum(data$NOXM_SUB_FLAG & data$HI_SUB_FLAG == 0 & data$NOXR_SUB_FLAG == 0)

#remove substitute data
data_temp <- data %>%
  filter(NOXR_SUB_FLAG != 1) %>%
  filter(HI_SUB_FLAG != 1) %>%
  filter(NOXM_SUB_FLAG != 1)

data_removed <- noxr_sub + hi_sub + noxm_sub
perc_removed <-
  round((data_removed / sum(data$NOXR_SUB_FLAG == 1 | data$NOXR_SUB_FLAG == 0) * 100), 2)
data_removed_text <-
  format(round(noxr_sub + hi_sub + noxm_sub), big.mark = ",")

#remove substitute data indicator columns
data <- data_temp %>%
  select(
    -HEAT_INPUT_MEASURE_FLG,
    -NOX_RATE_MEASURE_FLG,
    -NOX_MASS_MEASURE_FLG,
    -HI_MEASURED_FLAG:-NOXR_M_S_FLAG
  )

#clean up
rm(data_temp)

```

# Data Preparation

Partial hours of operation are excluded from this analysis. Hours in which substitute data were submitted for heat input, NO~X~ emissions rate, and/or NO~X~ mass emissions have also been removed from this analysis.^[Substitute data are required to be submitted in instances in which the monitor is not available or is not performing properly, or a quality assurance test of the monitoring equipment was missed. These substitute data values are usually similar to actual emissions; however in some rare cases, they may be significantly higher than actual emissions, sometimes representing the EGU’s maximum potential emissions.] In this analysis, ``r data_removed_text`` substitute data observations were removed, accounting for just ``r perc_removed``% of the data.

# NO~X~ Emission Rates and NO~X~ Mass Emissions by Operating Level for Individual EGUs 

This section of the analysis focuses on NO~X~ emission rates in the context of operating level, which can have an important impact on SCR performance and therefore NO~X~ emission rates and mass emissions. The patterns of operation of some coal-fired EGUs have been changing in recent years, and considering operating level allows any effects of these changes on NO~X~ emission rates and the operation of SCR controls to be taken into account. At lower levels of generation, the SCR may not be effective because temperatures may be too low to promote the catalytic reaction. Generally, the catalyst works more effectively as temperature in the unit increases with increasing heat input and correspondingly increasing steam and electrical output.

The plots in this section compare ozone season NO~X~ emission rates, NO~X~ mass emissions, and operating hours for individual coal-fired EGUs for either two or three years during the 2009-2019 period. Data for the year 2019 are always included for a recent snapshot. The other years shown are the years during the 2009-2019 period in which the unit achieved its lowest and third lowest average ozone season NO~X~ emission rates.^[Data for the ozone season with each unit’s third-lowest emission rate during the 2009-2019 period are included for comparability with the data sets EPA has used to identify emission reduction opportunities in the CSAPR Update and the proposed Revised CSAPR Update. Data for the ozone season with each unit’s lowest emission rate during the same period are also included for greater comparability with the data provided by the OTC to support its CAA section 184(c) recommendation. EPA notes that for some units the OTC has provided data for ozone seasons before 2009.] If plots are provided for only two years, then 2019 represents the year of either the unit’s lowest or third-lowest average emission rate during the 2009-2019 period.

This analysis uses heat input as a proxy for operating level. Specifically, this analysis uses a heat input factor. Heat input factor was calculated for each hour by dividing the heat input in a given hour by the individual EGU’s observed maximum hourly heat input during the ozone seasons in the 2009-2019 period. The hourly emission rate data are then divided into 10 bins according to the heat input factor for the hour (bin 1 contains data for all hours in which the heat input factor for the hour is between 0-10% of the unit’s maximum hourly heat input, bin 2 contains data for all hours in which the heat input factor is between 10-20% of the unit’s maximum hourly heat input, etc.).

The hourly emission rate data values for each heat input factor bin are displayed by means of a “box plot” or “box-and-whisker plot.” Each “box” represents the middle half of all the hourly data values in that heat input factor bin – that is, the hourly data values that fall in the “interquartile range” between the 25th percentile and 75th percentile hourly data values. The horizontal line in the box represents the median hourly data value. Vertical lines, or “whiskers,” extend to the highest and lowest hourly data values that fall above or below the top or bottom edges of the box within a distance of up to 1.5 times the interquartile range. Any outlying hourly data values that fall above or below the top or bottom edges of the box by a distance of more than 1.5 times the interquartile range are shown as individual dots. Thus, a lower median data value and lower overall placement of the box on the chart indicate generally lower hourly emission rates, while shorter vertical distances between the top and bottom edges of the box and between the top and bottom ends of the whiskers, as well as fewer outliers, indicate lower variability (or greater consistency) of a unit’s hourly emission rates at a given heat input factor bin. In this way, each box plot provides visual representations of both the magnitude and variability of a unit’s hourly NO~X~ emission rates at a given heat input factor bin in a single chart.

In all of the emission rate plots in this section, a horizontal dashed line showing the 0.12 lb/mmBtu emission rate limit that Pennsylvania’s SCR-equipped units are required to meet--under certain operating conditions, on a 30-day rolling average basis--by the state’s current RACT rules is included for comparison purposes.


```{r HI-bin-data-prep, include = FALSE, message = FALSE}
#calculate maximum heat input for each unit between 2010-2019 and join to main data
data_max_HI <- data %>%
  group_by(UNIT_ID) %>%
  summarise(HI_MAX = max(HEAT_INPUT))
data_HI_join <- full_join(data, data_max_HI)

#calculate heat input factor for each hour and sort into 10 even bins
data_binned <- data_HI_join %>%
  mutate(HI_FACTOR = HEAT_INPUT / HI_MAX) %>%
  mutate(HI_BIN = if_else(
    HI_FACTOR <= 0.1,
    "01",
    if_else(
      HI_FACTOR > 0.1 & HI_FACTOR <= 0.2,
      "02",
      if_else(
        HI_FACTOR > 0.2 & HI_FACTOR <= 0.3,
        "03",
        if_else(
          HI_FACTOR > 0.3 & HI_FACTOR <= 0.4,
          "04",
          if_else(
            HI_FACTOR > 0.4 & HI_FACTOR <= 0.5,
            "05",
            if_else(
              HI_FACTOR > 0.5 & HI_FACTOR <= 0.6,
              "06",
              if_else(
                HI_FACTOR > 0.6 & HI_FACTOR <= 0.7,
                "07",
                if_else(
                  HI_FACTOR > 0.7 & HI_FACTOR <= 0.8,
                  "08",
                  if_else(HI_FACTOR > 0.8 &
                            HI_FACTOR <= 0.9, "09", "10")
                )
              )
            )
          )
        )
      )
    )
  ))

#clean up
rm(data_max_HI)
rm(data_HI_join)

```

```{r min-min3-nox-rate-calc, include = FALSE, message = FALSE}
#data prep
data_binned_scr <- data_binned %>%
  filter(SCR_FLAG == "Y")

data_binned_sncr <- data_binned %>%
  filter(SNCR_FLAG == "Y")

data_binned_none <- data_binned %>%
  filter(is.na(SCR_FLAG) & is.na(SNCR_FLAG)) %>%
  filter(FACILITY_NAME != "Indian River") %>%
  filter(FACILITY_NAME != "Conemaugh") #remove facilities that have SCR later in operation

#function to calculate lowest and 3rd lowest ozone season rate and year for each unit
min_min3_calc <- function(df) {
  #calculate ozone season average NOX rate and find lowest and 3rd lowest rate
  data_noxr_os_avg <- df %>%
    group_by(UNIT_ID, FACILITY_NAME, UNITID, OP_YEAR) %>%
    summarise(NOX_RATE_OS_AVG = mean(NOX_RATE, na.rm = TRUE))
  
  data_os_min <- data_noxr_os_avg %>%
    group_by(UNIT_ID) %>%
    summarise(NOX_RATE_OS_MIN = min(NOX_RATE_OS_AVG))
  
  data_os_min3 <- data_noxr_os_avg %>%
    group_by(UNIT_ID) %>%
    arrange(NOX_RATE_OS_AVG) %>%
    summarise(NOX_RATE_OS_MIN3 = nth(NOX_RATE_OS_AVG, 3))
  
  #join lowest and 3rd lowest ozone season rate and year to main data for unit-level plots
  data_os_min_join <-
    left_join(
      data_os_min,
      data_noxr_os_avg,
      by = c("UNIT_ID" = "UNIT_ID", "NOX_RATE_OS_MIN" = "NOX_RATE_OS_AVG")
    ) %>%
    mutate(OP_YEAR_MIN = OP_YEAR) %>%
    select(-OP_YEAR)
  
  data_os_min3_join <-
    left_join(
      data_os_min3,
      data_noxr_os_avg,
      by = c("UNIT_ID" = "UNIT_ID", "NOX_RATE_OS_MIN3" = "NOX_RATE_OS_AVG")
    ) %>%
    mutate(OP_YEAR_MIN3 = OP_YEAR) %>%
    select(-OP_YEAR)
  
  df_mid <-
    left_join(
      df,
      data_os_min_join,
      by = c(
        "UNIT_ID" = "UNIT_ID",
        "FACILITY_NAME" = "FACILITY_NAME",
        "UNITID" = "UNITID"
      )
    )
  
  df_end <-
    left_join(
      df_mid,
      data_os_min3_join,
      by = c(
        "UNIT_ID" = "UNIT_ID",
        "FACILITY_NAME" = "FACILITY_NAME",
        "UNITID" = "UNITID"
      )
    )
  
  #sort data frame to print unit plots in order by state, facility name, and unit
  df_end <- df_end %>%
    arrange(STATE, FACILITY_NAME, UNITID)
  
}

#generate datasets
data_binned_scr <- min_min3_calc(data_binned_scr)
data_binned_sncr <- min_min3_calc(data_binned_sncr)
data_binned_none <- min_min3_calc(data_binned_none)

```

```{r function-unit-op-level, include = FALSE, message = FALSE}
#creates plots comparing NOX rates to unit total NOX emissions and operating hours in each heat input factor bin
op_level_nox_unit_plot <- function(df1, title) {
  #filter df1 for years to be plotted
  df2 <- df1 %>%
    filter(OP_YEAR == OP_YEAR_MIN |
             OP_YEAR == OP_YEAR_MIN3 |
             OP_YEAR == 2019)
  
  #generate column with labels for plots
  df2 <- df2 %>%
    mutate(OP_YEAR_label = if_else(
      OP_YEAR == OP_YEAR_MIN,
      paste(OP_YEAR_MIN, "(lowest)"),
      if_else(
        OP_YEAR == OP_YEAR_MIN3,
        paste(OP_YEAR_MIN3, "(3rd lowest)"),
        "2019"
      )
    ))
  
  #set y-axis maximum for plot p3a/p3b
  ylim <- max(df2$NOX_RATE)
  df3 <- df2 %>%
    group_by(OP_YEAR, OP_YEAR_label, HI_BIN) %>%
    summarise(
      OP_TIME = sum(OP_TIME, na.rm = TRUE),
      NOX_MASS = sum(NOX_MASS, na.rm = TRUE)
    ) %>%
    mutate(NOX_MASS = NOX_MASS / 2000)
  
  #generate operating hours bar chart
  p1 <-
    ggplot(df3, aes(x = HI_BIN, y = OP_TIME, fill = OP_YEAR_label)) +
    geom_col(position = position_dodge2(preserve = "single")) +
    theme_bw() +
    labs(
      x = 'Heat Input Factor Bin',
      y = 'Operating Hours',
      fill = 'Year',
      title = title
    ) +
    scale_y_continuous(labels = scales::comma_format(accuracy = 1)) +
    scale_fill_grey()
  p2 <- grid.arrange(p1, ncol = 1)
  
  #generate NOX emission rate box plot
  #p3a plots the rate data as-is, p3b plots the rate data with a limit on the y-axis scale because certain units were identified as having extreme outliers in years displayed on the plots, which interfered with the visualization
  #units needing this adjustment are specified when p5 is generated
  p3a <- ggplot(df2, aes(x = HI_BIN, y = NOX_RATE)) +
    geom_boxplot() +
    theme_bw() +
    labs(
      x = element_blank(),
      y = expression('NO'['X'] * ' Rate (lb/mmBtu)'),
      title = title
    ) +
    facet_grid(. ~ OP_YEAR_label) +
    geom_hline(yintercept = rate_line,
               linetype = "dashed",
               color = "black") +
    scale_y_continuous(
      labels = scales::number_format(accuracy = 0.1),
      breaks = seq(0, ylim, by = 0.1)
    )
  p3b <- ggplot(df2, aes(x = HI_BIN, y = NOX_RATE)) +
    geom_boxplot() +
    theme_bw() +
    labs(
      x = element_blank(),
      y = expression('NO'['X'] * ' Rate (lb/mmBtu)'),
      title = title
    ) +
    facet_grid(. ~ OP_YEAR_label) +
    geom_hline(yintercept = rate_line,
               linetype = "dashed",
               color = "black") +
    scale_y_continuous(
      labels = scales::number_format(accuracy = 0.1),
      breaks = seq(0, ylim, by = 0.1)
    ) +
    coord_cartesian(ylim = c(0, 0.45))
  
  #generate NOX mass emissions bar chart
  p4 <- ggplot(df3, aes(x = HI_BIN, y = NOX_MASS)) +
    geom_col() +
    theme_bw() +
    labs(x = 'Heat Input Factor Bin', y = expression('NO'['X'] * ' Emissions (tons)')) +
    facet_grid(. ~ OP_YEAR_label) +
    scale_y_continuous(labels = scales::comma_format(accuracy = 1))
  
  #generate NOX emission rate and mass emissions joint plot
  #if-else statement used to generate p5 because three units needed limit on y-axis scale for data viz purposes
  p5 <- grid.arrange(if (df2$UNIT_ID == 1066 |
                         df2$UNIT_ID == 3555 |
                         df2$UNIT_ID == 3529)
    p3b
    else
      p3a, p4, ncol = 1)
}

```

## Units with SCR Controls

```{r scr-unit-op-level-plots, echo = FALSE, message = FALSE, fig.width = 7, fig.height = 4.5}
#loop to create unit level plots for all units in dataset
for (unit in unique(data_binned_scr$UNIT_ID)) {
  df <- subset(data_binned_scr, UNIT_ID == unit)
  op_level_nox_unit_plot(df,
                         paste(df$FACILITY_NAME,
                               "Unit",
                               df$UNITID,
                               df$STATE))
}

```

## Units with SNCR Controls

```{r sncr-unit-op-level-plots, echo = FALSE, message = FALSE, fig.width = 7, fig.height = 4.5}
#loop to create unit level plots for all units in dataset
for (unit in unique(data_binned_sncr$UNIT_ID)) {
  df <- subset(data_binned_sncr, UNIT_ID == unit)
  op_level_nox_unit_plot(df,
                         paste(df$FACILITY_NAME,
                               "Unit",
                               df$UNITID,
                               df$STATE))
}

```

## Units Without SCR or SNCR Controls

```{r none-unit-op-level-plots, echo = FALSE, message = FALSE, fig.width = 7, fig.height = 4.5}
#loop to create unit level plots for all units in dataset
for (unit in unique(data_binned_none$UNIT_ID)) {
  df <- subset(data_binned_none, UNIT_ID == unit)
  op_level_nox_unit_plot(df,
                         paste(df$FACILITY_NAME,
                               "Unit",
                               df$UNITID,
                               df$STATE))
}

```


