#########################################################################
#												                                                #
#				                EJ National Analysis                            #
#												                                                #
#########################################################################

# Author: Alyssa Le
# Date Created: 09/17/21

rm(list=ls())
options(stringsAsFactors=FALSE)

library(dplyr)
library(readxl)
library(sf)

ORM <- "C:/Users/50526/ICF/WOTUS Reconsideration (ICF Internal Site) - ORM Data Processing/Outputs/"
EJ <- "C:/Users/50526/ICF/WOTUS Reconsideration (ICF Internal Site) - GIS - Environmental Justice/"
WQ <- "C:/Users/50526/ICF/WOTUS Reconsideration (ICF Internal Site) - Run Results/"
GIS <- "C:/Users/50526/OneDrive - ICF/Everest/GIS_Data/"

wet_area_break <- c(0.5, 1, 50)
water_break <- c(1, 5, 50)

race_var <- c("BLACK", "NATIVE", "ASIAN", "HAWAII", "SOMEOTHER", "TWORACES", "HISPANIC")

env_var <- c("DSLPM", "CANCER", "RESP", "PTRAF", "PWDIS", "PNPL", "PRMP", "PTSDF", "OZONE", "PM25")

watersheds <- c("0514_Analysis_Baseline-Policy_2022-08-22.xlsx", 
                "1027_Analysis_Baseline-Policy_2022-07-28.xlsx")

# PART 1: Read in ORM data and preprocess it. ---------------------------------

## ORM data
ORM_mitigation_temp <- read.csv(paste0(ORM,"Impacts_By_HUC12_2022.07.28.csv"),
                           colClasses=c("character",rep("numeric",18)))
ORM_mitigation_temp$TOTAL_CHANGE <- with(ORM_mitigation_temp, (PERM_ACRES_CHANGE*-1)*20 + 
                                           (TEMP_ACRES_CHANGE*-1)*15)

ORM_mitigation <- ORM_mitigation_temp %>% select(HUC12, TOTAL_CHANGE)

# write.csv(ORM_mitigation,paste0(EJ,"R_programs/Outputs/ORM_mitigation_GIS.csv"),row.names=F)

ORM_impactedwaters_temp <- read.csv(paste0(ORM,"ImpactedWaters_By_HUC12_2022.07.28.csv"),
                               colClasses=c("character",rep("numeric",9)))
ORM_impactedwaters_temp$WATERS_CHANGE <- with(ORM_impactedwaters_temp, Waters_Count_PERM_Change*20 + 
                                          Waters_Count_TEMP_Change*15)

ORM_impactedwaters <- ORM_impactedwaters_temp %>% select(HUC12, WATERS_CHANGE)

# write.csv(ORM_impactedwaters,paste0(EJ,"R_programs/Outputs/ORM_IW_GIS.csv"),row.names=F)

# PART 2: Read in socioeconomic data and preprocess ---------------------------

## ACS data
# Block group level
BG_HUC12_Intersect <- read.csv(paste0(EJ,"GIS_Outputs/AllStatesBG_HUC12_20211006.csv"),
                               colClasses = c("numeric","NULL","character",rep("integer",44),
                                              "character",rep("numeric",2)))
sort(unique(substr(BG_HUC12_Intersect$huc12,1,2)))

# HUC12s may be missing the leading zero; include this code as a safeguard
BG_HUC12_Intersect$huc12 <- with(BG_HUC12_Intersect, ifelse(nchar(huc12) < 12, paste0("0", huc12), huc12))
sort(unique(substr(BG_HUC12_Intersect$huc12,1,2)))
# QA (11/1) - confirmed have 19 HUC2s included in the dataset

# QA: Number of Unique HUC12s
length(unique(BG_HUC12_Intersect$huc12)) # 97225

# GEOIDs are also missing the leading zero
BG_HUC12_Intersect$GEOID <- with(BG_HUC12_Intersect, ifelse(nchar(GEOID) < 12, paste0("0", GEOID), GEOID))

# Consolidate variables
BG_HUC12_Intersect$UNDER_5 <- BG_HUC12_Intersect %>%
  select(B01001_003, B01001_027) %>%
  rowSums(na.rm=TRUE)
BG_HUC12_Intersect$OVER_64 <- BG_HUC12_Intersect %>%
  select(B01001_020:B01001_025,B01001_044:B01001_049) %>%
  rowSums(na.rm=TRUE)
BG_HUC12_Intersect$MINOR_POP <- BG_HUC12_Intersect %>%
  select(B03002_004:B03002_012) %>%
  rowSums(na.rm=TRUE)
BG_HUC12_Intersect$POV_POP <- BG_HUC12_Intersect %>%
  select(C17002_002:C17002_007) %>%
  rowSums(na.rm=TRUE)

BG_HUC12_Var <- BG_HUC12_Intersect %>%
  select(HUC12=huc12, GEOID, TOT_AGE_POP=B01001_001, UNDER_5=B01001_003, OVER_64,
         TOT_RACE_POP=B03002_001, MINOR_POP, BLACK=B03002_004, NATIVE=B03002_005,
         ASIAN=B03002_006, HAWAII=B03002_007, SOMEOTHER=B03002_008, TWORACES = B03002_009,
         HISPANIC=B03002_012, TOT_POV_POP=C17002_001, POV_POP, TOT_EDU_POP=B15002_001, 
         EDU_POP=estimate_edu, TOT_LINGISO_HH=C16002_001, LINGISO_HH=estimate_lingiso, 
         HUC12_SqKM, BG_SqKM=BG_Full_KM, BG_HUC12_SqKm)

## EJScreen data
# (11/1) Updated the variables to be read in based on EPA technical direction
# Had to read in the variables as character because got a "scan() expected 'a real', got 'None'" error
EJScreen_data <- read.csv(paste0(EJ,"R_programs/Inputs/EJSCREEN_2019_USPR.csv"), header=TRUE,
                          colClasses = c("NULL","character",rep("NULL",4),"numeric",rep("NULL",12),"numeric",
                          rep("NULL",7),rep("character",10),rep("NULL",330)))

for (i in 1:length(env_var)) {

  EJScreen_data[, paste0(env_var[i], "_num")] <- as.numeric(EJScreen_data[, env_var[i]])
    
}
# This introduces "NA's" where this variable was set to 'None' in the EJScreen dataset.
#*# KM (11/2): Confirmed correct selections of EJScreen data

BG_HUC12_EJ <- merge(BG_HUC12_Var, EJScreen_data %>% select(ID, ACSTOTHU, PRE1960, DSLPM_num:PM25_num), 
                     by.x="GEOID", by.y="ID", all=TRUE)

length(unique(BG_HUC12_EJ$HUC12)) #*# SE: 97226, new HUC12 added?
#*# AL: There isn't a perfect match between these two datasets. This is okay.

# rm(EJScreen_data, BG_HUC12_Intersect)

# QA
BG_HUC12_EJ_QA <- BG_HUC12_EJ[is.na(BG_HUC12_EJ$TOT_AGE_POP), ]
# ~3,500 records from EJScreen missing from ACS summary.
#*# AL (10/6): Includes GEOIDs that fall outside of the scope of this analysis (e.g., in Puerto Rico)

# Adjust variables by CBG/HUC12 overlap
BG_HUC12_EJ$BG_RATIO <- with(BG_HUC12_EJ,BG_HUC12_SqKm/BG_SqKM)
# QA: Confirmed that BG_RATIO values range from 0 to 1
min(BG_HUC12_EJ$BG_RATIO, na.rm=T) # 0
max(BG_HUC12_EJ$BG_RATIO, na.rm=T) # 1.000002 #*# SE: confirmed values

QA_BG_HUC12_EJ <- BG_HUC12_EJ %>% select(c("GEOID","BG_RATIO")) %>% 
  group_by(GEOID) %>% 
  summarise_all(sum) #*# SE: 5108 GEOIDs have a ratio > 1

#*# AL (10/6): The majority of BG_Ratios sum to 1 or nearly 1. 

length(QA_BG_HUC12_EJ$GEOID[QA_BG_HUC12_EJ$BG_RATIO < 0.99 & !is.na(QA_BG_HUC12_EJ$BG_RATIO)])/
  length(QA_BG_HUC12_EJ$GEOID[!is.na(QA_BG_HUC12_EJ$BG_RATIO)])

BG_HUC12_adj <- BG_HUC12_EJ %>%
  mutate_at(c("TOT_AGE_POP","UNDER_5","OVER_64", "TOT_RACE_POP", "MINOR_POP", "BLACK", 
              "NATIVE", "ASIAN", "HAWAII", "SOMEOTHER", "TWORACES", "HISPANIC", "TOT_POV_POP",
              "POV_POP", "TOT_EDU_POP", "EDU_POP", "TOT_LINGISO_HH", "LINGISO_HH", "ACSTOTHU", "PRE1960"), 
            ~ (.x * BG_RATIO)) %>%
  mutate_at(c("DSLPM_num", "CANCER_num", "RESP_num", "PTRAF_num", "PWDIS_num", "PNPL_num", "PRMP_num", "OZONE_num",
              "PM25_num"), ~ (.x * BG_HUC12_SqKm))

BG_HUC12_temp <- BG_HUC12_adj %>%
  select(HUC12, BG_HUC12_SqKm, TOT_AGE_POP:LINGISO_HH, ACSTOTHU:PM25_num) %>%
  group_by(HUC12) %>%
  summarise(across(TOT_AGE_POP:PRE1960, ~ sum(.x, na.rm=T)),
            across(DSLPM_num:PM25_num, ~ sum(.x, na.rm=T)/sum(BG_HUC12_SqKm)))

# Identify state with the largest overlap with the HUC12 boundary.
  # QA note: this output does not include Hawaii. EPA confirmed we can drop HI from the EJ analysis though.
BG_HUC12_adj$STATE_ID <- substr(BG_HUC12_adj$GEOID, 1, 2)

BG_HUC12_adj_sort <- BG_HUC12_adj[with(BG_HUC12_adj, order(HUC12, -BG_HUC12_SqKm)), ]

BG_HUC12_uniquestate <- BG_HUC12_adj_sort %>%
  select(HUC12, STATE_ID) %>%
  distinct(HUC12, .keep_all=T)

BG_HUC12_Final <- merge(BG_HUC12_temp, BG_HUC12_uniquestate, by="HUC12", all.x=T)

length(unique(BG_HUC12_Final$HUC12)) #*# SE: 97226 HUC12s

# PART 3: Join socioeconomic data with wetland or impacted waters changes -----

## Wetland Area Changes
ORM_mitig_ACS <- merge(BG_HUC12_Final, ORM_mitigation, by="HUC12", all=T)

length(unique(ORM_mitig_ACS$HUC12)) # 97282 HUC12s?

# rm(ORM_mitigation)

# QA
#*# SE: both confirmed to have 0 records
ORM_mitig_ACS_QA <- ORM_mitig_ACS[which(is.na(ORM_mitig_ACS$PERC_MIN)&
                                          ORM_mitig_ACS$TOTAL_CHANGE!=0&
                                          substr(ORM_mitig_ACS$HUC12,1,2)!="20"), ]
ORM_mitig_ACS_QATri <- ORM_mitig_ACS[which(is.na(ORM_mitig_ACS$PERC_TRIBAL)&
                                             ORM_mitig_ACS$TOTAL_CHANGE!=0&
                                             substr(ORM_mitig_ACS$HUC12,1,2)!="20"), ]
ORM_mitig_ACS_dup <- ORM_mitigation[duplicated(ORM_mitigation$HUC12), ]

## Impacted Waters Changes
ORM_IW_ACS <-  merge(BG_HUC12_Final, ORM_impactedwaters, by="HUC12", all=T)

length(unique(ORM_IW_ACS$HUC12)) # 97282 HUC12s?

# rm(ORM_impactedwaters)

# QA
ORM_IW_ACS_QA <- ORM_IW_ACS[which(is.na(ORM_IW_ACS$PERC_MIN)&ORM_IW_ACS$WATERS_CHANGE!=0&
                                    substr(ORM_IW_ACS$HUC12,1,2)!="20"), ] #*# SE: 0 observations

Natl_ORM_mitig_QA <- ORM_mitig_ACS %>%
  dplyr::filter(is.na(TOTAL_CHANGE))

length(unique(Natl_ORM_mitig_QA$HUC12)) #*# SE: 54035 HUC12s with no changes

Natl_ORM_IW_QA <- ORM_IW_ACS %>%
  dplyr::filter(is.na(WATERS_CHANGE)) #*# SE: 54035 HUC12s with no changes

# Since the two datasets cover the same universe, use just one to determine the population
#    "outside" of the affected watersheds.
identical(Natl_ORM_mitig_QA$HUC12, Natl_ORM_IW_QA$HUC12) # TRUE
#*# SE: confirmed True

# PART 4: Calculate HUC12 national and state statistics -----------------------

