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

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

options(stringsAsFactors=FALSE)

library(dplyr)
library(readxl)

ORM <- "C:/Users/50526/ICF/WOTUS Reconsideration (ICF Internal Site) - ORM Data Processing/Outputs/"
EJ <- "C:/Users/50526/ICF/WOTUS Reconsideration - GIS - GIS/Environmental Justice/"

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

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

## ORM data
ORM_mitigation <- read.csv(paste0(ORM,"Impacts_By_HUC12_GIS.txt"),
                           colClasses=c("NULL","character",rep("NULL",18),"numeric"))

ORM_impactedwaters <- read.csv(paste0(ORM,"ImpactedWaters_By_HUC12_GIS.txt"),
                               colClasses=c("NULL","character",rep("NULL",9),"numeric"))

#########################################################################
# 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

# 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, 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 PWDIS variable 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",29),"character","NULL","numeric",
                                         rep("NULL",333)))
EJScreen_data$PWDIS_num <- as.numeric(EJScreen_data$PWDIS)
# 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, by.x="GEOID", by.y="ID", all=TRUE)

# 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

QA_BG_HUC12_EJ <- BG_HUC12_EJ %>% select(c("GEOID","BG_RATIO")) %>% 
  group_by(GEOID) %>% 
  summarise_all(sum)

#*# AL (10/6): The majority of BG_Ratios sum to 1 or nearly 1. Less than 1% of GEOIDs sum to less than 0.99 of those
#*#   BG_Ratios that are not "NA".

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)])

#*# AL (11/3): no longer need to calculate a HUC12_RATIO variable
# BG_HUC12_EJ$HUC12_RATIO <- with(BG_HUC12_EJ,BG_HUC12_SqKm/HUC12_SqKM)
# # QA: Confirmed that HUC12_RATIO values range from 0 to 1
# min(BG_HUC12_EJ$HUC12_RATIO, na.rm=T) # 0
# max(BG_HUC12_EJ$HUC12_RATIO, na.rm=T) # 1
# 
# QA_BG_HUC12_EJ2 <- BG_HUC12_EJ %>% select(c("HUC12","HUC12_RATIO")) %>%
#   group_by(HUC12) %>%
#   summarise_all(sum)
# 
# length(QA_BG_HUC12_EJ2$HUC12[QA_BG_HUC12_EJ2$HUC12_RATIO < 0.99 & !is.na(QA_BG_HUC12_EJ2$HUC12_RATIO)])/
#   length(QA_BG_HUC12_EJ2$HUC12[!is.na(QA_BG_HUC12_EJ2$HUC12_RATIO)])
# 
# #*# AL (10/6): The majority of HUC12_Ratios sum to 1 or nearly 1. Less than 1% of GEOIDs sum to less than 0.99 of those
# #*#   BG_Ratios that are not "NA".

BG_HUC12_adj <- BG_HUC12_EJ

# rm(BG_HUC12_EJ)

BG_HUC12_adj[, c("TOT_AGE_POP","UNDER_5","OVER_64", "TOT_RACE_POP", "MINOR_POP", "TOT_POV_POP",
                 "POV_POP", "TOT_EDU_POP", "EDU_POP", "TOT_LINGISO_HH", "LINGISO_HH")] <-
  BG_HUC12_adj[, c("TOT_AGE_POP","UNDER_5","OVER_64", "TOT_RACE_POP", "MINOR_POP", "TOT_POV_POP",
                   "POV_POP", "TOT_EDU_POP", "EDU_POP", "TOT_LINGISO_HH", "LINGISO_HH")] *
  BG_HUC12_adj[["BG_RATIO"]]

#*# KM: Spot check hand calcs confirm proper multiplication
#*# AL (09/29): There is only "NA" HUC12_RATIO's where we have EJScreen data and not ACS data

BG_HUC12_adj$PRMP_mult_BGHUC12Area <- BG_HUC12_adj$PRMP * BG_HUC12_adj$BG_HUC12_SqKm
BG_HUC12_adj$PWDIS_mult_BGHUC12Area <- BG_HUC12_adj$PWDIS_num * BG_HUC12_adj$BG_HUC12_SqKm

BG_HUC12_Final <- BG_HUC12_adj %>%
  select(HUC12, TOT_AGE_POP:LINGISO_HH, BG_HUC12_SqKm, PRMP_mult_BGHUC12Area:PWDIS_mult_BGHUC12Area) %>%
  group_by(HUC12) %>%
  summarise(across(TOT_AGE_POP:LINGISO_HH, ~ sum(.x, na.rm=T)),
            PRMP_adj = sum(PRMP_mult_BGHUC12Area)/sum(BG_HUC12_SqKm),
            PWDIS_adj = sum(PWDIS_mult_BGHUC12Area)/sum(BG_HUC12_SqKm),
            )

#*# KM (11/2): Resulting PRMP_adj value for 031502010203 is 0.08649409, 
#*# which is within range of the original PRMP values associated with this HUC12
#*# AL (11/3): Looked at values for HUC 150501000809. Original PRMP values range from 0.05 to 0.32
#*#    Resulting PRMP_adj value is 0.1027613 which is within range of the original PRMP values.

#*# AL (11/1) - no longer need to calculate demographic information at the CBG level
# BG_HUC12_Final$PERC_MIN <- with(BG_HUC12_Final, ifelse(TOT_RACE_POP==0,0,MINOR_POP/TOT_RACE_POP*100))
# BG_HUC12_Final$PERC_POV <- with(BG_HUC12_Final, ifelse(TOT_POV_POP==0,0,POV_POP/TOT_POV_POP*100))
# BG_HUC12_Final$DEM_IND <- with(BG_HUC12_Final, (PERC_MIN+PERC_POV)/2)
# BG_HUC12_Final$PERC_LESSHS <- with(BG_HUC12_Final, ifelse(TOT_EDU_POP==0,0,EDU_POP/TOT_EDU_POP*100))
# BG_HUC12_Final$PERC_LING <- with(BG_HUC12_Final, ifelse(TOT_LINGISO_HH==0,0,LINGISO_HH/TOT_LINGISO_HH*100))
# BG_HUC12_Final$PERC_UNDER5 <- with(BG_HUC12_Final, ifelse(TOT_AGE_POP==0,0,UNDER_5/TOT_AGE_POP*100))
# BG_HUC12_Final$PERC_OVER64 <- with(BG_HUC12_Final, ifelse(TOT_AGE_POP==0,0,OVER_64/TOT_AGE_POP*100))