# Calculate national statistics
Natl_stats_tmp <- ORM_mitig_ACS %>%
  dplyr::filter(is.na(TOTAL_CHANGE)) %>%
  select(TOT_AGE_POP:PRE1960) %>%
  summarise_all(sum, na.rm=T)

Natl_stats_tmp$PERC_MIN <- with(Natl_stats_tmp, MINOR_POP/TOT_RACE_POP)
Natl_stats_tmp$PERC_POV <- with(Natl_stats_tmp, POV_POP/TOT_POV_POP)
Natl_stats_tmp$PERC_LESSHS <- with(Natl_stats_tmp, EDU_POP/TOT_EDU_POP)
Natl_stats_tmp$PERC_LING <- with(Natl_stats_tmp, LINGISO_HH/TOT_LINGISO_HH)
Natl_stats_tmp$PERC_UNDER5 <- with(Natl_stats_tmp, UNDER_5/TOT_AGE_POP)
Natl_stats_tmp$PERC_OVER64 <- with(Natl_stats_tmp, OVER_64/TOT_AGE_POP)
Natl_stats_tmp$PERC_LEAD <- with(Natl_stats_tmp, PRE1960/ACSTOTHU)

for (i in 1:length(race_var)) {
  
  Natl_stats_tmp[, paste0("PERC_",race_var[i])] <- 
    Natl_stats_tmp[, race_var[i]] / Natl_stats_tmp$TOT_RACE_POP
    
}

Natl_mean <- ORM_mitig_ACS %>%
  dplyr::filter(is.na(TOTAL_CHANGE)) %>%  
  select(DSLPM_num:PM25_num) %>%
  summarise_all(mean, na.rm=T)

Natl_stats <- cbind(Natl_stats_tmp, Natl_mean)

# write.csv(Natl_stats %>% select(PERC_MIN:PM25_num), paste0(EJ,"R_programs/Outputs/Natl_Stats.csv"),
#           row.names=F)

Natl_stats_tmp_env <- ORM_mitig_ACS %>%
  dplyr::filter(is.na(TOTAL_CHANGE) & 
                  !is.na(PRMP_num) & PRMP_num>0 & 
                  !is.na(PWDIS_num) & PWDIS_num>0) %>%
  select(TOT_AGE_POP:PRE1960) %>%
  summarise_all(sum, na.rm=T)

Natl_stats_tmp_env$PERC_MIN <- with(Natl_stats_tmp_env, MINOR_POP/TOT_RACE_POP)
Natl_stats_tmp_env$PERC_POV <- with(Natl_stats_tmp_env, POV_POP/TOT_POV_POP)
Natl_stats_tmp_env$PERC_LESSHS <- with(Natl_stats_tmp_env, EDU_POP/TOT_EDU_POP)
Natl_stats_tmp_env$PERC_LING <- with(Natl_stats_tmp_env, LINGISO_HH/TOT_LINGISO_HH)
Natl_stats_tmp_env$PERC_UNDER5 <- with(Natl_stats_tmp_env, UNDER_5/TOT_AGE_POP)
Natl_stats_tmp_env$PERC_OVER64 <- with(Natl_stats_tmp_env, OVER_64/TOT_AGE_POP)
Natl_stats_tmp_env$PERC_LEAD <- with(Natl_stats_tmp_env, PRE1960/ACSTOTHU)

for (i in 1:length(race_var)) {
  
  Natl_stats_tmp_env[, paste0("PERC_",race_var[i])] <- 
    Natl_stats_tmp_env[, race_var[i]] / Natl_stats_tmp_env$TOT_RACE_POP
  
}

Natl_mean_env <- ORM_mitig_ACS %>%
  dplyr::filter(is.na(TOTAL_CHANGE) & 
                  !is.na(PRMP_num) & PRMP_num>0 & 
                  !is.na(PWDIS_num) & PWDIS_num>0) %>%  
  select(DSLPM_num:PM25_num) %>%
  summarise_all(mean, na.rm=T)

Natl_stats_env <- cbind(Natl_stats_tmp_env, Natl_mean_env)

# write.csv(Natl_stats_env %>% select(PERC_MIN:PM25_num), paste0(EJ,"R_programs/Outputs/Natl_Stats_EJ.csv"),
#           row.names=F)

# Calculate state statistics
State_stats_tmp <- ORM_mitig_ACS %>%
  dplyr::filter(is.na(TOTAL_CHANGE)) %>%  
  select(STATE_ID, TOT_AGE_POP:PRE1960) %>%
  group_by(STATE_ID) %>%
  summarise_all(sum, na.rm=T)

State_stats_tmp$PERC_MIN <- with(State_stats_tmp, MINOR_POP/TOT_RACE_POP)
State_stats_tmp$PERC_POV <- with(State_stats_tmp, POV_POP/TOT_POV_POP)
State_stats_tmp$PERC_LESSHS <- with(State_stats_tmp, EDU_POP/TOT_EDU_POP)
State_stats_tmp$PERC_LING <- with(State_stats_tmp, LINGISO_HH/TOT_LINGISO_HH)
State_stats_tmp$PERC_UNDER5 <- with(State_stats_tmp, UNDER_5/TOT_AGE_POP)
State_stats_tmp$PERC_OVER64 <- with(State_stats_tmp, OVER_64/TOT_AGE_POP)
State_stats_tmp$PERC_LEAD <- with(State_stats_tmp, PRE1960/ACSTOTHU)

for (i in 1:length(race_var)) {
  
  State_stats_tmp[, paste0("PERC_",race_var[i])] <- 
    State_stats_tmp[, race_var[i]] / State_stats_tmp$TOT_RACE_POP
  
}

State_mean <- ORM_mitig_ACS %>%
  dplyr::filter(is.na(TOTAL_CHANGE)) %>%  
  select(STATE_ID, DSLPM_num:PM25_num) %>%
  group_by(STATE_ID) %>%
  summarise_all(mean, na.rm=T)

State_stats <- merge(State_stats_tmp, State_mean, by="STATE_ID", all.x=T)

# write.csv(State_stats %>% select(STATE_ID, PERC_MIN:PM25_num), paste0(EJ,"R_programs/Outputs/State_stats.csv"),
#           row.names=F)

State_stats_tmp_env <- ORM_mitig_ACS %>%
  dplyr::filter(is.na(TOTAL_CHANGE) & 
                !is.na(PRMP_num) & PRMP_num>0 & 
                  !is.na(PWDIS_num) & PWDIS_num>0) %>%  
  select(STATE_ID, TOT_AGE_POP:PRE1960) %>%
  group_by(STATE_ID) %>%
  summarise_all(sum, na.rm=T)

State_stats_tmp_env$PERC_MIN <- with(State_stats_tmp_env, MINOR_POP/TOT_RACE_POP)
State_stats_tmp_env$PERC_POV <- with(State_stats_tmp_env, POV_POP/TOT_POV_POP)
State_stats_tmp_env$PERC_LESSHS <- with(State_stats_tmp_env, EDU_POP/TOT_EDU_POP)
State_stats_tmp_env$PERC_LING <- with(State_stats_tmp_env, LINGISO_HH/TOT_LINGISO_HH)
State_stats_tmp_env$PERC_UNDER5 <- with(State_stats_tmp_env, UNDER_5/TOT_AGE_POP)
State_stats_tmp_env$PERC_OVER64 <- with(State_stats_tmp_env, OVER_64/TOT_AGE_POP)
State_stats_tmp_env$PERC_LEAD <- with(State_stats_tmp_env, PRE1960/ACSTOTHU)

for (i in 1:length(race_var)) {
  
  State_stats_tmp_env[, paste0("PERC_",race_var[i])] <- 
    State_stats_tmp_env[, race_var[i]] / State_stats_tmp_env$TOT_RACE_POP
  
}

#*# SE: The national version of the following statement does not include the
#*# "is.na(TOTAL_CHANGE)" in the filter statement
#*# AL: Corrected above.

State_mean_env <- ORM_mitig_ACS %>%
  dplyr::filter(is.na(TOTAL_CHANGE) & 
                !is.na(PRMP_num) & PRMP_num>0 & 
                  !is.na(PWDIS_num) & PWDIS_num>0) %>%  
  select(STATE_ID, DSLPM_num:PM25_num) %>%
  group_by(STATE_ID) %>%
  summarise_all(mean, na.rm=T)

State_stats_env <- merge(State_stats_tmp_env, State_mean_env, by="STATE_ID", all.x=T)

# write.csv(State_stats_env %>% select(STATE_ID, PERC_MIN:PM25_num), paste0(EJ,"R_programs/Outputs/State_stats_EJ.csv"),
#           row.names=F)
  
# PART 5: Bin the outputs based on wetland area/impacted waters changes -------

## Wetland Area Changes
Mitig_ACS_quant <- list()
Mitig_ACS_env <- list()
Mitig_ACS_totals <- list()
Mitig_ACS_env_states <- list()
Mitig_ACS_states <- list()
Mitig_ACS_states_fin <- list()

Mitig_ACS_EJ_quant <- list()
Mitig_ACS_EJ_totals <- list()
Mitig_ACS_EJ_states <- list()
Mitig_ACS_EJ_states_fin <- list()

Mitig_ACS_env_table <- data.frame(c(0:11))
names(Mitig_ACS_env_table) <- "Exceedances"

Mitig_ACS_env_states_table <- data.frame(c(0:11))
names(Mitig_ACS_env_states_table) <- "Exceedances"