# Calculate national statistics
HUC12_natl <- data.frame(t(colSums(BG_HUC12_Final%>%select(TOT_AGE_POP:LINGISO_HH),
                                   na.rm=T)))
HUC12_natl$PERC_MIN <- with(HUC12_natl, MINOR_POP/TOT_RACE_POP*100)
HUC12_natl$PERC_POV <- with(HUC12_natl, POV_POP/TOT_POV_POP*100)
#*# AL (11/1) no longer using demographic index
# HUC12_natl$DEM_IND <- with(HUC12_natl, (PERC_MIN+PERC_POV)/2)
HUC12_natl$PERC_LESSHS <- with(HUC12_natl, EDU_POP/TOT_EDU_POP*100)
HUC12_natl$PERC_LING <- with(HUC12_natl, LINGISO_HH/TOT_LINGISO_HH*100)
HUC12_natl$PERC_UNDER5 <- with(HUC12_natl, UNDER_5/TOT_AGE_POP*100)
HUC12_natl$PERC_OVER64 <- with(HUC12_natl, OVER_64/TOT_AGE_POP*100)

HUC12_natl$PWDIS <- mean(BG_HUC12_Final$PWDIS_adj, na.rm=T)
HUC12_natl$PRMP <- mean(BG_HUC12_Final$PRMP_adj, na.r=T)

HUC12_natl_EJ <- data.frame(t(colSums(BG_HUC12_Final[which(!is.na(BG_HUC12_Final$PRMP_adj) & 
                                                        BG_HUC12_Final$PRMP_adj>0 & 
                                                        !is.na(BG_HUC12_Final$PWDIS_adj) & 
                                                        BG_HUC12_Final$PWDIS_adj>0), ]%>%
                                     select(TOT_AGE_POP:LINGISO_HH), na.rm=T)))
HUC12_natl_EJ$PERC_MIN <- with(HUC12_natl_EJ, MINOR_POP/TOT_RACE_POP*100)
HUC12_natl_EJ$PERC_POV <- with(HUC12_natl_EJ, POV_POP/TOT_POV_POP*100)
HUC12_natl_EJ$PERC_LESSHS <- with(HUC12_natl_EJ, EDU_POP/TOT_EDU_POP*100)
HUC12_natl_EJ$PERC_LING <- with(HUC12_natl_EJ, LINGISO_HH/TOT_LINGISO_HH*100)
HUC12_natl_EJ$PERC_UNDER5 <- with(HUC12_natl_EJ, UNDER_5/TOT_AGE_POP*100)
HUC12_natl_EJ$PERC_OVER64 <- with(HUC12_natl_EJ, OVER_64/TOT_AGE_POP*100)

HUC12_natl_EJ$PWDIS <- mean(BG_HUC12_Final$PWDIS_adj[which(!is.na(BG_HUC12_Final$PRMP_adj) & 
                                                             BG_HUC12_Final$PRMP_adj>0 & 
                                                             !is.na(BG_HUC12_Final$PWDIS_adj) & 
                                                             BG_HUC12_Final$PWDIS_adj>0)], na.rm=T)
HUC12_natl_EJ$PRMP <- mean(BG_HUC12_Final$PRMP_adj[which(!is.na(BG_HUC12_Final$PRMP_adj) & 
                                                           BG_HUC12_Final$PRMP_adj>0 & 
                                                           !is.na(BG_HUC12_Final$PWDIS_adj) & 
                                                           BG_HUC12_Final$PWDIS_adj>0)], na.r=T)

write.csv(HUC12_natl_EJ, paste0(EJ,"R_programs/Outputs/Natl_Stats_EJ.csv"),
                                row.names=F)

#*# KM (11/2): Confirmed that values match prior results, where appropriate

#*# AL (11/1) - not basing the tables off of percentile exceedances
# # Calculate percentiles for each HUC12 (80th, 90th, 95th percentiles)
# # Using the 80th percentile for this analysis
# BG_P80 <- apply(BG_HUC12_Final%>%select(PERC_MIN:PERC_OVER64), 2,
#                 function(x) quantile(x, probs=.8))
# BG_P90 <- apply(BG_HUC12_Final%>%select(PERC_MIN:PERC_OVER64), 2,
#                 function(x) quantile(x, probs=.9))
# BG_P95 <- apply(BG_HUC12_Final%>%select(PERC_MIN:PERC_OVER64), 2,
#                 function(x) quantile(x, probs=.95))
# 
# # Identify HUC12s that exceed the 80th percentile
# BG_HUC12_Final$P80_MIN <- with(BG_HUC12_Final, ifelse(PERC_MIN>BG_P80["PERC_MIN"],1,0))
# BG_HUC12_Final$P80_POV <- with(BG_HUC12_Final, ifelse(PERC_POV>BG_P80["PERC_POV"],1,0))
# BG_HUC12_Final$P80_DEM_IND <- with(BG_HUC12_Final, ifelse(DEM_IND>BG_P80["DEM_IND"],1,0))
# BG_HUC12_Final$P80_LESSHS <- with(BG_HUC12_Final, ifelse(PERC_LESSHS>BG_P80["PERC_LESSHS"],1,0))
# BG_HUC12_Final$P80_LING <- with(BG_HUC12_Final, ifelse(PERC_LING>BG_P80["PERC_LING"],1,0))
# BG_HUC12_Final$P80_UNDER5 <- with(BG_HUC12_Final, ifelse(PERC_UNDER5>BG_P80["PERC_UNDER5"],1,0))
# BG_HUC12_Final$P80_OVER64 <- with(BG_HUC12_Final, ifelse(PERC_OVER64>BG_P80["PERC_OVER64"],1,0))
# 
# BG_HUC12_Final$P80_DEM_RSEI <- with(BG_HUC12_Final, ifelse(DEM_IND>BG_P80["DEM_IND"]&P_RSEI_adj>80,1,0))
# BG_HUC12_Final$P80_DEM_RMP <- with(BG_HUC12_Final, ifelse(DEM_IND>BG_P80["DEM_IND"]&P_RMP_adj>80,1,0))
# BG_HUC12_Final$P80_DEM_RSEI_RMP <- with(BG_HUC12_Final, ifelse(DEM_IND>BG_P80["DEM_IND"]&P_RSEI_adj>80&P_RMP_adj>80,
#                                                                1,0))