for (i in 1:4) {
  
  # Splitting up the data into the bins agreed upon with EPA
  if (i==1) {
    Mitig_ACS_quant[[i]] <- ORM_mitig_ACS %>%
      dplyr::filter(!is.na(TOTAL_CHANGE) & TOTAL_CHANGE<=wet_area_break[i])
    
    Mitig_ACS_EJ_quant[[i]] <- ORM_mitig_ACS %>%
      dplyr::filter(!is.na(TOTAL_CHANGE) & TOTAL_CHANGE<=wet_area_break[i] & 
                      !is.na(PRMP_num) & PRMP_num>0 &
                      !is.na(PWDIS_num) & PWDIS_num>0)
    
  } else if (i<4) {
    Mitig_ACS_quant[[i]] <- ORM_mitig_ACS %>%
      dplyr::filter(!is.na(TOTAL_CHANGE) & wet_area_break[i-1]<TOTAL_CHANGE &
                      TOTAL_CHANGE<=wet_area_break[i])
    
    Mitig_ACS_EJ_quant[[i]] <- ORM_mitig_ACS %>%
      dplyr::filter(!is.na(TOTAL_CHANGE) & wet_area_break[i-1]<TOTAL_CHANGE & 
                      TOTAL_CHANGE<=wet_area_break[i] &
                      !is.na(PRMP_num) & PRMP_num>0 & 
                      !is.na(PWDIS_num) & PWDIS_num>0)
    
  } else if (i==4) {
    Mitig_ACS_quant[[i]] <- ORM_mitig_ACS %>%
      dplyr::filter(!is.na(TOTAL_CHANGE) & TOTAL_CHANGE>wet_area_break[i-1])
    
    Mitig_ACS_EJ_quant[[i]] <- ORM_mitig_ACS %>%
      dplyr::filter(!is.na(TOTAL_CHANGE) & TOTAL_CHANGE>wet_area_break[i-1] & 
                      !is.na(PRMP_num) & PRMP_num>0 & 
                      !is.na(PWDIS_num) & PWDIS_num>0)
    
  }
  
  for (j in 1:length(env_var)) {
    
    Mitig_ACS_quant[[i]][, paste0("EXCEED_",env_var[j])] <- 
      ifelse(Mitig_ACS_quant[[i]][, paste0(env_var[j],"_num")] > 
               Natl_stats[, paste0(env_var[j],"_num")], 1, 0)
    
    jpeg(paste0(EJ,"R_programs/Outputs/Wet_Natl_Exceedances_",env_var[j],"_",i,".jpg"))
    hist(Mitig_ACS_quant[[i]][, paste0(env_var[j],"_num")])
    dev.off()
    
  }
  
  Mitig_ACS_quant[[i]]$PERC_LEAD_ENV <- with(Mitig_ACS_quant[[i]], PRE1960/ACSTOTHU)
  
  jpeg(paste0(EJ,"R_programs/Outputs/Wet_Natl_Exceedances_PERC_LEAD_",i,".jpg"))
  hist(Mitig_ACS_quant[[i]]$PERC_LEAD_ENV)
  dev.off()
  
  Mitig_ACS_quant[[i]]$EXCEED_LEAD <- with(Mitig_ACS_quant[[i]], 
                                           ifelse(PERC_LEAD_ENV > 
                                                    Natl_stats$PERC_LEAD, 1, 0))
  
  exceed_colname <- paste0("EXCEED_CNT_",i)
  
  Mitig_ACS_env[[i]] <- Mitig_ACS_quant[[i]] %>%
    mutate(EXCEED_TOTAL = rowSums(across(starts_with("EXCEED")), na.rm=T)) %>%
    group_by(EXCEED_TOTAL) %>%
    summarize(!!exceed_colname := n_distinct(HUC12))
  
  Mitig_ACS_env_table <- merge(Mitig_ACS_env_table, Mitig_ACS_env[[i]], by.x="Exceedances",
                               by.y="EXCEED_TOTAL", all.x=T)
  
  Mitig_ACS_totals[[i]] <- data.frame(cbind(t(colSums(Mitig_ACS_quant[[i]]%>%select(TOT_AGE_POP:PRE1960), na.rm=T)),
                                             t(colMeans(Mitig_ACS_quant[[i]]%>%select(DSLPM_num:PM25_num), na.rm=T))))
  
  Mitig_ACS_totals[[i]]$PERC_MIN <- with(Mitig_ACS_totals[[i]], MINOR_POP/TOT_RACE_POP)
  Mitig_ACS_totals[[i]]$PERC_POV <- with(Mitig_ACS_totals[[i]], POV_POP/TOT_POV_POP)
  Mitig_ACS_totals[[i]]$PERC_LESSHS <- with(Mitig_ACS_totals[[i]], EDU_POP/TOT_EDU_POP)
  Mitig_ACS_totals[[i]]$PERC_LING <- with(Mitig_ACS_totals[[i]], LINGISO_HH/TOT_LINGISO_HH)
  Mitig_ACS_totals[[i]]$PERC_UNDER5 <- with(Mitig_ACS_totals[[i]], UNDER_5/TOT_AGE_POP)
  Mitig_ACS_totals[[i]]$PERC_OVER64 <- with(Mitig_ACS_totals[[i]], OVER_64/TOT_AGE_POP)
  Mitig_ACS_totals[[i]]$PERC_LEAD <- with(Mitig_ACS_totals[[i]], PRE1960/ACSTOTHU)
  
  for (j in 1:length(race_var)) {
    
    Mitig_ACS_totals[[i]][, paste0("PERC_",race_var[j])] <- 
      Mitig_ACS_totals[[i]][, race_var[j]] / Mitig_ACS_totals[[i]]$TOT_RACE_POP
    
  }
  
  Mitig_ACS_totals[[i]]$HUC12_Num <- length(unique(Mitig_ACS_quant[[i]]$HUC12))
  Mitig_ACS_totals[[i]]$ORM_Total <- length(unique(ORM_mitigation$HUC12))
  Mitig_ACS_totals[[i]]$HUC12_Perc <- with(Mitig_ACS_totals[[i]], HUC12_Num/ORM_Total)
  
  if(exists("Mitig_ACS_totals_table")) {
    Mitig_ACS_totals_table <- rbind(Mitig_ACS_totals_table, Mitig_ACS_totals[[i]])
  } else {
    Mitig_ACS_totals_table <- Mitig_ACS_totals[[i]]
  }
  
  Mitig_ACS_states[[i]] <- Mitig_ACS_quant[[i]]
  
  Mitig_ACS_states[[i]]$PERC_MIN_WET <- with(Mitig_ACS_states[[i]], MINOR_POP/TOT_RACE_POP)
  Mitig_ACS_states[[i]]$PERC_POV_WET <- with(Mitig_ACS_states[[i]], POV_POP/TOT_POV_POP)
  Mitig_ACS_states[[i]]$PERC_LESSHS_WET <- with(Mitig_ACS_states[[i]], EDU_POP/TOT_EDU_POP)
  Mitig_ACS_states[[i]]$PERC_LING_WET <- with(Mitig_ACS_states[[i]], LINGISO_HH/TOT_LINGISO_HH)
  Mitig_ACS_states[[i]]$PERC_UNDER5_WET <- with(Mitig_ACS_states[[i]], UNDER_5/TOT_AGE_POP)
  Mitig_ACS_states[[i]]$PERC_OVER64_WET <- with(Mitig_ACS_states[[i]], OVER_64/TOT_AGE_POP)
  Mitig_ACS_states[[i]]$PERC_LEAD_WET <- with(Mitig_ACS_states[[i]], PRE1960/ACSTOTHU)
  
  for (j in 1:length(race_var)) {
    
    Mitig_ACS_states[[i]][, paste0("PERC_",race_var[j],"_WET")] <- 
      Mitig_ACS_states[[i]][, race_var[j]] / Mitig_ACS_states[[i]]$TOT_RACE_POP
    
  }
  
  Mitig_ACS_states[[i]] <- merge(Mitig_ACS_states[[i]], 
                                 State_stats %>% select(STATE_ID, PERC_MIN:PM25_num), by="STATE_ID", all.x=T)
  
  Mitig_ACS_states[[i]]$EXCEED_MIN <- with(Mitig_ACS_states[[i]], ifelse(!is.na(PERC_MIN) & !is.na(PERC_MIN_WET) & 
                                                                           PERC_MIN_WET>PERC_MIN, 1, 0))
  Mitig_ACS_states[[i]]$EXCEED_POV <- with(Mitig_ACS_states[[i]], ifelse(!is.na(PERC_POV) & !is.na(PERC_POV_WET) &
                                                                           PERC_POV_WET>PERC_POV, 1, 0))
  Mitig_ACS_states[[i]]$EXCEED_LESSHS <- with(Mitig_ACS_states[[i]], ifelse(!is.na(PERC_LESSHS) & !is.na(PERC_LESSHS_WET) &
                                                                              PERC_LESSHS_WET>PERC_LESSHS, 1, 0))
  Mitig_ACS_states[[i]]$EXCEED_LING <- with(Mitig_ACS_states[[i]], ifelse(!is.na(PERC_LING) & !is.na(PERC_LING_WET) &
                                                                            PERC_LING_WET>PERC_LING, 1, 0))
  Mitig_ACS_states[[i]]$EXCEED_UNDER5 <- with(Mitig_ACS_states[[i]], ifelse(!is.na(PERC_UNDER5) & !is.na(PERC_UNDER5_WET) &
                                                                              PERC_UNDER5_WET>PERC_UNDER5, 1, 0))
  Mitig_ACS_states[[i]]$EXCEED_OVER64 <- with(Mitig_ACS_states[[i]], ifelse(!is.na(PERC_OVER64) & !is.na(PERC_OVER64_WET) &
                                                                              PERC_OVER64_WET>PERC_OVER64, 1, 0))
  Mitig_ACS_states[[i]]$EXCEED_LEAD <- with(Mitig_ACS_states[[i]], ifelse(!is.na(PERC_LEAD) & !is.na(PERC_LEAD_WET) &
                                                                            PERC_LEAD_WET>PERC_LEAD, 1, 0))
  
  for (j in 1:length(race_var)) {
    
    Mitig_ACS_states[[i]][, paste0("EXCEED_",race_var[j])] <- 
      ifelse(!is.na(Mitig_ACS_states[[i]][, paste0("PERC_",race_var[j],"_WET")]) & 
               !is.na(Mitig_ACS_states[[i]][, paste0("PERC_",race_var[j])]) & 
               Mitig_ACS_states[[i]][, paste0("PERC_",race_var[j],"_WET")] > 
               Mitig_ACS_states[[i]][, paste0("PERC_",race_var[j])], 1, 0)
    
  }
  
  Mitig_ACS_states_fin[[i]] <- Mitig_ACS_states[[i]] %>%
    select(starts_with("EXCEED")) %>%
    summarise_all(sum, na.rm=T)
  
  Mitig_ACS_states_fin[[i]]$HUC12_Num <- length(unique(Mitig_ACS_quant[[i]]$HUC12))
  Mitig_ACS_states_fin[[i]]$ORM_Total <- length(unique(ORM_mitigation$HUC12))
  Mitig_ACS_states_fin[[i]]$HUC12_Perc <- with(Mitig_ACS_states_fin[[i]], HUC12_Num/ORM_Total)
  
  if(exists("Mitig_ACS_states_table")) {
    Mitig_ACS_states_table <- rbind(Mitig_ACS_states_table, Mitig_ACS_states_fin[[i]])
  } else {
    Mitig_ACS_states_table <- Mitig_ACS_states_fin[[i]]
  }
  
  for (j in 1:length(env_var)) {
    
    Mitig_ACS_states[[i]][, paste0("EXCEED_ENV_",env_var[j])] <- 
      ifelse(Mitig_ACS_states[[i]][, paste0(env_var[j],"_num.x")] > 
               Mitig_ACS_states[[i]][, paste0(env_var[j],"_num.y")], 1, 0)
    
  }
  
  exceed_colname <- paste0("EXCEED_CNT_",i)
  
  Mitig_ACS_env_states[[i]] <- Mitig_ACS_states[[i]] %>%
    mutate(EXCEED_TOTAL = rowSums(across(c(starts_with("EXCEED_ENV"),"EXCEED_LEAD")), na.rm=T)) %>%
    group_by(EXCEED_TOTAL) %>%
    summarize(!!exceed_colname := n_distinct(HUC12))
  
  Mitig_ACS_env_states_table <- merge(Mitig_ACS_env_states_table, Mitig_ACS_env_states[[i]], by.x="Exceedances",
                               by.y="EXCEED_TOTAL", all.x=T)
  
  Mitig_ACS_EJ_totals[[i]] <- data.frame(cbind(t(colSums(Mitig_ACS_EJ_quant[[i]]%>%select(TOT_AGE_POP:PRE1960), na.rm=T)),
                                               t(colMeans(Mitig_ACS_EJ_quant[[i]]%>%select(DSLPM_num:PM25_num), na.rm=T))))
  
  Mitig_ACS_EJ_totals[[i]]$PERC_MIN <- with(Mitig_ACS_EJ_totals[[i]], MINOR_POP/TOT_RACE_POP)
  Mitig_ACS_EJ_totals[[i]]$PERC_POV <- with(Mitig_ACS_EJ_totals[[i]], POV_POP/TOT_POV_POP)
  Mitig_ACS_EJ_totals[[i]]$PERC_LESSHS <- with(Mitig_ACS_EJ_totals[[i]], EDU_POP/TOT_EDU_POP)
  Mitig_ACS_EJ_totals[[i]]$PERC_LING <- with(Mitig_ACS_EJ_totals[[i]], LINGISO_HH/TOT_LINGISO_HH)
  Mitig_ACS_EJ_totals[[i]]$PERC_UNDER5 <- with(Mitig_ACS_EJ_totals[[i]], UNDER_5/TOT_AGE_POP)
  Mitig_ACS_EJ_totals[[i]]$PERC_OVER64 <- with(Mitig_ACS_EJ_totals[[i]], OVER_64/TOT_AGE_POP)
  Mitig_ACS_EJ_totals[[i]]$PERC_LEAD <- with(Mitig_ACS_EJ_totals[[i]], PRE1960/ACSTOTHU)
  
  for (j in 1:length(race_var)) {
    
    Mitig_ACS_EJ_totals[[i]][, paste0("PERC_",race_var[j])] <- 
      Mitig_ACS_EJ_totals[[i]][, race_var[j]] / Mitig_ACS_EJ_totals[[i]]$TOT_RACE_POP
    
  }
  
  Mitig_ACS_EJ_totals[[i]]$HUC12_Num <- length(unique(Mitig_ACS_EJ_quant[[i]]$HUC12))
  Mitig_ACS_EJ_totals[[i]]$ORM_Total <- length(unique(ORM_mitigation$HUC12))
  Mitig_ACS_EJ_totals[[i]]$HUC12_Perc <- with(Mitig_ACS_EJ_totals[[i]], HUC12_Num/ORM_Total)
  
  if(exists("Mitig_ACS_EJ_table")) {
    Mitig_ACS_EJ_table <- rbind(Mitig_ACS_EJ_table, Mitig_ACS_EJ_totals[[i]])
  } else {
    Mitig_ACS_EJ_table <- Mitig_ACS_EJ_totals[[i]]
  }
  
  Mitig_ACS_EJ_states[[i]] <- Mitig_ACS_quant[[i]]
  
  Mitig_ACS_EJ_states[[i]]$PERC_MIN_WET <- with(Mitig_ACS_EJ_states[[i]], MINOR_POP/TOT_RACE_POP)
  Mitig_ACS_EJ_states[[i]]$PERC_POV_WET <- with(Mitig_ACS_EJ_states[[i]], POV_POP/TOT_POV_POP)
  Mitig_ACS_EJ_states[[i]]$PERC_LESSHS_WET <- with(Mitig_ACS_EJ_states[[i]], EDU_POP/TOT_EDU_POP)
  Mitig_ACS_EJ_states[[i]]$PERC_LING_WET <- with(Mitig_ACS_EJ_states[[i]], LINGISO_HH/TOT_LINGISO_HH)
  Mitig_ACS_EJ_states[[i]]$PERC_UNDER5_WET <- with(Mitig_ACS_EJ_states[[i]], UNDER_5/TOT_AGE_POP)
  Mitig_ACS_EJ_states[[i]]$PERC_OVER64_WET <- with(Mitig_ACS_EJ_states[[i]], OVER_64/TOT_AGE_POP)
  Mitig_ACS_EJ_states[[i]]$PERC_LEAD_WET <- with(Mitig_ACS_EJ_states[[i]], PRE1960/ACSTOTHU)
  
  for (j in 1:length(race_var)) {
    
    Mitig_ACS_EJ_states[[i]][, paste0("PERC_",race_var[j],"_WET")] <- 
      Mitig_ACS_EJ_states[[i]][, race_var[j]] / Mitig_ACS_EJ_states[[i]]$TOT_RACE_POP
    
  }
  
  Mitig_ACS_EJ_states[[i]] <- merge(Mitig_ACS_EJ_states[[i]], 
                                 State_stats_env %>% 
                                   select(STATE_ID, PERC_MIN:PERC_HISPANIC, 
                                          PRMP_state=PRMP_num, PWDIS_state=PWDIS_num),
                                 by="STATE_ID", all.x=T)
  
  Mitig_ACS_EJ_states[[i]]$EXCEED_MIN <- with(Mitig_ACS_EJ_states[[i]], ifelse(!is.na(PERC_MIN) & 
                                                                                 !is.na(PERC_MIN_WET) & 
                                                                                 PERC_MIN_WET>PERC_MIN, 1, 0))
  Mitig_ACS_EJ_states[[i]]$EXCEED_POV <- with(Mitig_ACS_EJ_states[[i]], ifelse(!is.na(PERC_POV) & 
                                                                                 !is.na(PERC_POV_WET) & 
                                                                                 PERC_POV_WET>PERC_POV, 1, 0))
  Mitig_ACS_EJ_states[[i]]$EXCEED_LESSHS <- with(Mitig_ACS_EJ_states[[i]], ifelse(!is.na(PERC_LESSHS) & 
                                                                                    !is.na(PERC_LESSHS_WET) & 
                                                                                    PERC_LESSHS_WET>PERC_LESSHS, 1, 0))
  Mitig_ACS_EJ_states[[i]]$EXCEED_LING <- with(Mitig_ACS_EJ_states[[i]], ifelse(!is.na(PERC_LING) & 
                                                                                  !is.na(PERC_LING_WET) & 
                                                                                  PERC_LING_WET>PERC_LING, 1, 0))
  Mitig_ACS_EJ_states[[i]]$EXCEED_UNDER5 <- with(Mitig_ACS_EJ_states[[i]], ifelse(!is.na(PERC_UNDER5) & 
                                                                                    !is.na(PERC_UNDER5_WET) & 
                                                                                    PERC_UNDER5_WET>PERC_UNDER5, 1, 0))
  Mitig_ACS_EJ_states[[i]]$EXCEED_OVER64 <- with(Mitig_ACS_EJ_states[[i]], ifelse(!is.na(PERC_OVER64) & 
                                                                                    !is.na(PERC_OVER64_WET) & 
                                                                                    PERC_OVER64_WET>PERC_OVER64, 1, 0))
  Mitig_ACS_EJ_states[[i]]$EXCEED_LEAD <- with(Mitig_ACS_EJ_states[[i]], ifelse(!is.na(PERC_LEAD) & 
                                                                                  !is.na(PERC_LEAD_WET) & 
                                                                                  PERC_LEAD_WET>PERC_LEAD, 1, 0))
  
  for (j in 1:length(race_var)) {
    
    Mitig_ACS_EJ_states[[i]][, paste0("EXCEED_",race_var[j])] <- 
      ifelse(!is.na(Mitig_ACS_EJ_states[[i]][, paste0("PERC_",race_var[j],"_WET")]) & 
               !is.na(Mitig_ACS_EJ_states[[i]][, paste0("PERC_",race_var[j])]) & 
               Mitig_ACS_EJ_states[[i]][, paste0("PERC_",race_var[j],"_WET")] > 
               Mitig_ACS_EJ_states[[i]][, paste0("PERC_",race_var[j])], 1, 0)
    
  }
  
  Mitig_ACS_EJ_states[[i]]$EXCEED_PRMP <- with(Mitig_ACS_EJ_states[[i]], ifelse(!is.na(PRMP_num) & !is.na(PRMP_state) & 
                                                                                  PRMP_num>PRMP_state, 1, 0))
  Mitig_ACS_EJ_states[[i]]$EXCEED_PWDIS <- with(Mitig_ACS_EJ_states[[i]], ifelse(!is.na(PWDIS_num) & !is.na(PWDIS_state) & 
                                                                                   PWDIS_num>PWDIS_state, 1, 0))
  
  Mitig_ACS_EJ_states_fin[[i]] <- Mitig_ACS_EJ_states[[i]] %>%
    select(starts_with("EXCEED")) %>%
    summarise_all(sum, na.rm=T)
  
  Mitig_ACS_EJ_states_fin[[i]]$HUC12_Num <- length(unique(Mitig_ACS_EJ_quant[[i]]$HUC12))
  Mitig_ACS_EJ_states_fin[[i]]$ORM_Total <- length(unique(ORM_mitigation$HUC12))
  Mitig_ACS_EJ_states_fin[[i]]$HUC12_Perc <- with(Mitig_ACS_EJ_states_fin[[i]], HUC12_Num/ORM_Total)
  
  if(exists("Mitig_ACS_EJ_states_table")) {
    Mitig_ACS_EJ_states_table <- rbind(Mitig_ACS_EJ_states_table, Mitig_ACS_EJ_states_fin[[i]])
  } else {
    Mitig_ACS_EJ_states_table <- Mitig_ACS_EJ_states_fin[[i]]
  }
  
}