# Census tract level
Tract_HUC12_Intersect <- read.csv(paste0(EJ,"GIS_Outputs/AllStatesCT_HUC12_TableToFile.csv"),
                                  colClasses = c("NULL","character","integer","numeric","integer",
                                                 rep("character",2),rep("numeric",2)))
sort(unique(substr(Tract_HUC12_Intersect$huc12_num,1,2)))

# HUC12s missing the leading zero
Tract_HUC12_Intersect$huc12 <- with(Tract_HUC12_Intersect, ifelse(nchar(huc12_num) < 12, 
                                                                  paste0("0", huc12_num), huc12_num))
sort(unique(substr(Tract_HUC12_Intersect$huc12,1,2)))

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

#*# KM (09/30): Any idea why some CTs (e.g., 56039967704) are in the census tract file (Tract_HUC12_Intersect), 
#*# but not the CBG file (BG_HUC12_Intersect; has 56039967701X but not 56039967704X and 56039967703X)?
#*# AL (10/6): The two still don't exist in both. However, in the CT file, it doesn't have any population information.

# Aggregate the total race population counts from the block group to census tract level
# B02014_001 variable is a count of people who are American Indian and Alaska Native alone and people with no tribe
# which is a subset of the people sampled for race.
BG_HUC12_Var$GEOID_Tract <- substr(BG_HUC12_Var$GEOID,1,11)
Tract_TOT_RACE <- BG_HUC12_Var %>%
  select(GEOID_Tract, TOT_RACE_POP) %>%
  group_by(GEOID_Tract) %>%
  summarise(TOT_RACE_POP_CT=sum(TOT_RACE_POP,na.rm=T))

Tract_HUC12_Race <- merge(Tract_HUC12_Intersect, Tract_TOT_RACE, by.x="GEOID", by.y="GEOID_Tract",
                          all.x=T)

Tract_HUC12_Var <- Tract_HUC12_Race %>%
  select(HUC12=huc12, GEOID, TOT_RACE_POP_CT, TRIBE_POP=B02014_002, CT_SqKM=CT_Full_SqKM,
         CT_HUC12_SqKm)

# Adjust variables by CT/HUC12 overlap
Tract_HUC12_Var$CT_RATIO <- with(Tract_HUC12_Var,CT_HUC12_SqKm/CT_SqKM)
# QA: Confirmed that CT_RATIO values range from 0 to 1
min(Tract_HUC12_Var$CT_RATIO, na.rm=T) # 0
max(Tract_HUC12_Var$CT_RATIO, na.rm=T) # 1

QA_Tract_HUC12_Var <- Tract_HUC12_Var %>% select(c("GEOID","CT_RATIO")) %>% 
  group_by(GEOID) %>% 
  summarise_all(sum)

length(QA_Tract_HUC12_Var$GEOID[QA_Tract_HUC12_Var$CT_RATIO<0.99 & !is.na(QA_Tract_HUC12_Var$CT_RATIO)])/
  length(QA_Tract_HUC12_Var$GEOID[!is.na(QA_Tract_HUC12_Var$CT_RATIO)])

#*# AL (10/6): The majority of census tract ratios add up to nearly 1. Less than 2% of census tracts have ratios that add up
#*#   to less than 1 and are not "NA"

CT_HUC12_adj <- Tract_HUC12_Var

CT_HUC12_adj[, c("TOT_RACE_POP_CT","TRIBE_POP")] <-
  CT_HUC12_adj[, c("TOT_RACE_POP_CT","TRIBE_POP")] * CT_HUC12_adj[["CT_RATIO"]]

CT_HUC12_Final <- CT_HUC12_adj %>%
  select(HUC12, TOT_RACE_POP_CT, TRIBE_POP) %>%
  group_by(HUC12) %>%
  summarise(across(TOT_RACE_POP_CT:TRIBE_POP, ~ sum(.x, na.rm=T)))

CT_HUC12_Final$PERC_TRIBAL <- with(CT_HUC12_Final, ifelse(TOT_RACE_POP_CT==0|is.na(TOT_RACE_POP_CT),0,
                                                          TRIBE_POP/TOT_RACE_POP_CT*100))

# Calculate national statistics
HUC12_CT_natl <- data.frame(t(colSums(CT_HUC12_Final%>%select(TOT_RACE_POP_CT:TRIBE_POP),
                                      na.rm=T)))
HUC12_CT_natl$PERC_TRIBAL <- with(HUC12_CT_natl, TRIBE_POP/TOT_RACE_POP_CT*100)

HUC12_natl_all <- cbind(HUC12_natl, HUC12_CT_natl)

write.csv(HUC12_natl_all,paste0(EJ,"R_programs/Outputs/Natl_stats.csv"),row.names=F)

#*# AL (11/1) - not basing the tables off of percentile exceedances
# # Calculate percentiles for each HUC12 (80th, 90th, 95th percentiles)
# # Using the 80th percentile for this analysis
# CT_P80 <- apply(CT_HUC12_Final%>%select(PERC_TRIBAL), 2,
#                 function(x) quantile(x, probs=.8))
# CT_P90 <- apply(CT_HUC12_Final%>%select(PERC_TRIBAL), 2,
#                 function(x) quantile(x, probs=.9))
# CT_P95 <- apply(CT_HUC12_Final%>%select(PERC_TRIBAL), 2,
#                 function(x) quantile(x, probs=.95))
# 
# # Identify HUC12s that exceed the 80th percentile
# CT_HUC12_Final$P80_TRIBAL <- with(CT_HUC12_Final, ifelse(PERC_TRIBAL>CT_P80["PERC_TRIBAL"],1,0))

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

## Wetland Area Changes
ORM_mitig_ACS <- Reduce(function(x,y) merge(x=x, y=y, by="HUC12", all=T),
                        list(ORM_mitigation, BG_HUC12_Final, CT_HUC12_Final))

# rm(ORM_mitigation)

# QA
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"), ]

# write.csv(ORM_mitig_BG_QA, paste0(EJ,"R_programs/Outputs/QA/BG_Intersect_QA.csv"))

## Impacted Waters Changes
ORM_IW_ACS <-  Reduce(function(x,y) merge(x=x, y=y, by="HUC12", all=T),
                      list(ORM_impactedwaters, BG_HUC12_Final, CT_HUC12_Final))

# 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"), ]

# write.csv(ORM_IW_BG_QA, paste0(EJ,"R_programs/Outputs/QA/IW_BG_Intersect_QA.csv"))

#########################################################################
# PART 4: Bin the outputs based on wetland area/impacted waters changes
#########################################################################

## Wetland Area Changes
Mitig_ACS_quant <- list()
Mitig_ACS_totals <- list()

Mitig_ACS_PWDIS_quant <- list()
Mitig_ACS_PWDIS_totals <- list()

Mitig_ACS_PRMP_quant <- list()
Mitig_ACS_PRMP_totals <- list()

Mitig_ACS_EJ_quant <- list()
Mitig_ACS_EJ_totals <- list()

for (i in 1:5) {
  
  # Splitting up the data into the bins agreed upon with EPA on 09/30
  if (i==1) {
    Mitig_ACS_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                  ORM_mitig_ACS$TOTAL_CHANGE<=wet_area_break[i]), ]
    
    Mitig_ACS_PWDIS_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                        ORM_mitig_ACS$TOTAL_CHANGE<=wet_area_break[i] & 
                                                        !is.na(ORM_mitig_ACS$PWDIS_adj) & 
                                                        ORM_mitig_ACS$PWDIS_adj>0), ]
    
    Mitig_ACS_PRMP_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                       ORM_mitig_ACS$TOTAL_CHANGE<=wet_area_break[i] & 
                                                       !is.na(ORM_mitig_ACS$PRMP_adj) & 
                                                       ORM_mitig_ACS$PRMP_adj>0), ]
    
    Mitig_ACS_EJ_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                     ORM_mitig_ACS$TOTAL_CHANGE<=wet_area_break[i] & 
                                                     !is.na(ORM_mitig_ACS$PRMP_adj) & 
                                                     ORM_mitig_ACS$PRMP_adj>0 & 
                                                     !is.na(ORM_mitig_ACS$PWDIS_adj) & 
                                                     ORM_mitig_ACS$PWDIS_adj>0), ]
    
  } else if (i<4) {
    Mitig_ACS_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                  wet_area_break[i-1]<ORM_mitig_ACS$TOTAL_CHANGE&
                                                  ORM_mitig_ACS$TOTAL_CHANGE<=wet_area_break[i]), ]
    
    Mitig_ACS_PWDIS_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                        wet_area_break[i-1]<ORM_mitig_ACS$TOTAL_CHANGE&
                                                        ORM_mitig_ACS$TOTAL_CHANGE<=wet_area_break[i] & 
                                                        !is.na(ORM_mitig_ACS$PWDIS_adj) & 
                                                        ORM_mitig_ACS$PWDIS_adj>0), ]
    
    Mitig_ACS_PRMP_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                       wet_area_break[i-1]<ORM_mitig_ACS$TOTAL_CHANGE&
                                                       ORM_mitig_ACS$TOTAL_CHANGE<=wet_area_break[i] & 
                                                       !is.na(ORM_mitig_ACS$PRMP_adj) & 
                                                       ORM_mitig_ACS$PRMP_adj>0), ]
    
    Mitig_ACS_EJ_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                     wet_area_break[i-1]<ORM_mitig_ACS$TOTAL_CHANGE&
                                                     ORM_mitig_ACS$TOTAL_CHANGE<=wet_area_break[i] & 
                                                     !is.na(ORM_mitig_ACS$PRMP_adj) & 
                                                     ORM_mitig_ACS$PRMP_adj>0 & 
                                                     !is.na(ORM_mitig_ACS$PWDIS_adj) & 
                                                     ORM_mitig_ACS$PWDIS_adj>0), ]
    
  } else if (i==4) {
    Mitig_ACS_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                  ORM_mitig_ACS$TOTAL_CHANGE>wet_area_break[i-1]), ]
    
    Mitig_ACS_PWDIS_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                        ORM_mitig_ACS$TOTAL_CHANGE>wet_area_break[i-1] & 
                                                        !is.na(ORM_mitig_ACS$PWDIS_adj) & 
                                                        ORM_mitig_ACS$PWDIS_adj>0), ]
    
    Mitig_ACS_PRMP_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                       ORM_mitig_ACS$TOTAL_CHANGE>wet_area_break[i-1] & 
                                                       !is.na(ORM_mitig_ACS$PRMP_adj) & 
                                                       ORM_mitig_ACS$PRMP_adj>0), ]
    
    Mitig_ACS_EJ_quant[[i]] <- ORM_mitig_ACS[which(!is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                     ORM_mitig_ACS$TOTAL_CHANGE>wet_area_break[i-1] & 
                                                     !is.na(ORM_mitig_ACS$PRMP_adj) & 
                                                     ORM_mitig_ACS$PRMP_adj>0 & 
                                                     !is.na(ORM_mitig_ACS$PWDIS_adj) & 
                                                     ORM_mitig_ACS$PWDIS_adj), ]
    
  } else if (i==5) {
    Mitig_ACS_quant[[i]] <- ORM_mitig_ACS[which(is.na(ORM_mitig_ACS$TOTAL_CHANGE)), ]
    
    Mitig_ACS_PWDIS_quant[[i]] <- ORM_mitig_ACS[which(is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                        !is.na(ORM_mitig_ACS$PWDIS_adj) & 
                                                        ORM_mitig_ACS$PWDIS_adj>0), ]
    
    Mitig_ACS_PRMP_quant[[i]] <- ORM_mitig_ACS[which(is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                       !is.na(ORM_mitig_ACS$PRMP_adj) & 
                                                       ORM_mitig_ACS$PRMP_adj>0), ]
    
    Mitig_ACS_EJ_quant[[i]] <- ORM_mitig_ACS[which(is.na(ORM_mitig_ACS$TOTAL_CHANGE) & 
                                                     !is.na(ORM_mitig_ACS$PRMP_adj) & 
                                                     ORM_mitig_ACS$PRMP_adj>0 & 
                                                     !is.na(ORM_mitig_ACS$PWDIS_adj) & 
                                                     ORM_mitig_ACS$PWDIS_adj>0), ]
    
  }
  
  #*# AL (11/1) - no longer based on percentiles
  # Mitig_ACS_totals[[i]] <- data.frame(t(colSums(Mitig_ACS_quant[[i]]%>%select(P80_MIN:P80_DEM_RSEI_RMP,P80_TRIBAL), na.rm=TRUE)))
  Mitig_ACS_totals[[i]] <- data.frame(t(colSums(Mitig_ACS_quant[[i]]%>%select(TOT_AGE_POP:LINGISO_HH,TOT_RACE_POP_CT:TRIBE_POP),
                                                na.rm=T)))
  Mitig_ACS_totals[[i]]$PERC_MIN <- with(Mitig_ACS_totals[[i]], MINOR_POP/TOT_RACE_POP*100)
  Mitig_ACS_totals[[i]]$PERC_POV <- with(Mitig_ACS_totals[[i]], POV_POP/TOT_POV_POP*100)
  Mitig_ACS_totals[[i]]$PERC_LESSHS <- with(Mitig_ACS_totals[[i]], EDU_POP/TOT_EDU_POP*100)
  Mitig_ACS_totals[[i]]$PERC_LING <- with(Mitig_ACS_totals[[i]], LINGISO_HH/TOT_LINGISO_HH*100)
  Mitig_ACS_totals[[i]]$PERC_UNDER5 <- with(Mitig_ACS_totals[[i]], UNDER_5/TOT_AGE_POP*100)
  Mitig_ACS_totals[[i]]$PERC_OVER64 <- with(Mitig_ACS_totals[[i]], OVER_64/TOT_AGE_POP*100)
  Mitig_ACS_totals[[i]]$PERC_TRIBAL <- with(Mitig_ACS_totals[[i]], TRIBE_POP/TOT_RACE_POP_CT*100)
  
  Mitig_ACS_totals[[i]]$PWDIS <- mean(Mitig_ACS_quant[[i]]$PWDIS_adj, na.rm=T)
  Mitig_ACS_totals[[i]]$PRMP <- mean(Mitig_ACS_quant[[i]]$PRMP_adj, na.r=T)
  
  Mitig_ACS_totals[[i]]$HUC12_Num <- length(unique(Mitig_ACS_quant[[i]]$HUC12))
  Mitig_ACS_totals[[i]]$ORM_Total <- length(unique(ORM_mitigation$HUC12))
  
  write.csv(Mitig_ACS_totals[[i]],paste0(EJ,"R_programs/Outputs/WetArea_",i,"_ACS.csv"),row.names=FALSE)
  
  Mitig_ACS_PWDIS_totals[[i]] <- data.frame(t(colSums(Mitig_ACS_PWDIS_quant[[i]]%>%select(TOT_AGE_POP:LINGISO_HH,TOT_RACE_POP_CT:TRIBE_POP),
                                                      na.rm=T)))
  Mitig_ACS_PWDIS_totals[[i]]$PERC_MIN <- with(Mitig_ACS_PWDIS_totals[[i]], MINOR_POP/TOT_RACE_POP*100)
  Mitig_ACS_PWDIS_totals[[i]]$PERC_POV <- with(Mitig_ACS_PWDIS_totals[[i]], POV_POP/TOT_POV_POP*100)
  Mitig_ACS_PWDIS_totals[[i]]$PERC_LESSHS <- with(Mitig_ACS_PWDIS_totals[[i]], EDU_POP/TOT_EDU_POP*100)
  Mitig_ACS_PWDIS_totals[[i]]$PERC_LING <- with(Mitig_ACS_PWDIS_totals[[i]], LINGISO_HH/TOT_LINGISO_HH*100)
  Mitig_ACS_PWDIS_totals[[i]]$PERC_UNDER5 <- with(Mitig_ACS_PWDIS_totals[[i]], UNDER_5/TOT_AGE_POP*100)
  Mitig_ACS_PWDIS_totals[[i]]$PERC_OVER64 <- with(Mitig_ACS_PWDIS_totals[[i]], OVER_64/TOT_AGE_POP*100)
  Mitig_ACS_PWDIS_totals[[i]]$PERC_TRIBAL <- with(Mitig_ACS_PWDIS_totals[[i]], TRIBE_POP/TOT_RACE_POP_CT*100)
  
  Mitig_ACS_PWDIS_totals[[i]]$PWDIS <- mean(Mitig_ACS_PWDIS_quant[[i]]$PWDIS_adj, na.rm=T)
  Mitig_ACS_PWDIS_totals[[i]]$PRMP <- mean(Mitig_ACS_PWDIS_quant[[i]]$PRMP_adj, na.r=T)
  
  Mitig_ACS_PWDIS_totals[[i]]$HUC12_Num <- length(unique(Mitig_ACS_PWDIS_quant[[i]]$HUC12))
  Mitig_ACS_PWDIS_totals[[i]]$ORM_Total <- length(unique(ORM_mitigation$HUC12))
  
  write.csv(Mitig_ACS_PWDIS_totals[[i]],paste0(EJ,"R_programs/Outputs/WetArea_PWDIS_",i,"_ACS.csv"),row.names=FALSE)
  
  Mitig_ACS_PRMP_totals[[i]] <- data.frame(t(colSums(Mitig_ACS_PRMP_quant[[i]]%>%select(TOT_AGE_POP:LINGISO_HH,TOT_RACE_POP_CT:TRIBE_POP),
                                                     na.rm=T)))
  Mitig_ACS_PRMP_totals[[i]]$PERC_MIN <- with(Mitig_ACS_PRMP_totals[[i]], MINOR_POP/TOT_RACE_POP*100)
  Mitig_ACS_PRMP_totals[[i]]$PERC_POV <- with(Mitig_ACS_PRMP_totals[[i]], POV_POP/TOT_POV_POP*100)
  Mitig_ACS_PRMP_totals[[i]]$PERC_LESSHS <- with(Mitig_ACS_PRMP_totals[[i]], EDU_POP/TOT_EDU_POP*100)
  Mitig_ACS_PRMP_totals[[i]]$PERC_LING <- with(Mitig_ACS_PRMP_totals[[i]], LINGISO_HH/TOT_LINGISO_HH*100)
  Mitig_ACS_PRMP_totals[[i]]$PERC_UNDER5 <- with(Mitig_ACS_PRMP_totals[[i]], UNDER_5/TOT_AGE_POP*100)
  Mitig_ACS_PRMP_totals[[i]]$PERC_OVER64 <- with(Mitig_ACS_PRMP_totals[[i]], OVER_64/TOT_AGE_POP*100)
  Mitig_ACS_PRMP_totals[[i]]$PERC_TRIBAL <- with(Mitig_ACS_PRMP_totals[[i]], TRIBE_POP/TOT_RACE_POP_CT*100)
  
  Mitig_ACS_PRMP_totals[[i]]$PWDIS <- mean(Mitig_ACS_PRMP_quant[[i]]$PWDIS_adj, na.rm=T)
  Mitig_ACS_PRMP_totals[[i]]$PRMP <- mean(Mitig_ACS_PRMP_quant[[i]]$PRMP_adj, na.r=T)
  
  Mitig_ACS_PRMP_totals[[i]]$HUC12_Num <- length(unique(Mitig_ACS_PRMP_quant[[i]]$HUC12))
  Mitig_ACS_PRMP_totals[[i]]$ORM_Total <- length(unique(ORM_mitigation$HUC12))
  
  write.csv(Mitig_ACS_PRMP_totals[[i]],paste0(EJ,"R_programs/Outputs/WetArea_PRMP_",i,"_ACS.csv"),row.names=FALSE)
  
  Mitig_ACS_EJ_totals[[i]] <- data.frame(t(colSums(Mitig_ACS_EJ_quant[[i]]%>%select(TOT_AGE_POP:LINGISO_HH,TOT_RACE_POP_CT:TRIBE_POP),
                                                   na.rm=T)))
  Mitig_ACS_EJ_totals[[i]]$PERC_MIN <- with(Mitig_ACS_EJ_totals[[i]], MINOR_POP/TOT_RACE_POP*100)
  Mitig_ACS_EJ_totals[[i]]$PERC_POV <- with(Mitig_ACS_EJ_totals[[i]], POV_POP/TOT_POV_POP*100)
  Mitig_ACS_EJ_totals[[i]]$PERC_LESSHS <- with(Mitig_ACS_EJ_totals[[i]], EDU_POP/TOT_EDU_POP*100)
  Mitig_ACS_EJ_totals[[i]]$PERC_LING <- with(Mitig_ACS_EJ_totals[[i]], LINGISO_HH/TOT_LINGISO_HH*100)
  Mitig_ACS_EJ_totals[[i]]$PERC_UNDER5 <- with(Mitig_ACS_EJ_totals[[i]], UNDER_5/TOT_AGE_POP*100)
  Mitig_ACS_EJ_totals[[i]]$PERC_OVER64 <- with(Mitig_ACS_EJ_totals[[i]], OVER_64/TOT_AGE_POP*100)
  Mitig_ACS_EJ_totals[[i]]$PERC_TRIBAL <- with(Mitig_ACS_EJ_totals[[i]], TRIBE_POP/TOT_RACE_POP_CT*100)
  
  Mitig_ACS_EJ_totals[[i]]$PWDIS <- mean(Mitig_ACS_EJ_quant[[i]]$PWDIS_adj, na.rm=T)
  Mitig_ACS_EJ_totals[[i]]$PRMP <- mean(Mitig_ACS_EJ_quant[[i]]$PRMP_adj, na.r=T)
  
  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))
  
  write.csv(Mitig_ACS_EJ_totals[[i]],paste0(EJ,"R_programs/Outputs/WetArea_EJ_",i,"_ACS.csv"),row.names=FALSE)
  
}