# Cumulative impacts table - national
write.csv(Mitig_ACS_env_table,paste0(EJ,"R_programs/Outputs/WetArea_CumImpacts_Natl.csv"),row.names=FALSE)

# Cumulative impacts table - state
write.csv(Mitig_ACS_env_states_table,paste0(EJ,"R_programs/Outputs/WetArea_CumImpacts_State.csv"),row.names=FALSE)

# National comparison table
write.csv(Mitig_ACS_totals_table,paste0(EJ,"R_programs/Outputs/WetArea_Natl.csv"),row.names=FALSE)

# State comparison table
write.csv(Mitig_ACS_states_table,paste0(EJ,"R_programs/Outputs/WetArea_States.csv"),row.names=FALSE)

# National comparison table - Env considerations
write.csv(Mitig_ACS_EJ_table,paste0(EJ,"R_programs/Outputs/WetArea_Natl_Env.csv"),row.names=FALSE)

# State comparison table - Env considerations
write.csv(Mitig_ACS_EJ_states_table,paste0(EJ,"R_programs/Outputs/WetArea_State_Env.csv"),row.names=FALSE)

#*# SE: Mitig_ACS_totals do not sum to ORM_mitig_ACS totals
sum(Mitig_ACS_totals[[1]]$HUC12_Num, Mitig_ACS_totals[[2]]$HUC12_Num, Mitig_ACS_totals[[3]]$HUC12_Num,
Mitig_ACS_totals[[4]]$HUC12_Num) # No [[5]] list
# 43247, matches IW_ACS_Totals number

length(unique(ORM_mitig_ACS$HUC12))
# 97282

#*# AL: 43247 are just those with wetland changes or affected waters. This is okay.

## Impacted waters
IW_ACS_quant <- list()
IW_ACS_env <- list()
IW_ACS_totals <- list()
IW_ACS_env_states <- list()
IW_ACS_states <- list()
IW_ACS_states_fin <- list()

IW_ACS_EJ_quant <- list()
IW_ACS_EJ_totals <- list()
IW_ACS_EJ_states <- list()
IW_ACS_EJ_states_fin <- list()

IW_ACS_env_table <- data.frame(c(0:11))
names(IW_ACS_env_table) <- "Exceedances"

IW_ACS_env_states_table <- data.frame(c(0:11))
names(IW_ACS_env_states_table) <- "Exceedances"

for (i in 1:4) {
  # Splitting up the data into the bins agreed upon with EPA
  if (i==1) {
    IW_ACS_quant[[i]] <- ORM_IW_ACS %>%
      dplyr::filter(!is.na(WATERS_CHANGE) & WATERS_CHANGE<=water_break[i])
    
    IW_ACS_EJ_quant[[i]] <- ORM_IW_ACS %>%
      dplyr::filter(!is.na(WATERS_CHANGE) & WATERS_CHANGE<=water_break[i] & 
                      !is.na(PRMP_num) & PRMP_num>0 & 
                      !is.na(PWDIS_num) & PWDIS_num>0)
    
  } else if (i<4) {
    IW_ACS_quant[[i]] <- ORM_IW_ACS %>%
      dplyr::filter(!is.na(WATERS_CHANGE) & water_break[i-1]<WATERS_CHANGE &
                      WATERS_CHANGE<=water_break[i])
    
    IW_ACS_EJ_quant[[i]] <- ORM_IW_ACS %>%
      dplyr::filter(!is.na(WATERS_CHANGE) & water_break[i-1]<WATERS_CHANGE &
                      WATERS_CHANGE<=water_break[i] & 
                      !is.na(PRMP_num) & PRMP_num>0 & 
                      !is.na(PWDIS_num) & PWDIS_num>0)
    
  } else if (i==4) {
    IW_ACS_quant[[i]] <- ORM_IW_ACS %>%
      dplyr::filter(!is.na(WATERS_CHANGE) & WATERS_CHANGE>water_break[i-1])
    
    IW_ACS_EJ_quant[[i]] <- ORM_IW_ACS %>%
      dplyr::filter(!is.na(WATERS_CHANGE) & WATERS_CHANGE>water_break[i-1] & 
                      !is.na(PRMP_num) & PRMP_num>0 & 
                      !is.na(PWDIS_num) & PWDIS_num>0)
    
  }
  
  for (j in 1:length(env_var)) {
    
    IW_ACS_quant[[i]][, paste0("EXCEED_",env_var[j])] <- 
      ifelse(IW_ACS_quant[[i]][, paste0(env_var[j],"_num")] > 
               Natl_stats[, paste0(env_var[j],"_num")], 1, 0)
    
    jpeg(paste0(EJ,"R_programs/Outputs/IW_Natl_Exceedances_",env_var[j],"_",i,".jpg"))
    hist(IW_ACS_quant[[i]][, paste0(env_var[j],"_num")])
    dev.off()
    
  }
  
  IW_ACS_quant[[i]]$PERC_LEAD_ENV <- with(IW_ACS_quant[[i]], PRE1960/ACSTOTHU)
  
  jpeg(paste0(EJ,"R_programs/Outputs/IW_Natl_Exceedances_PERC_LEAD_",i,".jpg"))
  hist(IW_ACS_quant[[i]]$PERC_LEAD_ENV)
  dev.off()
  
  IW_ACS_quant[[i]]$EXCEED_LEAD <- with(IW_ACS_quant[[i]], 
                                        ifelse(PERC_LEAD_ENV > Natl_stats$PERC_LEAD, 1, 0))
  
  exceed_colname <- paste0("EXCEED_CNT_",i)
  
  IW_ACS_env[[i]] <- IW_ACS_quant[[i]] %>%
    mutate(EXCEED_TOTAL = rowSums(across(starts_with("EXCEED")), na.rm=T)) %>%
    group_by(EXCEED_TOTAL) %>%
    summarize(!!exceed_colname := n_distinct(HUC12))
  
  IW_ACS_env_table <- merge(IW_ACS_env_table, IW_ACS_env[[i]], by.x="Exceedances",
                            by.y="EXCEED_TOTAL", all.x=T)
  IW_ACS_totals[[i]] <- data.frame(cbind(t(colSums(IW_ACS_quant[[i]]%>%select(TOT_AGE_POP:PRE1960), na.rm=T)),
                                         t(colMeans(IW_ACS_quant[[i]]%>%select(DSLPM_num:PM25_num), na.rm=T))))
  
  IW_ACS_totals[[i]]$PERC_MIN <- with(IW_ACS_totals[[i]], MINOR_POP/TOT_RACE_POP)
  IW_ACS_totals[[i]]$PERC_POV <- with(IW_ACS_totals[[i]], POV_POP/TOT_POV_POP)
  IW_ACS_totals[[i]]$PERC_LESSHS <- with(IW_ACS_totals[[i]], EDU_POP/TOT_EDU_POP)
  IW_ACS_totals[[i]]$PERC_LING <- with(IW_ACS_totals[[i]], LINGISO_HH/TOT_LINGISO_HH)
  IW_ACS_totals[[i]]$PERC_UNDER5 <- with(IW_ACS_totals[[i]], UNDER_5/TOT_AGE_POP)
  IW_ACS_totals[[i]]$PERC_OVER64 <- with(IW_ACS_totals[[i]], OVER_64/TOT_AGE_POP)
  IW_ACS_totals[[i]]$PERC_LEAD <- with(IW_ACS_totals[[i]], PRE1960/ACSTOTHU)
  
  for (j in 1:length(race_var)) {
    
    IW_ACS_totals[[i]][, paste0("PERC_",race_var[j])] <- 
      IW_ACS_totals[[i]][, race_var[j]] / IW_ACS_totals[[i]]$TOT_RACE_POP
    
  }
  
  IW_ACS_totals[[i]]$HUC12_Num <- length(unique(IW_ACS_quant[[i]]$HUC12))
  IW_ACS_totals[[i]]$ORM_Total <- length(unique(ORM_impactedwaters$HUC12))
  IW_ACS_totals[[i]]$HUC12_Perc <- with(IW_ACS_totals[[i]], HUC12_Num/ORM_Total)
  
  if(exists("IW_ACS_totals_table")) {
    IW_ACS_totals_table <- rbind(IW_ACS_totals_table, IW_ACS_totals[[i]])
  } else {
    IW_ACS_totals_table <- IW_ACS_totals[[i]]
  }
  
  IW_ACS_states[[i]] <- IW_ACS_quant[[i]]
  
  IW_ACS_states[[i]]$PERC_MIN_WET <- with(IW_ACS_states[[i]], MINOR_POP/TOT_RACE_POP)
  IW_ACS_states[[i]]$PERC_POV_WET <- with(IW_ACS_states[[i]], POV_POP/TOT_POV_POP)
  IW_ACS_states[[i]]$PERC_LESSHS_WET <- with(IW_ACS_states[[i]], EDU_POP/TOT_EDU_POP)
  IW_ACS_states[[i]]$PERC_LING_WET <- with(IW_ACS_states[[i]], LINGISO_HH/TOT_LINGISO_HH)
  IW_ACS_states[[i]]$PERC_UNDER5_WET <- with(IW_ACS_states[[i]], UNDER_5/TOT_AGE_POP)
  IW_ACS_states[[i]]$PERC_OVER64_WET <- with(IW_ACS_states[[i]], OVER_64/TOT_AGE_POP)
  IW_ACS_states[[i]]$PERC_LEAD_WET <- with(IW_ACS_states[[i]], PRE1960/ACSTOTHU)
  
  for (j in 1:length(race_var)) {
    
    IW_ACS_states[[i]][, paste0("PERC_",race_var[j],"_WET")] <- 
      IW_ACS_states[[i]][, race_var[j]] / IW_ACS_states[[i]]$TOT_RACE_POP
    
  }
  
  IW_ACS_states[[i]] <- merge(IW_ACS_states[[i]], 
                              State_stats %>% select(STATE_ID, PERC_MIN:PM25_num), by="STATE_ID", all.x=T)
  
  IW_ACS_states[[i]]$EXCEED_MIN <- with(IW_ACS_states[[i]], ifelse(!is.na(PERC_MIN) & !is.na(PERC_MIN_WET) & 
                                                                     PERC_MIN_WET>PERC_MIN, 1, 0))
  IW_ACS_states[[i]]$EXCEED_POV <- with(IW_ACS_states[[i]], ifelse(!is.na(PERC_POV) & !is.na(PERC_POV_WET) &
                                                                     PERC_POV_WET>PERC_POV, 1, 0))
  IW_ACS_states[[i]]$EXCEED_LESSHS <- with(IW_ACS_states[[i]], ifelse(!is.na(PERC_LESSHS) & !is.na(PERC_LESSHS_WET) &
                                                                        PERC_LESSHS_WET>PERC_LESSHS, 1, 0))
  IW_ACS_states[[i]]$EXCEED_LING <- with(IW_ACS_states[[i]], ifelse(!is.na(PERC_LING) & !is.na(PERC_LING_WET) &
                                                                      PERC_LING_WET>PERC_LING, 1, 0))
  IW_ACS_states[[i]]$EXCEED_UNDER5 <- with(IW_ACS_states[[i]], ifelse(!is.na(PERC_UNDER5) & !is.na(PERC_UNDER5_WET) &
                                                                        PERC_UNDER5_WET>PERC_UNDER5, 1, 0))
  IW_ACS_states[[i]]$EXCEED_OVER64 <- with(IW_ACS_states[[i]], ifelse(!is.na(PERC_OVER64) & !is.na(PERC_OVER64_WET) &
                                                                        PERC_OVER64_WET>PERC_OVER64, 1, 0))
  IW_ACS_states[[i]]$EXCEED_LEAD <- with(IW_ACS_states[[i]], ifelse(!is.na(PERC_LEAD) & !is.na(PERC_LEAD_WET) &
                                                                      PERC_LEAD_WET>PERC_LEAD, 1, 0))
  
  for (j in 1:length(race_var)) {
    
    IW_ACS_states[[i]][, paste0("EXCEED_",race_var[j])] <- 
      ifelse(!is.na(IW_ACS_states[[i]][, paste0("PERC_",race_var[j],"_WET")]) & 
               !is.na(IW_ACS_states[[i]][, paste0("PERC_",race_var[j])]) & 
               IW_ACS_states[[i]][, paste0("PERC_",race_var[j],"_WET")] > 
               IW_ACS_states[[i]][, paste0("PERC_",race_var[j])], 1, 0)
    
  }
  
  IW_ACS_states_fin[[i]] <- IW_ACS_states[[i]] %>%
    select(starts_with("EXCEED")) %>%
    summarise_all(sum, na.rm=T)
  
  IW_ACS_states_fin[[i]]$HUC12_Num <- length(unique(IW_ACS_quant[[i]]$HUC12))
  IW_ACS_states_fin[[i]]$ORM_Total <- length(unique(ORM_impactedwaters$HUC12))
  IW_ACS_states_fin[[i]]$HUC12_Perc <- with(IW_ACS_states_fin[[i]], HUC12_Num/ORM_Total)
  
  if(exists("IW_ACS_states_table")) {
    IW_ACS_states_table <- rbind(IW_ACS_states_table, IW_ACS_states_fin[[i]])
  } else {
    IW_ACS_states_table <- IW_ACS_states_fin[[i]]
  }
  
  for (j in 1:length(env_var)) {
    
    IW_ACS_states[[i]][, paste0("EXCEED_ENV_",env_var[j])] <- 
      ifelse(IW_ACS_states[[i]][, paste0(env_var[j],"_num.x")] > 
               IW_ACS_states[[i]][, paste0(env_var[j],"_num.y")], 1, 0)
    
  }
  
  exceed_colname <- paste0("EXCEED_CNT_",i)
  
  IW_ACS_env_states[[i]] <- IW_ACS_states[[i]] %>%
    mutate(EXCEED_TOTAL = rowSums(across(c(starts_with("EXCEED_ENV"),"EXCEED_LEAD")), na.rm=T)) %>%
    group_by(EXCEED_TOTAL) %>%
    summarize(!!exceed_colname := n_distinct(HUC12))
  
  IW_ACS_env_states_table <- merge(IW_ACS_env_states_table, IW_ACS_env_states[[i]], by.x="Exceedances",
                                   by.y="EXCEED_TOTAL", all.x=T)
  
  IW_ACS_EJ_totals[[i]] <- data.frame(cbind(t(colSums(IW_ACS_EJ_quant[[i]]%>%select(TOT_AGE_POP:PRE1960), na.rm=T)),
                                            t(colMeans(IW_ACS_EJ_quant[[i]]%>%select(DSLPM_num:PM25_num), na.rm=T))))
  
  IW_ACS_EJ_totals[[i]]$PERC_MIN <- with(IW_ACS_EJ_totals[[i]], MINOR_POP/TOT_RACE_POP)
  IW_ACS_EJ_totals[[i]]$PERC_POV <- with(IW_ACS_EJ_totals[[i]], POV_POP/TOT_POV_POP)
  IW_ACS_EJ_totals[[i]]$PERC_LESSHS <- with(IW_ACS_EJ_totals[[i]], EDU_POP/TOT_EDU_POP)
  IW_ACS_EJ_totals[[i]]$PERC_LING <- with(IW_ACS_EJ_totals[[i]], LINGISO_HH/TOT_LINGISO_HH)
  IW_ACS_EJ_totals[[i]]$PERC_UNDER5 <- with(IW_ACS_EJ_totals[[i]], UNDER_5/TOT_AGE_POP)
  IW_ACS_EJ_totals[[i]]$PERC_OVER64 <- with(IW_ACS_EJ_totals[[i]], OVER_64/TOT_AGE_POP)
  IW_ACS_EJ_totals[[i]]$PERC_LEAD <- with(IW_ACS_EJ_totals[[i]], PRE1960/ACSTOTHU)
  
  for (j in 1:length(race_var)) {
    
    IW_ACS_EJ_totals[[i]][, paste0("PERC_",race_var[j])] <- 
      IW_ACS_EJ_totals[[i]][, race_var[j]] / IW_ACS_EJ_totals[[i]]$TOT_RACE_POP
    
  }
  
  IW_ACS_EJ_totals[[i]]$HUC12_Num <- length(unique(IW_ACS_EJ_quant[[i]]$HUC12))
  IW_ACS_EJ_totals[[i]]$ORM_Total <- length(unique(ORM_impactedwaters$HUC12))
  IW_ACS_EJ_totals[[i]]$HUC12_Perc <- with(IW_ACS_EJ_totals[[i]], HUC12_Num/ORM_Total)
  
  if(exists("IW_ACS_EJ_table")) {
    IW_ACS_EJ_table <- rbind(IW_ACS_EJ_table, IW_ACS_EJ_totals[[i]])
  } else {
    IW_ACS_EJ_table <- IW_ACS_EJ_totals[[i]]
  }
  
  IW_ACS_EJ_states[[i]] <- IW_ACS_quant[[i]]
  
  IW_ACS_EJ_states[[i]]$PERC_MIN_WET <- with(IW_ACS_EJ_states[[i]], MINOR_POP/TOT_RACE_POP)
  IW_ACS_EJ_states[[i]]$PERC_POV_WET <- with(IW_ACS_EJ_states[[i]], POV_POP/TOT_POV_POP)
  IW_ACS_EJ_states[[i]]$PERC_LESSHS_WET <- with(IW_ACS_EJ_states[[i]], EDU_POP/TOT_EDU_POP)
  IW_ACS_EJ_states[[i]]$PERC_LING_WET <- with(IW_ACS_EJ_states[[i]], LINGISO_HH/TOT_LINGISO_HH)
  IW_ACS_EJ_states[[i]]$PERC_UNDER5_WET <- with(IW_ACS_EJ_states[[i]], UNDER_5/TOT_AGE_POP)
  IW_ACS_EJ_states[[i]]$PERC_OVER64_WET <- with(IW_ACS_EJ_states[[i]], OVER_64/TOT_AGE_POP)
  IW_ACS_EJ_states[[i]]$PERC_LEAD_WET <- with(IW_ACS_EJ_states[[i]], PRE1960/ACSTOTHU)
  
  for (j in 1:length(race_var)) {
    
    IW_ACS_EJ_states[[i]][, paste0("PERC_",race_var[j],"_WET")] <- 
      IW_ACS_EJ_states[[i]][, race_var[j]] / IW_ACS_EJ_states[[i]]$TOT_RACE_POP
    
  }
  
  IW_ACS_EJ_states[[i]] <- merge(IW_ACS_EJ_states[[i]], 
                                 State_stats_env %>% 
                                   select(STATE_ID, PERC_MIN:PERC_HISPANIC, 
                                          PRMP_state=PRMP_num, PWDIS_state=PWDIS_num),
                                 by="STATE_ID", all.x=T)
  
  IW_ACS_EJ_states[[i]]$EXCEED_MIN <- with(IW_ACS_EJ_states[[i]], ifelse(!is.na(PERC_MIN) & 
                                                                           !is.na(PERC_MIN_WET) & 
                                                                           PERC_MIN_WET>PERC_MIN, 1, 0))
  IW_ACS_EJ_states[[i]]$EXCEED_POV <- with(IW_ACS_EJ_states[[i]], ifelse(!is.na(PERC_POV) & 
                                                                           !is.na(PERC_POV_WET) & 
                                                                           PERC_POV_WET>PERC_POV, 1, 0))
  IW_ACS_EJ_states[[i]]$EXCEED_LESSHS <- with(IW_ACS_EJ_states[[i]], ifelse(!is.na(PERC_LESSHS) & 
                                                                              !is.na(PERC_LESSHS_WET) & 
                                                                              PERC_LESSHS_WET>PERC_LESSHS, 1, 0))
  IW_ACS_EJ_states[[i]]$EXCEED_LING <- with(IW_ACS_EJ_states[[i]], ifelse(!is.na(PERC_LING) & 
                                                                            !is.na(PERC_LING_WET) & 
                                                                            PERC_LING_WET>PERC_LING, 1, 0))
  IW_ACS_EJ_states[[i]]$EXCEED_UNDER5 <- with(IW_ACS_EJ_states[[i]], ifelse(!is.na(PERC_UNDER5) & 
                                                                              !is.na(PERC_UNDER5_WET) & 
                                                                              PERC_UNDER5_WET>PERC_UNDER5, 1, 0))
  IW_ACS_EJ_states[[i]]$EXCEED_OVER64 <- with(IW_ACS_EJ_states[[i]], ifelse(!is.na(PERC_OVER64) & 
                                                                              !is.na(PERC_OVER64_WET) & 
                                                                              PERC_OVER64_WET>PERC_OVER64, 1, 0))
  IW_ACS_EJ_states[[i]]$EXCEED_LEAD <- with(IW_ACS_EJ_states[[i]], ifelse(!is.na(PERC_LEAD) & 
                                                                            !is.na(PERC_LEAD_WET) & 
                                                                            PERC_LEAD_WET>PERC_LEAD, 1, 0))
  
  for (j in 1:length(race_var)) {
    
    IW_ACS_EJ_states[[i]][, paste0("EXCEED_",race_var[j])] <- 
      ifelse(!is.na(IW_ACS_EJ_states[[i]][, paste0("PERC_",race_var[j],"_WET")]) & 
               !is.na(IW_ACS_EJ_states[[i]][, paste0("PERC_",race_var[j])]) & 
               IW_ACS_EJ_states[[i]][, paste0("PERC_",race_var[j],"_WET")] > 
               IW_ACS_EJ_states[[i]][, paste0("PERC_",race_var[j])], 1, 0)
    
  }
  
  IW_ACS_EJ_states[[i]]$EXCEED_PRMP <- with(IW_ACS_EJ_states[[i]], ifelse(!is.na(PRMP_num) & !is.na(PRMP_state) & 
                                                                            PRMP_num>PRMP_state, 1, 0))
  IW_ACS_EJ_states[[i]]$EXCEED_PWDIS <- with(IW_ACS_EJ_states[[i]], ifelse(!is.na(PWDIS_num) & !is.na(PWDIS_state) & 
                                                                             PWDIS_num>PWDIS_state, 1, 0))
  
  IW_ACS_EJ_states_fin[[i]] <- IW_ACS_EJ_states[[i]] %>%
    select(starts_with("EXCEED")) %>%
    summarise_all(sum, na.rm=T)
  
  IW_ACS_EJ_states_fin[[i]]$HUC12_Num <- length(unique(IW_ACS_EJ_quant[[i]]$HUC12))
  IW_ACS_EJ_states_fin[[i]]$ORM_Total <- length(unique(ORM_impactedwaters$HUC12))
  IW_ACS_EJ_states_fin[[i]]$HUC12_Perc <- with(IW_ACS_EJ_states_fin[[i]], HUC12_Num/ORM_Total)
  
  if(exists("IW_ACS_EJ_states_table")) {
    IW_ACS_EJ_states_table <- rbind(IW_ACS_EJ_states_table, IW_ACS_EJ_states_fin[[i]])
  } else {
    IW_ACS_EJ_states_table <- IW_ACS_EJ_states_fin[[i]]
  }
  
}