#*# KM: Confirmed that sum of all HUC12s in Mitig_ACS_totals[[1:5]] = length(unique(ORM_mitig_ACS$HUC12))

## Impacted waters
IW_ACS_quant <- list()
IW_ACS_totals <- list()

IW_ACS_PWDIS_quant <- list()
IW_ACS_PWDIS_totals <- list()

IW_ACS_PRMP_quant <- list()
IW_ACS_PRMP_totals <- list()

IW_ACS_EJ_quant <- list()
IW_ACS_EJ_totals <- list()

for (i in 1:5) {
  
  # Splitting up the data into the first quantile, IQR, and last quantile
  if (i==1) {
    IW_ACS_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                            ORM_IW_ACS$WATERS_CHANGE<=water_break[i]), ]
    
    IW_ACS_PWDIS_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                                  ORM_IW_ACS$WATERS_CHANGE<=water_break[i] & 
                                                  !is.na(ORM_IW_ACS$PWDIS_adj) & 
                                                  ORM_IW_ACS$PWDIS_adj>0), ]
    
    IW_ACS_PRMP_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                                 ORM_IW_ACS$WATERS_CHANGE<=water_break[i] & 
                                                 !is.na(ORM_IW_ACS$PRMP_adj) & 
                                                 ORM_IW_ACS$PRMP_adj>0), ]
    
    IW_ACS_EJ_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                               ORM_IW_ACS$WATERS_CHANGE<=water_break[i] & 
                                               !is.na(ORM_IW_ACS$PWDIS_adj) & 
                                               ORM_IW_ACS$PWDIS_adj>0 & 
                                               !is.na(ORM_IW_ACS$PRMP_adj) & 
                                               ORM_IW_ACS$PRMP_adj>0), ]
    
  } else if (i<4) {
    IW_ACS_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                            water_break[i-1]<ORM_IW_ACS$WATERS_CHANGE&
                                            ORM_IW_ACS$WATERS_CHANGE<=water_break[i]), ]
    
    IW_ACS_PWDIS_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                                  water_break[i-1]<ORM_IW_ACS$WATERS_CHANGE&
                                                  ORM_IW_ACS$WATERS_CHANGE<=water_break[i] & 
                                                  !is.na(ORM_IW_ACS$PWDIS_adj) & 
                                                  ORM_IW_ACS$PWDIS_adj>0), ]
    
    IW_ACS_PRMP_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                                 water_break[i-1]<ORM_IW_ACS$WATERS_CHANGE&
                                                 ORM_IW_ACS$WATERS_CHANGE<=water_break[i] & 
                                                 !is.na(ORM_IW_ACS$PRMP_adj) & 
                                                 ORM_IW_ACS$PRMP_adj>0), ]
    
    IW_ACS_EJ_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                               water_break[i-1]<ORM_IW_ACS$WATERS_CHANGE&
                                               ORM_IW_ACS$WATERS_CHANGE<=water_break[i] & 
                                               !is.na(ORM_IW_ACS$PWDIS_adj) & 
                                               ORM_IW_ACS$PWDIS_adj>0 & 
                                               !is.na(ORM_IW_ACS$PRMP_adj) & 
                                               ORM_IW_ACS$PRMP_adj>0), ]
    
  } else if (i==4) {
    IW_ACS_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                            ORM_IW_ACS$WATERS_CHANGE>water_break[i-1]), ]
    
    IW_ACS_PWDIS_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                                  ORM_IW_ACS$WATERS_CHANGE>water_break[i-1] & 
                                                  !is.na(ORM_IW_ACS$PWDIS_adj) & 
                                                  ORM_IW_ACS$PWDIS_adj>0), ]
    
    IW_ACS_PRMP_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                                 ORM_IW_ACS$WATERS_CHANGE>water_break[i-1] & 
                                                 !is.na(ORM_IW_ACS$PRMP_adj) & 
                                                 ORM_IW_ACS$PRMP_adj>0), ]
    
    IW_ACS_EJ_quant[[i]] <- ORM_IW_ACS[which(!is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                               ORM_IW_ACS$WATERS_CHANGE>water_break[i-1] & 
                                               !is.na(ORM_IW_ACS$PWDIS_adj) & 
                                               ORM_IW_ACS$PWDIS_adj>0 & 
                                               !is.na(ORM_IW_ACS$PRMP_adj) & 
                                               ORM_IW_ACS$PRMP_adj>0), ]
    
  } else if (i==5) {
    IW_ACS_quant[[i]] <- ORM_IW_ACS[which(is.na(ORM_IW_ACS$WATERS_CHANGE)), ]
    
    IW_ACS_PWDIS_quant[[i]] <- ORM_IW_ACS[which(is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                                  !is.na(ORM_IW_ACS$PWDIS_adj) & 
                                                  ORM_IW_ACS$PWDIS_adj>0), ]
    
    IW_ACS_PRMP_quant[[i]] <- ORM_IW_ACS[which(is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                                 !is.na(ORM_IW_ACS$PRMP_adj) & 
                                                 ORM_IW_ACS$PRMP_adj>0), ]
    
    IW_ACS_EJ_quant[[i]] <- ORM_IW_ACS[which(is.na(ORM_IW_ACS$WATERS_CHANGE) & 
                                               !is.na(ORM_IW_ACS$PRMP_adj) & 
                                               ORM_IW_ACS$PRMP_adj>0 & 
                                               !is.na(ORM_IW_ACS$PWDIS_adj) & 
                                               ORM_IW_ACS$PWDIS_adj>0), ]
    
  }
  
  #*# AL (11/1) - no longer based on percentiles
  # IW_ACS_totals[[i]] <- data.frame(t(colSums(IW_ACS_quant[[i]]%>%select(P80_MIN:P80_DEM_RSEI_RMP,P80_TRIBAL), na.rm=TRUE)))
  IW_ACS_totals[[i]] <- data.frame(t(colSums(IW_ACS_quant[[i]]%>%select(TOT_AGE_POP:LINGISO_HH,TOT_RACE_POP_CT:TRIBE_POP),
                                             na.rm=T)))
  IW_ACS_totals[[i]]$PERC_MIN <- with(IW_ACS_totals[[i]], MINOR_POP/TOT_RACE_POP*100)
  IW_ACS_totals[[i]]$PERC_POV <- with(IW_ACS_totals[[i]], POV_POP/TOT_POV_POP*100)
  IW_ACS_totals[[i]]$PERC_LESSHS <- with(IW_ACS_totals[[i]], EDU_POP/TOT_EDU_POP*100)
  IW_ACS_totals[[i]]$PERC_LING <- with(IW_ACS_totals[[i]], LINGISO_HH/TOT_LINGISO_HH*100)
  IW_ACS_totals[[i]]$PERC_UNDER5 <- with(IW_ACS_totals[[i]], UNDER_5/TOT_AGE_POP*100)
  IW_ACS_totals[[i]]$PERC_OVER64 <- with(IW_ACS_totals[[i]], OVER_64/TOT_AGE_POP*100)
  IW_ACS_totals[[i]]$PERC_TRIBAL <- with(IW_ACS_totals[[i]], TRIBE_POP/TOT_RACE_POP_CT*100)
  
  IW_ACS_totals[[i]]$PWDIS <- mean(IW_ACS_quant[[i]]$PWDIS_adj, na.rm=T)
  IW_ACS_totals[[i]]$PRMP <- mean(IW_ACS_quant[[i]]$PRMP_adj, na.r=T)
  
  IW_ACS_totals[[i]]$HUC12_Num <- length(unique(IW_ACS_quant[[i]]$HUC12))
  IW_ACS_totals[[i]]$ORM_Total <- length(unique(ORM_impactedwaters$HUC12))
  
  write.csv(IW_ACS_totals[[i]],paste0(EJ,"R_programs/Outputs/IW_",i,"_ACS.csv"),row.names=FALSE)
  
  IW_ACS_PWDIS_totals[[i]] <- data.frame(t(colSums(IW_ACS_PWDIS_quant[[i]]%>%select(TOT_AGE_POP:LINGISO_HH,TOT_RACE_POP_CT:TRIBE_POP),
                                                   na.rm=T)))
  IW_ACS_PWDIS_totals[[i]]$PERC_MIN <- with(IW_ACS_PWDIS_totals[[i]], MINOR_POP/TOT_RACE_POP*100)
  IW_ACS_PWDIS_totals[[i]]$PERC_POV <- with(IW_ACS_PWDIS_totals[[i]], POV_POP/TOT_POV_POP*100)
  IW_ACS_PWDIS_totals[[i]]$PERC_LESSHS <- with(IW_ACS_PWDIS_totals[[i]], EDU_POP/TOT_EDU_POP*100)
  IW_ACS_PWDIS_totals[[i]]$PERC_LING <- with(IW_ACS_PWDIS_totals[[i]], LINGISO_HH/TOT_LINGISO_HH*100)
  IW_ACS_PWDIS_totals[[i]]$PERC_UNDER5 <- with(IW_ACS_PWDIS_totals[[i]], UNDER_5/TOT_AGE_POP*100)
  IW_ACS_PWDIS_totals[[i]]$PERC_OVER64 <- with(IW_ACS_PWDIS_totals[[i]], OVER_64/TOT_AGE_POP*100)
  IW_ACS_PWDIS_totals[[i]]$PERC_TRIBAL <- with(IW_ACS_PWDIS_totals[[i]], TRIBE_POP/TOT_RACE_POP_CT*100)
  
  IW_ACS_PWDIS_totals[[i]]$PWDIS <- mean(IW_ACS_PWDIS_quant[[i]]$PWDIS_adj, na.rm=T)
  IW_ACS_PWDIS_totals[[i]]$PRMP <- mean(IW_ACS_PWDIS_quant[[i]]$PRMP_adj, na.r=T)
  
  IW_ACS_PWDIS_totals[[i]]$HUC12_Num <- length(unique(IW_ACS_PWDIS_quant[[i]]$HUC12))
  IW_ACS_PWDIS_totals[[i]]$ORM_Total <- length(unique(ORM_impactedwaters$HUC12))
  
  write.csv(IW_ACS_PWDIS_totals[[i]],paste0(EJ,"R_programs/Outputs/IW_PWDIS_",i,"_ACS.csv"),row.names=FALSE)
  
  IW_ACS_PRMP_totals[[i]] <- data.frame(t(colSums(IW_ACS_PRMP_quant[[i]]%>%select(TOT_AGE_POP:LINGISO_HH,TOT_RACE_POP_CT:TRIBE_POP),
                                                  na.rm=T)))
  IW_ACS_PRMP_totals[[i]]$PERC_MIN <- with(IW_ACS_PRMP_totals[[i]], MINOR_POP/TOT_RACE_POP*100)
  IW_ACS_PRMP_totals[[i]]$PERC_POV <- with(IW_ACS_PRMP_totals[[i]], POV_POP/TOT_POV_POP*100)
  IW_ACS_PRMP_totals[[i]]$PERC_LESSHS <- with(IW_ACS_PRMP_totals[[i]], EDU_POP/TOT_EDU_POP*100)
  IW_ACS_PRMP_totals[[i]]$PERC_LING <- with(IW_ACS_PRMP_totals[[i]], LINGISO_HH/TOT_LINGISO_HH*100)
  IW_ACS_PRMP_totals[[i]]$PERC_UNDER5 <- with(IW_ACS_PRMP_totals[[i]], UNDER_5/TOT_AGE_POP*100)
  IW_ACS_PRMP_totals[[i]]$PERC_OVER64 <- with(IW_ACS_PRMP_totals[[i]], OVER_64/TOT_AGE_POP*100)
  IW_ACS_PRMP_totals[[i]]$PERC_TRIBAL <- with(IW_ACS_PRMP_totals[[i]], TRIBE_POP/TOT_RACE_POP_CT*100)
  
  IW_ACS_PRMP_totals[[i]]$PWDIS <- mean(IW_ACS_PRMP_quant[[i]]$PWDIS_adj, na.rm=T)
  IW_ACS_PRMP_totals[[i]]$PRMP <- mean(IW_ACS_PRMP_quant[[i]]$PRMP_adj, na.r=T)
  
  IW_ACS_PRMP_totals[[i]]$HUC12_Num <- length(unique(IW_ACS_PRMP_quant[[i]]$HUC12))
  IW_ACS_PRMP_totals[[i]]$ORM_Total <- length(unique(ORM_impactedwaters$HUC12))
  
  write.csv(IW_ACS_PRMP_totals[[i]],paste0(EJ,"R_programs/Outputs/IW_PRMP_",i,"_ACS.csv"),row.names=FALSE)
  
  IW_ACS_EJ_totals[[i]] <- data.frame(t(colSums(IW_ACS_EJ_quant[[i]]%>%select(TOT_AGE_POP:LINGISO_HH,TOT_RACE_POP_CT:TRIBE_POP),
                                                na.rm=T)))
  IW_ACS_EJ_totals[[i]]$PERC_MIN <- with(IW_ACS_EJ_totals[[i]], MINOR_POP/TOT_RACE_POP*100)
  IW_ACS_EJ_totals[[i]]$PERC_POV <- with(IW_ACS_EJ_totals[[i]], POV_POP/TOT_POV_POP*100)
  IW_ACS_EJ_totals[[i]]$PERC_LESSHS <- with(IW_ACS_EJ_totals[[i]], EDU_POP/TOT_EDU_POP*100)
  IW_ACS_EJ_totals[[i]]$PERC_LING <- with(IW_ACS_EJ_totals[[i]], LINGISO_HH/TOT_LINGISO_HH*100)
  IW_ACS_EJ_totals[[i]]$PERC_UNDER5 <- with(IW_ACS_EJ_totals[[i]], UNDER_5/TOT_AGE_POP*100)
  IW_ACS_EJ_totals[[i]]$PERC_OVER64 <- with(IW_ACS_EJ_totals[[i]], OVER_64/TOT_AGE_POP*100)
  IW_ACS_EJ_totals[[i]]$PERC_TRIBAL <- with(IW_ACS_EJ_totals[[i]], TRIBE_POP/TOT_RACE_POP_CT*100)
  
  IW_ACS_EJ_totals[[i]]$PWDIS <- mean(IW_ACS_EJ_quant[[i]]$PWDIS_adj, na.rm=T)
  IW_ACS_EJ_totals[[i]]$PRMP <- mean(IW_ACS_EJ_quant[[i]]$PRMP_adj, na.r=T)
  
  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))
  
  write.csv(IW_ACS_EJ_totals[[i]],paste0(EJ,"R_programs/Outputs/IW_EJ_",i,"_ACS.csv"),row.names=FALSE)
  
}

#*# KM: Confirmed that sum of all HUC12s in IW_ACS_totals[[1:5]] = length(unique(ORM_IW_ACS$HUC12))