# Cumulative impacts table - national
write.csv(IW_ACS_env_table,paste0(EJ,"R_programs/Outputs/IW_CumImpacts_Natl.csv"),row.names=FALSE)

# Cumulative impacts table - state
write.csv(IW_ACS_env_states_table,paste0(EJ,"R_programs/Outputs/IW_CumImpacts_State.csv"),row.names=FALSE)

# National comparison table
write.csv(IW_ACS_totals_table,paste0(EJ,"R_programs/Outputs/IW_Natl.csv"),row.names=FALSE)

# State comparison table
write.csv(IW_ACS_states_table,paste0(EJ,"R_programs/Outputs/IW_States.csv"),row.names=FALSE)

# National comparison table - Env considerations
write.csv(IW_ACS_EJ_table,paste0(EJ,"R_programs/Outputs/IW_Natl_Env.csv"),row.names=FALSE)

# State comparison table - Env considerations
write.csv(IW_ACS_EJ_states_table,paste0(EJ,"R_programs/Outputs/IW_State_Env.csv"),row.names=FALSE)

#*# SE: There is no IW_ACS_totals[[5]], totals do not match
sum(IW_ACS_totals[[1]]$HUC12_Num, IW_ACS_totals[[2]]$HUC12_Num, IW_ACS_totals[[3]]$HUC12_Num, IW_ACS_totals[[4]]$HUC12_Num)
# 43247, matches Mitigs_ACS_Totals number

length(unique(ORM_IW_ACS$HUC12))
# 97282

#*# AL: 43247 are just those with wetland changes or affected waters. This is okay.

# PART 6: Read in and merge with water quality changes -------------------------

for (i in 1:length(watersheds)) {
  
  WQ_changes <- read_xlsx(paste0(WQ,substr(watersheds[i],1,4),"/",watersheds[i]),
                          col_types = c("text","numeric","text",rep("numeric",6)))
  WQ_changes$Subbasin <- with(WQ_changes, ifelse(nchar(Subbasin)==11,paste0("0",Subbasin),Subbasin))
  
  print(summary(WQ_changes$WQI))
  
  WQ_changes_ACS <- merge(WQ_changes, BG_HUC12_Final, by.x="Subbasin", by.y="HUC12", all.x=T)
  
  # PART 7: Bin the outputs based on water quality changes ---------------------
  
  WQ_changes_ACS_sub <- list()
  WQ_changes_ACS_env <- list()
  WQ_changes_ACS_env_states <- list()
  WQ_changes_ACS_totals <- list()
  WQ_changes_ACS_states <- list()
  WQ_changes_ACS_states_fin <- list()
  
  WQ_changes_ACS_env_table <- data.frame(c(0:11))
  names(WQ_changes_ACS_env_table) <- "Exceedances"
  
  WQ_changes_ACS_env_states_table <- data.frame(c(0:11))
  names(WQ_changes_ACS_env_states_table) <- "Exceedances"
  
  for (k in 1:3) {  
    
    if (k==1) {
      WQ_changes_ACS_sub[[k]] <- WQ_changes_ACS %>%
        dplyr::filter(WQI<0)
      
    } else if (k==2) {
      WQ_changes_ACS_sub[[k]] <- WQ_changes_ACS %>%
        dplyr::filter(WQI==0)
      
    } else if (k==3) {
      WQ_changes_ACS_sub[[k]] <- WQ_changes_ACS %>%
        dplyr::filter(WQI>0)
      
    }
    
    for (m in 1:length(env_var)) {
      
      WQ_changes_ACS_sub[[k]][, paste0("EXCEED_",env_var[m])] <- 
        ifelse(WQ_changes_ACS_sub[[k]][, paste0(env_var[m],"_num")] > 
                 Natl_stats[, paste0(env_var[m],"_num")], 1, 0)
      
      jpeg(paste0(EJ,"R_programs/Outputs/WQchanges_Natl_Exceedances_",substr(watersheds[i],1,4),
                  env_var[j],"_",i,".jpg"))
      hist(WQ_changes_ACS_sub[[k]][, paste0(env_var[j],"_num")])
      dev.off()
      
    }
    
    WQ_changes_ACS_sub[[k]]$PERC_LEAD_ENV <- with(WQ_changes_ACS_sub[[k]], PRE1960/ACSTOTHU)
    
    jpeg(paste0(EJ,"R_programs/Outputs/WQchanges_Natl_Exceedances_",substr(watersheds[i],1,4),
                "PERC_LEAD_",i,".jpg"))
    hist(WQ_changes_ACS_sub[[k]]$PERC_LEAD_ENV)
    dev.off()
    
    WQ_changes_ACS_sub[[k]]$EXCEED_LEAD <- with(WQ_changes_ACS_sub[[k]], 
                                          ifelse(PERC_LEAD_ENV > Natl_stats$PERC_LEAD, 1, 0))
    
    exceed_colname <- paste0("EXCEED_CNT_",k)
    
    WQ_changes_ACS_env[[k]] <- WQ_changes_ACS_sub[[k]] %>%
      mutate(EXCEED_TOTAL = rowSums(across(starts_with("EXCEED")), na.rm=T)) %>%
      group_by(EXCEED_TOTAL) %>%
      summarize(!!exceed_colname := n_distinct(Subbasin))
    
    WQ_changes_ACS_env_table <- merge(WQ_changes_ACS_env_table, WQ_changes_ACS_env[[k]], by.x="Exceedances",
                                      by.y="EXCEED_TOTAL", all.x=T)
    
    WQ_changes_ACS_totals[[k]] <- data.frame(cbind(t(colSums(WQ_changes_ACS_sub[[k]]%>%select(TOT_AGE_POP:PRE1960), na.rm=T)),
                                                   t(colMeans(WQ_changes_ACS_sub[[k]]%>%select(DSLPM_num:PM25_num), na.rm=T))))
    
    WQ_changes_ACS_totals[[k]]$PERC_MIN <- with(WQ_changes_ACS_totals[[k]], MINOR_POP/TOT_RACE_POP)
    WQ_changes_ACS_totals[[k]]$PERC_POV <- with(WQ_changes_ACS_totals[[k]], POV_POP/TOT_POV_POP)
    WQ_changes_ACS_totals[[k]]$PERC_LESSHS <- with(WQ_changes_ACS_totals[[k]], EDU_POP/TOT_EDU_POP)
    WQ_changes_ACS_totals[[k]]$PERC_LING <- with(WQ_changes_ACS_totals[[k]], LINGISO_HH/TOT_LINGISO_HH)
    WQ_changes_ACS_totals[[k]]$PERC_UNDER5 <- with(WQ_changes_ACS_totals[[k]], UNDER_5/TOT_AGE_POP)
    WQ_changes_ACS_totals[[k]]$PERC_OVER64 <- with(WQ_changes_ACS_totals[[k]], OVER_64/TOT_AGE_POP)
    WQ_changes_ACS_totals[[k]]$PERC_LEAD <- with(WQ_changes_ACS_totals[[k]], PRE1960/ACSTOTHU)
    
    for (m in 1:length(race_var)) {
      
      WQ_changes_ACS_totals[[k]][, paste0("PERC_",race_var[m])] <- 
        WQ_changes_ACS_totals[[k]][, race_var[m]] / WQ_changes_ACS_totals[[k]]$TOT_RACE_POP
      
    }
    
    WQ_changes_ACS_totals[[k]]$HUC12_Num <- length(unique(WQ_changes_ACS_sub[[k]]$Subbasin))
    WQ_changes_ACS_totals[[k]]$HUC4_Total <- length(unique(WQ_changes_ACS$Subbasin))
    WQ_changes_ACS_totals[[k]]$HUC12_Perc <- with(WQ_changes_ACS_totals[[k]], HUC12_Num/HUC4_Total)
    
    if(exists("WQ_changes_ACS_totals_table")) {
      WQ_changes_ACS_totals_table <- rbind(WQ_changes_ACS_totals_table, WQ_changes_ACS_totals[[k]])
    } else {
      WQ_changes_ACS_totals_table <- WQ_changes_ACS_totals[[k]]
    }
    
    WQ_changes_ACS_states[[k]] <- WQ_changes_ACS_sub[[k]]
    
    WQ_changes_ACS_states[[k]]$PERC_MIN_WET <- with(WQ_changes_ACS_states[[k]], MINOR_POP/TOT_RACE_POP)
    WQ_changes_ACS_states[[k]]$PERC_POV_WET <- with(WQ_changes_ACS_states[[k]], POV_POP/TOT_POV_POP)
    WQ_changes_ACS_states[[k]]$PERC_LESSHS_WET <- with(WQ_changes_ACS_states[[k]], EDU_POP/TOT_EDU_POP)
    WQ_changes_ACS_states[[k]]$PERC_LING_WET <- with(WQ_changes_ACS_states[[k]], LINGISO_HH/TOT_LINGISO_HH)
    WQ_changes_ACS_states[[k]]$PERC_UNDER5_WET <- with(WQ_changes_ACS_states[[k]], UNDER_5/TOT_AGE_POP)
    WQ_changes_ACS_states[[k]]$PERC_OVER64_WET <- with(WQ_changes_ACS_states[[k]], OVER_64/TOT_AGE_POP)
    WQ_changes_ACS_states[[k]]$PERC_LEAD_WET <- with(WQ_changes_ACS_states[[k]], PRE1960/ACSTOTHU)
    
    for (m in 1:length(race_var)) {
      
      WQ_changes_ACS_states[[k]][, paste0("PERC_",race_var[m],"_WET")] <- 
        WQ_changes_ACS_states[[k]][, race_var[m]] / WQ_changes_ACS_states[[k]]$TOT_RACE_POP
      
    }
    
    WQ_changes_ACS_states[[k]] <- merge(WQ_changes_ACS_states[[k]], 
                                        State_stats %>% select(STATE_ID, PERC_MIN:PM25_num), by="STATE_ID", all.x=T)
    
    WQ_changes_ACS_states[[k]]$EXCEED_MIN <- with(WQ_changes_ACS_states[[k]], ifelse(!is.na(PERC_MIN) & !is.na(PERC_MIN_WET) & 
                                                                                       PERC_MIN_WET>PERC_MIN, 1, 0))
    WQ_changes_ACS_states[[k]]$EXCEED_POV <- with(WQ_changes_ACS_states[[k]], ifelse(!is.na(PERC_POV) & !is.na(PERC_POV_WET) &
                                                                                       PERC_POV_WET>PERC_POV, 1, 0))
    WQ_changes_ACS_states[[k]]$EXCEED_LESSHS <- with(WQ_changes_ACS_states[[k]], ifelse(!is.na(PERC_LESSHS) & !is.na(PERC_LESSHS_WET) &
                                                                                          PERC_LESSHS_WET>PERC_LESSHS, 1, 0))
    WQ_changes_ACS_states[[k]]$EXCEED_LING <- with(WQ_changes_ACS_states[[k]], ifelse(!is.na(PERC_LING) & !is.na(PERC_LING_WET) &
                                                                                        PERC_LING_WET>PERC_LING, 1, 0))
    WQ_changes_ACS_states[[k]]$EXCEED_UNDER5 <- with(WQ_changes_ACS_states[[k]], ifelse(!is.na(PERC_UNDER5) & !is.na(PERC_UNDER5_WET) &
                                                                                          PERC_UNDER5_WET>PERC_UNDER5, 1, 0))
    WQ_changes_ACS_states[[k]]$EXCEED_OVER64 <- with(WQ_changes_ACS_states[[k]], ifelse(!is.na(PERC_OVER64) & !is.na(PERC_OVER64_WET) &
                                                                                          PERC_OVER64_WET>PERC_OVER64, 1, 0))
    WQ_changes_ACS_states[[k]]$EXCEED_LEAD <- with(WQ_changes_ACS_states[[k]], ifelse(!is.na(PERC_LEAD) & !is.na(PERC_LEAD_WET) &
                                                                                        PERC_LEAD_WET>PERC_LEAD, 1, 0))
    
    for (m in 1:length(race_var)) {
      
      WQ_changes_ACS_states[[k]][, paste0("EXCEED_",race_var[m])] <- 
        ifelse(!is.na(WQ_changes_ACS_states[[k]][, paste0("PERC_",race_var[m],"_WET")]) & 
                 !is.na(WQ_changes_ACS_states[[k]][, paste0("PERC_",race_var[m])]) & 
                 WQ_changes_ACS_states[[k]][, paste0("PERC_",race_var[m],"_WET")] > 
                 WQ_changes_ACS_states[[k]][, paste0("PERC_",race_var[m])], 1, 0)
      
    }
    
    WQ_changes_ACS_states_fin[[k]] <- WQ_changes_ACS_states[[k]] %>%
      select(starts_with("EXCEED")) %>%
      summarise_all(sum, na.rm=T)
    
    WQ_changes_ACS_states_fin[[k]]$HUC12_Num <- length(unique(WQ_changes_ACS_sub[[k]]$Subbasin))
    WQ_changes_ACS_states_fin[[k]]$HUC4_Total <- length(unique(WQ_changes_ACS$Subbasin))
    WQ_changes_ACS_states_fin[[k]]$HUC12_Perc <- with(WQ_changes_ACS_states_fin[[k]], HUC12_Num/HUC4_Total)
    
    if(exists("WQ_changes_ACS_states_table")) {
      WQ_changes_ACS_states_table <- rbind(WQ_changes_ACS_states_table, WQ_changes_ACS_states_fin[[k]])
    } else {
      WQ_changes_ACS_states_table <- WQ_changes_ACS_states_fin[[k]]
    }
    
    for (j in 1:length(env_var)) {
      
      WQ_changes_ACS_states[[k]][, paste0("EXCEED_ENV_",env_var[j])] <- 
        ifelse(WQ_changes_ACS_states[[k]][, paste0(env_var[j],"_num.x")] > 
                 WQ_changes_ACS_states[[k]][, paste0(env_var[j],"_num.y")], 1, 0)
      
    }
    
    exceed_colname <- paste0("EXCEED_CNT_",k)
    
    WQ_changes_ACS_env_states[[k]] <- WQ_changes_ACS_states[[k]] %>%
      mutate(EXCEED_TOTAL = rowSums(across(c(starts_with("EXCEED_ENV"),"EXCEED_LEAD")), na.rm=T)) %>%
      group_by(EXCEED_TOTAL) %>%
      summarize(!!exceed_colname := n_distinct(Subbasin))
    
    WQ_changes_ACS_env_states_table <- merge(WQ_changes_ACS_env_states_table, WQ_changes_ACS_env_states[[k]], by.x="Exceedances",
                                     by.y="EXCEED_TOTAL", all.x=T)
    
  }
  
  # Cumulative impacts table - national
  write.csv(WQ_changes_ACS_env_table,
            paste0(EJ,"R_programs/Outputs/WetArea_WQI_",
                   substr(watersheds[i],1,4),"_CumImpacts_Natl.csv"),row.names=FALSE)
  
  # Cumulative impacts table - state
  write.csv(WQ_changes_ACS_env_states_table,
            paste0(EJ, "R_programs/Outputs/WetArea_WQI_",
                   substr(watersheds[i],1,4),"_CumImpacts_States.csv"),row.names=F)
  
  # National comparison table
  write.csv(WQ_changes_ACS_totals_table,
            paste0(EJ,"R_programs/Outputs/WetArea_WQI_",
                   substr(watersheds[i],1,4),"_Natl.csv"),row.names=FALSE)
  
  # State comparison table
  write.csv(WQ_changes_ACS_states_table,
            paste0(EJ,"R_programs/Outputs/WetArea_WQI_",
                   substr(watersheds[i],1,4),"_States.csv"),row.names=FALSE)
  
  rm(WQ_changes_ACS_env_table, WQ_changes_ACS_env_states_table,
     WQ_changes_ACS_totals_table, WQ_changes_ACS_states_table)
  
}

# PART 8: Read in and preprocess tribal area data -----------------------------

US_tribal <- sf::st_read(paste0(GIS,"TIGER/2020/Tribal_Areas/tl_2020_us_aiannh/tl_2020_us_aiannh.shp"))
US_tribal_proj <- st_transform(US_tribal, 5072)

US_tribal_df <- st_drop_geometry(US_tribal_proj)
US_tribal_df$AREA <- with(US_tribal_df, ALAND+AWATER)

HAWQS_WBD_merge <- sf::st_read(paste0(GIS,"00_WOTUS/Analysis/2022/WOTUS_2022.gdb"),layer="HAWQS_WBD_Merge")

HAWQS_tribal <- st_join(HAWQS_WBD_merge, US_tribal_proj)

HAWQS_tribal_df <- st_drop_geometry(HAWQS_tribal)

  # QA - how many HUC12s does the tribal area overlap with?
  length(HAWQS_tribal_df$huc12[is.na(HAWQS_tribal$AIANNHCE)])
  length(HAWQS_tribal_df$huc12[!is.na(HAWQS_tribal$AIANNHCE)])
  
Tribal_Wet <- merge(HAWQS_tribal_df %>% dplyr::filter(!is.na(AIANNHCE)),
                    ORM_mitigation, by.x="huc12", by.y="HUC12", all.x=T)

Tribal_aff_area <- Tribal_Wet %>%
  dplyr::filter(!is.na(TOTAL_CHANGE)) %>%
  distinct(AIANNHCE, .keep_all = T) %>%
  mutate(AREA = ALAND+AWATER)

Perc_tribal <- sum(Tribal_aff_area$AREA)/sum(US_tribal_df$AREA)

rm(HAWQS_tribal, HAWQS_WBD_merge)
gc()

# PART 9: Inside and outside of tribal areas ----------------------------------

Tribal_ACS_env_table <- data.frame(c(0:11))
names(Tribal_ACS_env_table) <- "EXCEED_TOTAL"

# INSIDE
Tribal_in_ACS <- merge(HAWQS_tribal_df %>% dplyr::filter(!is.na(AIANNHCE)), 
                       BG_HUC12_Final, by.x="huc12", by.y="HUC12", all.x=T)
    
for (m in 1:length(env_var)) {
  
  Tribal_in_ACS[, paste0("EXCEED_",env_var[m])] <- 
    ifelse(Tribal_in_ACS[, paste0(env_var[m],"_num")] > 
             Natl_stats[, paste0(env_var[m],"_num")], 1, 0)
  
  jpeg(paste0(EJ,"R_programs/Outputs/Tribalin_Natl_Exceedances_",env_var[m],".jpg"))
  hist(Tribal_in_ACS[, paste0(env_var[m],"_num")])
  dev.off()
  
}

Tribal_in_ACS$PERC_LEAD_ENV <- with(Tribal_in_ACS, PRE1960/ACSTOTHU)

jpeg(paste0(EJ,"R_programs/Outputs/Tribalin_Natl_Exceedances_PERC_LEAD.jpg"))
hist(Tribal_in_ACS$PERC_LEAD_ENV)
dev.off()

Tribal_in_ACS$EXCEED_LEAD <- with(Tribal_in_ACS, 
                                      ifelse(PERC_LEAD_ENV > Natl_stats$PERC_LEAD, 1, 0))

# Make sure duplicate HUC12s are removed before summing.
Tribal_in_ACS_env <- Tribal_in_ACS %>%
  distinct(huc12, .keep_all=T) %>%
  mutate(EXCEED_TOTAL = rowSums(across(starts_with("EXCEED")), na.rm=T)) %>%
  group_by(EXCEED_TOTAL) %>%
  summarize(Exceedances_in = n_distinct(huc12))

Tribal_in_ACS_temp <- Tribal_in_ACS %>%
  distinct(huc12, .keep_all=T)

Tribal_in_ACS_totals <- data.frame(cbind(t(colSums(Tribal_in_ACS%>%select(TOT_AGE_POP:PRE1960), na.rm=T)),
                                              t(colMeans(Tribal_in_ACS%>%select(DSLPM_num:PM25_num), na.rm=T))))

Tribal_in_ACS_totals$PERC_MIN <- with(Tribal_in_ACS_totals, MINOR_POP/TOT_RACE_POP)
Tribal_in_ACS_totals$PERC_POV <- with(Tribal_in_ACS_totals, POV_POP/TOT_POV_POP)
Tribal_in_ACS_totals$PERC_LESSHS <- with(Tribal_in_ACS_totals, EDU_POP/TOT_EDU_POP)
Tribal_in_ACS_totals$PERC_LING <- with(Tribal_in_ACS_totals, LINGISO_HH/TOT_LINGISO_HH)
Tribal_in_ACS_totals$PERC_UNDER5 <- with(Tribal_in_ACS_totals, UNDER_5/TOT_AGE_POP)
Tribal_in_ACS_totals$PERC_OVER64 <- with(Tribal_in_ACS_totals, OVER_64/TOT_AGE_POP)
Tribal_in_ACS_totals$PERC_LEAD <- with(Tribal_in_ACS_totals, PRE1960/ACSTOTHU)
    
for (m in 1:length(race_var)) {
  
  Tribal_in_ACS_totals[, paste0("PERC_",race_var[m])] <- 
    Tribal_in_ACS_totals[, race_var[m]] / Tribal_in_ACS_totals$TOT_RACE_POP
  
}

Tribal_in_ACS_totals$HUC12_Num <- length(unique(Tribal_in_ACS$huc12))
Tribal_in_ACS_totals$HUC12_Perc <- Tribal_in_ACS_totals$HUC12_Num/length(unique(ORM_mitigation$HUC12))

Tribal_in_ACS_states <- Tribal_in_ACS %>%
  distinct(huc12, .keep_all=T)

Tribal_in_ACS_states$PERC_MIN_WET <- with(Tribal_in_ACS_states, MINOR_POP/TOT_RACE_POP)
Tribal_in_ACS_states$PERC_POV_WET <- with(Tribal_in_ACS_states, POV_POP/TOT_POV_POP)
Tribal_in_ACS_states$PERC_LESSHS_WET <- with(Tribal_in_ACS_states, EDU_POP/TOT_EDU_POP)
Tribal_in_ACS_states$PERC_LING_WET <- with(Tribal_in_ACS_states, LINGISO_HH/TOT_LINGISO_HH)
Tribal_in_ACS_states$PERC_UNDER5_WET <- with(Tribal_in_ACS_states, UNDER_5/TOT_AGE_POP)
Tribal_in_ACS_states$PERC_OVER64_WET <- with(Tribal_in_ACS_states, OVER_64/TOT_AGE_POP)
Tribal_in_ACS_states$PERC_LEAD_WET <- with(Tribal_in_ACS_states, PRE1960/ACSTOTHU)

for (m in 1:length(race_var)) {
  
  Tribal_in_ACS_states[, paste0("PERC_",race_var[m],"_WET")] <- 
    Tribal_in_ACS_states[, race_var[m]] / Tribal_in_ACS_states$TOT_RACE_POP
  
}

Tribal_in_ACS_states <- merge(Tribal_in_ACS_states, 
                                    State_stats %>% select(STATE_ID, PERC_MIN:PERC_HISPANIC), by="STATE_ID", all.x=T)

Tribal_in_ACS_states$EXCEED_MIN <- with(Tribal_in_ACS_states, ifelse(!is.na(PERC_MIN) & !is.na(PERC_MIN_WET) & 
                                                                                   PERC_MIN_WET>PERC_MIN, 1, 0))
Tribal_in_ACS_states$EXCEED_POV <- with(Tribal_in_ACS_states, ifelse(!is.na(PERC_POV) & !is.na(PERC_POV_WET) &
                                                                                   PERC_POV_WET>PERC_POV, 1, 0))
Tribal_in_ACS_states$EXCEED_LESSHS <- with(Tribal_in_ACS_states, ifelse(!is.na(PERC_LESSHS) & !is.na(PERC_LESSHS_WET) &
                                                                                      PERC_LESSHS_WET>PERC_LESSHS, 1, 0))
Tribal_in_ACS_states$EXCEED_LING <- with(Tribal_in_ACS_states, ifelse(!is.na(PERC_LING) & !is.na(PERC_LING_WET) &
                                                                                    PERC_LING_WET>PERC_LING, 1, 0))
Tribal_in_ACS_states$EXCEED_UNDER5 <- with(Tribal_in_ACS_states, ifelse(!is.na(PERC_UNDER5) & !is.na(PERC_UNDER5_WET) &
                                                                                      PERC_UNDER5_WET>PERC_UNDER5, 1, 0))
Tribal_in_ACS_states$EXCEED_OVER64 <- with(Tribal_in_ACS_states, ifelse(!is.na(PERC_OVER64) & !is.na(PERC_OVER64_WET) &
                                                                                      PERC_OVER64_WET>PERC_OVER64, 1, 0))
Tribal_in_ACS_states$EXCEED_LEAD <- with(Tribal_in_ACS_states, ifelse(!is.na(PERC_LEAD) & !is.na(PERC_LEAD_WET) &
                                                                                    PERC_LEAD_WET>PERC_LEAD, 1, 0))

for (m in 1:length(race_var)) {
  
  Tribal_in_ACS_states[, paste0("EXCEED_",race_var[m])] <- 
    ifelse(!is.na(Tribal_in_ACS_states[, paste0("PERC_",race_var[m],"_WET")]) & 
             !is.na(Tribal_in_ACS_states[, paste0("PERC_",race_var[m])]) & 
             Tribal_in_ACS_states[, paste0("PERC_",race_var[m],"_WET")] > 
             Tribal_in_ACS_states[, paste0("PERC_",race_var[m])], 1, 0)
  
}

Tribal_in_ACS_states_fin <- Tribal_in_ACS_states %>%
  select(starts_with("EXCEED")) %>%
  summarise_all(sum, na.rm=T)

Tribal_in_ACS_states_fin$HUC12_Num <- length(unique(Tribal_in_ACS$huc12))
Tribal_in_ACS_states_fin$HUC12_Perc <- Tribal_in_ACS_states_fin$HUC12_Num/length(unique(HAWQS_tribal_df$HUC12))

# OUTSIDE
Tribal_out_ACS <- merge(HAWQS_tribal_df %>% dplyr::filter(is.na(AIANNHCE)), 
                       BG_HUC12_Final, by.x="huc12", by.y="HUC12", all.x=T)

for (m in 1:length(env_var)) {
  
  Tribal_out_ACS[, paste0("EXCEED_",env_var[m])] <- 
    ifelse(Tribal_out_ACS[, paste0(env_var[m],"_num")] > 
             Natl_stats[, paste0(env_var[m],"_num")], 1, 0)
  
  jpeg(paste0(EJ,"R_programs/Outputs/Tribalout_Natl_Exceedances_",env_var[m],".jpg"))
  hist(Tribal_out_ACS[, paste0(env_var[m],"_num")])
  dev.off()
  
}

Tribal_out_ACS$PERC_LEAD_ENV <- with(Tribal_out_ACS, PRE1960/ACSTOTHU)

jpeg(paste0(EJ,"R_programs/Outputs/Tribalout_Natl_Exceedances_PERC_LEAD.jpg"))
hist(Tribal_out_ACS$PERC_LEAD_ENV)
dev.off()

Tribal_out_ACS$EXCEED_LEAD <- with(Tribal_out_ACS, 
                                  ifelse(PERC_LEAD_ENV > Natl_stats$PERC_LEAD, 1, 0))

# Make sure duplicate HUC12s are removed before summing.
Tribal_out_ACS_env <- Tribal_out_ACS %>%
  distinct(huc12, .keep_all=T) %>%
  mutate(EXCEED_TOTAL = rowSums(across(starts_with("EXCEED")), na.rm=T)) %>%
  group_by(EXCEED_TOTAL) %>%
  summarize(Exceedances_out = n_distinct(huc12))

Tribal_out_ACS_temp <- Tribal_out_ACS %>%
  distinct(huc12, .keep_all=T)

Tribal_out_ACS_totals <- data.frame(cbind(t(colSums(Tribal_out_ACS%>%select(TOT_AGE_POP:PRE1960), na.rm=T)),
                                         t(colMeans(Tribal_out_ACS%>%select(DSLPM_num:PM25_num), na.rm=T))))

Tribal_out_ACS_totals$PERC_MIN <- with(Tribal_out_ACS_totals, MINOR_POP/TOT_RACE_POP)
Tribal_out_ACS_totals$PERC_POV <- with(Tribal_out_ACS_totals, POV_POP/TOT_POV_POP)
Tribal_out_ACS_totals$PERC_LESSHS <- with(Tribal_out_ACS_totals, EDU_POP/TOT_EDU_POP)
Tribal_out_ACS_totals$PERC_LING <- with(Tribal_out_ACS_totals, LINGISO_HH/TOT_LINGISO_HH)
Tribal_out_ACS_totals$PERC_UNDER5 <- with(Tribal_out_ACS_totals, UNDER_5/TOT_AGE_POP)
Tribal_out_ACS_totals$PERC_OVER64 <- with(Tribal_out_ACS_totals, OVER_64/TOT_AGE_POP)
Tribal_out_ACS_totals$PERC_LEAD <- with(Tribal_out_ACS_totals, PRE1960/ACSTOTHU)

for (m in 1:length(race_var)) {
  
  Tribal_out_ACS_totals[, paste0("PERC_",race_var[m])] <- 
    Tribal_out_ACS_totals[, race_var[m]] / Tribal_out_ACS_totals$TOT_RACE_POP
  
}

Tribal_out_ACS_totals$HUC12_Num <- length(unique(Tribal_out_ACS$huc12))
Tribal_out_ACS_totals$HUC12_Perc <- Tribal_out_ACS_totals$HUC12_Num/length(unique(ORM_mitigation$HUC12))

Tribal_out_ACS_states <- Tribal_out_ACS %>%
  distinct(huc12, .keep_all=T)

Tribal_out_ACS_states$PERC_MIN_WET <- with(Tribal_out_ACS_states, MINOR_POP/TOT_RACE_POP)
Tribal_out_ACS_states$PERC_POV_WET <- with(Tribal_out_ACS_states, POV_POP/TOT_POV_POP)
Tribal_out_ACS_states$PERC_LESSHS_WET <- with(Tribal_out_ACS_states, EDU_POP/TOT_EDU_POP)
Tribal_out_ACS_states$PERC_LING_WET <- with(Tribal_out_ACS_states, LINGISO_HH/TOT_LINGISO_HH)
Tribal_out_ACS_states$PERC_UNDER5_WET <- with(Tribal_out_ACS_states, UNDER_5/TOT_AGE_POP)
Tribal_out_ACS_states$PERC_OVER64_WET <- with(Tribal_out_ACS_states, OVER_64/TOT_AGE_POP)
Tribal_out_ACS_states$PERC_LEAD_WET <- with(Tribal_out_ACS_states, PRE1960/ACSTOTHU)

for (m in 1:length(race_var)) {
  
  Tribal_out_ACS_states[, paste0("PERC_",race_var[m],"_WET")] <- 
    Tribal_out_ACS_states[, race_var[m]] / Tribal_out_ACS_states$TOT_RACE_POP
  
}

Tribal_out_ACS_states <- merge(Tribal_out_ACS_states, 
                              State_stats %>% select(STATE_ID, PERC_MIN:PERC_HISPANIC), by="STATE_ID", all.x=T)

Tribal_out_ACS_states$EXCEED_MIN <- with(Tribal_out_ACS_states, ifelse(!is.na(PERC_MIN) & !is.na(PERC_MIN_WET) & 
                                                                       PERC_MIN_WET>PERC_MIN, 1, 0))
Tribal_out_ACS_states$EXCEED_POV <- with(Tribal_out_ACS_states, ifelse(!is.na(PERC_POV) & !is.na(PERC_POV_WET) &
                                                                       PERC_POV_WET>PERC_POV, 1, 0))
Tribal_out_ACS_states$EXCEED_LESSHS <- with(Tribal_out_ACS_states, ifelse(!is.na(PERC_LESSHS) & !is.na(PERC_LESSHS_WET) &
                                                                          PERC_LESSHS_WET>PERC_LESSHS, 1, 0))
Tribal_out_ACS_states$EXCEED_LING <- with(Tribal_out_ACS_states, ifelse(!is.na(PERC_LING) & !is.na(PERC_LING_WET) &
                                                                        PERC_LING_WET>PERC_LING, 1, 0))
Tribal_out_ACS_states$EXCEED_UNDER5 <- with(Tribal_out_ACS_states, ifelse(!is.na(PERC_UNDER5) & !is.na(PERC_UNDER5_WET) &
                                                                          PERC_UNDER5_WET>PERC_UNDER5, 1, 0))
Tribal_out_ACS_states$EXCEED_OVER64 <- with(Tribal_out_ACS_states, ifelse(!is.na(PERC_OVER64) & !is.na(PERC_OVER64_WET) &
                                                                          PERC_OVER64_WET>PERC_OVER64, 1, 0))
Tribal_out_ACS_states$EXCEED_LEAD <- with(Tribal_out_ACS_states, ifelse(!is.na(PERC_LEAD) & !is.na(PERC_LEAD_WET) &
                                                                        PERC_LEAD_WET>PERC_LEAD, 1, 0))

for (m in 1:length(race_var)) {
  
  Tribal_out_ACS_states[, paste0("EXCEED_",race_var[m])] <- 
    ifelse(!is.na(Tribal_out_ACS_states[, paste0("PERC_",race_var[m],"_WET")]) & 
             !is.na(Tribal_out_ACS_states[, paste0("PERC_",race_var[m])]) & 
             Tribal_out_ACS_states[, paste0("PERC_",race_var[m],"_WET")] > 
             Tribal_out_ACS_states[, paste0("PERC_",race_var[m])], 1, 0)
  
}

Tribal_out_ACS_states_fin <- Tribal_out_ACS_states %>%
  select(starts_with("EXCEED")) %>%
  summarise_all(sum, na.rm=T)

Tribal_out_ACS_states_fin$HUC12_Num <- length(unique(Tribal_out_ACS$huc12))
Tribal_out_ACS_states_fin$HUC12_Perc <- Tribal_out_ACS_states_fin$HUC12_Num/length(unique(HAWQS_tribal_df$HUC12))

Tribal_ACS_env_table <- Reduce(function(x, y) merge(x=x, y=y, by="EXCEED_TOTAL", all.x=T),
                               list(Tribal_ACS_env_table, Tribal_in_ACS_env, Tribal_out_ACS_env))

# Cumulative impacts table
write.csv(Tribal_ACS_env_table,
          paste0(EJ,"R_programs/Outputs/Tribal_CumImpacts.csv"),row.names=FALSE)

Tribal_ACS_totals_table <- rbind(Tribal_in_ACS_totals, Tribal_out_ACS_totals)

# National comparison table
write.csv(Tribal_ACS_totals_table,
          paste0(EJ,"R_programs/Outputs/Tribal_Natl.csv"),row.names=FALSE)

Tribal_ACS_states_table <- rbind(Tribal_in_ACS_states_fin, Tribal_out_ACS_states_fin)

# State comparison table
write.csv(Tribal_ACS_states_table,
          paste0(EJ,"R_programs/Outputs/Tribal_States.csv"),row.names=FALSE)

rm(Tribal_ACS_env_table, Tribal_ACS_totals_table,
   Tribal_ACS_states_table)
