# Jessica Balukas, ICF, September 2021
# Script to process ORM data (2010-2019), calculate section 404 impacts with mitigation requirements and single-water permit counts for Rapanos, NWPR, and the change
# Federalism Scenario: removes permits issued in states that "regulate waters more broadly than the CWA requires."
# NOTE: TO ADJUST FOR DIFFERENT EMPLOYEE ID/FOLDER SYNCING, UPDATE CODE SIGNIFIED WITH #*# (search for combo to get all occurrences)
# NOTE: File paths assume syncing of entire WOTUS Reconsideration general folder

# Install Packages
#install.packages("tidyverse")
#install.packages("stringr")

# Load Packages
library(tidyverse)
library(stringr)

# Folder locations
var <- "40179" #*#
IN_DIR <- (paste0("C:/Users/",var,"/ICF/WOTUS Reconsideration - General/ORM Data Processing/Inputs")) #*#
OUT_DIR <- (paste0("C:/Users/",var,"/ICF/WOTUS Reconsideration - General/ORM Data Processing/Outputs")) #*#
QA_DIR <- (paste0("C:/Users/",var,"/ICF/WOTUS Reconsideration - General/ORM Data Processing/Outputs/QA")) #*#

# Read input files
ORM_Dataset <- read.csv(file.path(IN_DIR,"ORM_Data_2010-2019_ReducedFields.csv"))
CWR_Filter <- read.csv(file.path(IN_DIR,"filter_for_CWR_ICF_Edits.csv"))
Perc_Permitted_Waters <- read.csv(file.path(IN_DIR,"Perc_Permitted_Waters.csv"))
HUC12 <- read.delim(file.path(IN_DIR,"ORMPermit_Georef_0907.txt"))
HUC12_AK <- read.delim(file.path(IN_DIR,"ORMPermit_AK_Georef_0907.txt"))
Prop_10YrPeriod_NonCWR <- read.csv(file.path(IN_DIR,"Proportion10YrPeriod_NonCWR.csv"))
StateMatrix <- read.csv(file.path(IN_DIR,"State Response Matrix.csv"))

# QA - HUC12 obs. lower than ORM_Dataset obs.
HUC12_ID <- HUC12[c("PERMIT_ID")]
HUC12_ID$HUC12_Flag <- 1
ID_test <- merge(x=ORM_Dataset,y=HUC12_ID,by="PERMIT_ID",all.x=TRUE)
ID_test_QA <- ID_test[is.na(ID_test$HUC12_Flag), ]
write.csv(ID_test_QA,paste(QA_DIR,"HUC12_MissingIDs.csv",sep="/"),row.names=FALSE)
rm(list = c('HUC12_ID','ID_test','ID_test_QA'))

# Ensure that HUC12 values are always 12 digits (add leading zeros when necessary), add HUC4 variable
names(HUC12)[names(HUC12) == "WBD_HUC12"] <- "HUC12_Main"
HUC12$HUC12 <- str_pad(HUC12$HUC12_Main,12, pad = "0")
HUC12$HUC4 <- substr(HUC12$HUC12,start=1,stop=4)

# Convert HUC12 and HUC4 variables to character
HUC12$HUC12 <- as.character(HUC12$HUC12)
HUC12$HUC4 <- as.character(HUC12$HUC4)

# Merge HUC12 to ORM_Dataset
HUC12 <- HUC12[c("PERMIT_ID", "HUC12", "HUC4", "STUSPS")]
ORM_Dataset_Final <- merge(x=ORM_Dataset,y=HUC12,by="PERMIT_ID",all.x=TRUE)

# Create combined STATE variable (prioritize spatial state, then STATE_AR, then STATE)
ORM_Dataset_Final$STATE_FINAL <- ifelse(ORM_Dataset_Final$STUSPS=='' | is.na(ORM_Dataset_Final$STUSPS),ORM_Dataset_Final$STATE_AR,ORM_Dataset_Final$STUSPS)
ORM_Dataset_Final$STATE_FINAL <- ifelse(ORM_Dataset_Final$STATE_FINAL=='' | is.na(ORM_Dataset_Final$STATE_FINAL),ORM_Dataset_Final$STATE,ORM_Dataset_Final$STATE_FINAL)

# Check for remaining missing values in STATE_FINAL variable
STATE_QA  <- ORM_Dataset_Final[ which(ORM_Dataset_Final$STATE_FINAL ==''), ]
write.csv(STATE_QA,paste(QA_DIR,"STATE_MISSING.csv",sep="/"),row.names=FALSE)

# Merge HUC12_AK to ORM_Dataset_Final, replace HUC12 with AK value if STATE_FINAL=AK (Alaska HUC12 merge processed separately)
names(HUC12_AK)[names(HUC12_AK) == "WBD_HUC12"] <- "HUC12_AK"
HUC12_AK$HUC12_AK <- str_pad(HUC12_AK$HUC12_AK,12, pad = "0")
HUC12_AK$HUC4_AK <- substr(HUC12_AK$HUC12_AK,start=1,stop=4)
HUC12_AK$HUC12_AK <- as.character(HUC12_AK$HUC12_AK)
HUC12_AK$HUC4_AK <- as.character(HUC12_AK$HUC4_AK)
HUC12_AK <- HUC12_AK[c("PERMIT_ID", "HUC12_AK","HUC4_AK")]
ORM_Dataset_Final <- merge(x=ORM_Dataset_Final,y=HUC12_AK,by="PERMIT_ID",all.x=TRUE)
ORM_Dataset_Final$HUC12 <- ifelse(ORM_Dataset_Final$STATE_FINAL=='AK',ORM_Dataset_Final$HUC12_AK,ORM_Dataset_Final$HUC12)
ORM_Dataset_Final$HUC4 <- ifelse(ORM_Dataset_Final$STATE_FINAL=='AK',ORM_Dataset_Final$HUC4_AK,ORM_Dataset_Final$HUC4)

# QA - obs. with missing HUC12
HUC12_QA <- ORM_Dataset_Final[is.na(ORM_Dataset_Final$HUC12), ]
write.csv(HUC12_QA,paste(QA_DIR,"HUC12_MissingValue.csv",sep="/"),row.names=FALSE)
ORM_Dataset_Final[ ,c('HUC12_AK', 'HUC4_AK')] <- list(NULL)

### Merge remaining datasets to ORM_Dataset
# CWR Filter (to remove permits issued under Clean Water Rule Guidance)
CWR_Filter$State_Date <- paste(CWR_Filter$state,"_",CWR_Filter$END_DATE)
CWR_Filter <- CWR_Filter[c("State_Date","CWR")]
ORM_Dataset_Final$State_Date <- paste(ORM_Dataset_Final$STATE_FINAL,"_",ORM_Dataset_Final$END_DATE)
ORM_Dataset_Final <- merge(x=ORM_Dataset_Final,y=CWR_Filter,by="State_Date",all.x=TRUE)

# Check CWR Filter Merge
#write.csv(ORM_Dataset_Final,paste(QA_DIR,"CWR_Filter_Merge.csv",sep="/"),row.names=FALSE)
ORM_Dataset_Final$State_Date <- NULL

# Perc_Permitted_Waters (jurisdictional determination corrections)
ORM_Dataset_Final <- merge(x=ORM_Dataset_Final,y=Perc_Permitted_Waters,by="COWARDIN_FINAL",all.x=TRUE)
#write.csv(ORM_Dataset_Final,paste(QA_DIR,"Perc_Permitted_Waters_Merge.csv",sep="/"),row.names=FALSE)

# Prop_10YrPeriod_NonCWR
ORM_Dataset_Final <- merge(x=ORM_Dataset_Final,y=Prop_10YrPeriod_NonCWR,by="STATE_FINAL",all.x=TRUE)
#write.csv(ORM_Dataset_Final,paste(QA_DIR,"Prop_10YrPeriod_NonCWR_Merge.csv",sep="/"),row.names=FALSE)

# StateMatrix
ORM_Dataset_Final <- merge(x=ORM_Dataset_Final,y=StateMatrix,by.x="STATE_FINAL",by.y="State",all.x=TRUE)
write.csv(ORM_Dataset_Final,paste(QA_DIR,"StateMatrix_Merge.csv",sep="/"),row.names=FALSE)

# Drop intermediate datasets
rm(list = c('ORM_Dataset','CWR_Filter','Perc_Permitted_Waters','HUC12','HUC12_AK','HUC12_QA','Prop_10YrPeriod_NonCWR','STATE_QA','StateMatrix'))

# Remove CWR permits and projects with a primarily mitigation-based purpose from analysis
ORM_Dataset_Final_NoCWR  <- ORM_Dataset_Final[ which(ORM_Dataset_Final$CWR==0), ]
ORM_Dataset_Final_NoCWRNoMit  <- ORM_Dataset_Final_NoCWR[ which(ORM_Dataset_Final_NoCWR$NonMitigation_Work_Type==1), ]

# Remove permits in states that "regulate waters more broadly than the CWA requires" (StateMatrix)
ORM_Dataset_Final_NoCWRNoMit  <- ORM_Dataset_Final_NoCWRNoMit[ which(ORM_Dataset_Final_NoCWRNoMit$X404_BroaderThanFedRegs != 1), ]

# Create binary variable to indicate permit type (1=individual permt, 0=general permit)
ORM_Dataset_Final_NoCWRNoMit %>% count(ACTION_TYPE)
ORM_Dataset_Final_NoCWRNoMit$Ind_Permit <- ifelse(ORM_Dataset_Final_NoCWRNoMit$ACTION_TYPE=='LOP',1,
                                           ifelse(ORM_Dataset_Final_NoCWRNoMit$ACTION_TYPE=='SP',1,0))

# Create binary variables to indicate single water permits (1=single water affected, 0=multiple waters affected)
PermitID_Count <- ORM_Dataset_Final_NoCWRNoMit %>% count(PERMIT_ID2)
ORM_Dataset_Final_NoCWRNoMit <- merge(x=ORM_Dataset_Final_NoCWRNoMit,y=PermitID_Count,by="PERMIT_ID2",all.x=TRUE)
rm(list = c('ORM_Dataset_Final','ORM_Dataset_Final_NoCWR','PermitID_Count'))
ORM_Dataset_Final_NoCWRNoMit$SingleWater <- ifelse(ORM_Dataset_Final_NoCWRNoMit$n>1,0,1)
ORM_Dataset_Final_NoCWRNoMit$n <- NULL

# Create SI, SG, MI, MG binary variables (described on EPA's technical direction)
ORM_Dataset_Final_NoCWRNoMit$SI <- ifelse(ORM_Dataset_Final_NoCWRNoMit$Ind_Permit==1 & ORM_Dataset_Final_NoCWRNoMit$SingleWater==1,1,0)
ORM_Dataset_Final_NoCWRNoMit$SG <- ifelse(ORM_Dataset_Final_NoCWRNoMit$Ind_Permit==0 & ORM_Dataset_Final_NoCWRNoMit$SingleWater==1,1,0)
ORM_Dataset_Final_NoCWRNoMit$MI <- ifelse(ORM_Dataset_Final_NoCWRNoMit$Ind_Permit==1 & ORM_Dataset_Final_NoCWRNoMit$SingleWater==0,1,0)
ORM_Dataset_Final_NoCWRNoMit$MG <- ifelse(ORM_Dataset_Final_NoCWRNoMit$Ind_Permit==0 & ORM_Dataset_Final_NoCWRNoMit$SingleWater==0,1,0)

# Calculate SI, SG, MI, MG percentages by HUC4/COWARDIN_FINAL (divide by # affected waters)
# Determine unique number of permits by HUC4/COWARDIN_FINAL for later permit analysis
count_SI_SG_MI_MG <- ORM_Dataset_Final_NoCWRNoMit %>%  
  group_by(HUC4,COWARDIN_FINAL) %>%
  summarise(Permit_Count = n_distinct(PERMIT_ID2),
            Waters_Count = n(),
            SI_Count = sum(SI),
            SG_Count = sum(SG),
            MI_Count = sum(MI),
            MG_Count = sum(MG)) %>%
  mutate(SI_Perc = SI_Count/Waters_Count,
         SG_Perc = SG_Count/Waters_Count,
         MI_Perc = MI_Count/Waters_Count,
         MG_Perc = MG_Count/Waters_Count)
count_SI_SG_MI_MG$HUC4_COWARDIN <- paste(count_SI_SG_MI_MG$HUC4,"_",count_SI_SG_MI_MG$COWARDIN_FINAL)

count_SI_SG_MI_MG_2 <- count_SI_SG_MI_MG[c("HUC4_COWARDIN","Permit_Count","SI_Perc","SG_Perc","MI_Perc","MG_Perc")]

# Merge SI, SG, MI, MG percentages with ORM_Dataset_Final_NoCWRNoMit dataset
ORM_Dataset_Final_NoCWRNoMit$HUC4_COWARDIN <- paste(ORM_Dataset_Final_NoCWRNoMit$HUC4,"_",ORM_Dataset_Final_NoCWRNoMit$COWARDIN_FINAL)
ORM_Dataset_Final_NoCWRNoMit <- merge(x=ORM_Dataset_Final_NoCWRNoMit,y=count_SI_SG_MI_MG_2,by="HUC4_COWARDIN",all.x=TRUE)
rm(list = c('count_SI_SG_MI_MG','count_SI_SG_MI_MG_2'))

# Create dummy variable for mitigation (1=mitigation required, 0=mitigation NOT required)
ORM_Dataset_Final_NoCWRNoMit$MitReq <- ifelse(ORM_Dataset_Final_NoCWRNoMit$MITIGATION_TYPES!='',1,0)

# Create dummy variables for permanent and temporary impacts (1=has impacts, 0=no impacts)
ORM_Dataset_Final_NoCWRNoMit$Perm_Impact <- ifelse(ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_ACRES_PERM>0 | ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_PERM>0,1,0)
ORM_Dataset_Final_NoCWRNoMit$Temp_Impact <- ifelse(ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_ACRES_TEMP>0 | ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_TEMP>0,1,0)
ORM_Dataset_Final_NoCWRNoMit$PermTemp_Impact <- ifelse(ORM_Dataset_Final_NoCWRNoMit$Perm_Impact==1 | ORM_Dataset_Final_NoCWRNoMit$Temp_Impact==1,1,0)
ORM_Dataset_Final_NoCWRNoMit$Perm_Impact <- ifelse(is.na(ORM_Dataset_Final_NoCWRNoMit$Perm_Impact),0,ORM_Dataset_Final_NoCWRNoMit$Perm_Impact)
ORM_Dataset_Final_NoCWRNoMit$Temp_Impact <- ifelse(is.na(ORM_Dataset_Final_NoCWRNoMit$Temp_Impact),0,ORM_Dataset_Final_NoCWRNoMit$Temp_Impact)
ORM_Dataset_Final_NoCWRNoMit$PermTemp_Impact <- ifelse(is.na(ORM_Dataset_Final_NoCWRNoMit$PermTemp_Impact),0,ORM_Dataset_Final_NoCWRNoMit$PermTemp_Impact)
#write.csv(ORM_Dataset_Final_NoCWRNoMit,paste(QA_DIR,"PermTempImpactDummies_Federalism.csv",sep="/"),row.names=FALSE)

# Create PM, TM, PN, TN binary variables (described on EPA's technical direction)
ORM_Dataset_Final_NoCWRNoMit$PM <- ifelse(ORM_Dataset_Final_NoCWRNoMit$Perm_Impact==1 & ORM_Dataset_Final_NoCWRNoMit$MitReq==1,1,0)
ORM_Dataset_Final_NoCWRNoMit$TM <- ifelse(ORM_Dataset_Final_NoCWRNoMit$Temp_Impact==1 & ORM_Dataset_Final_NoCWRNoMit$MitReq==1,1,0)
ORM_Dataset_Final_NoCWRNoMit$PN <- ifelse(ORM_Dataset_Final_NoCWRNoMit$Perm_Impact==1 & ORM_Dataset_Final_NoCWRNoMit$MitReq==0,1,0)
ORM_Dataset_Final_NoCWRNoMit$TN <- ifelse(ORM_Dataset_Final_NoCWRNoMit$Temp_Impact==1 & ORM_Dataset_Final_NoCWRNoMit$MitReq==0,1,0)
#write.csv(ORM_Dataset_Final_NoCWRNoMit,paste(QA_DIR,"PM_TM_PN_TN_Federalism.csv",sep="/"),row.names=FALSE)

# Calculate PM, TM, PN, TN percentages by HUC4/COWARDIN_FINAL (divide by # affected waters)
count_PM_TM_PN_TN <- ORM_Dataset_Final_NoCWRNoMit %>%  
  group_by(HUC4,COWARDIN_FINAL) %>%
  summarise(Waters_Count = n(),
            PM_Count = sum(PM),
            TM_Count = sum(TM),
            PN_Count = sum(PN),
            TN_Count = sum(TN)) %>%
  mutate(PM_Perc = PM_Count/Waters_Count,
         TM_Perc = TM_Count/Waters_Count,
         PN_Perc = PN_Count/Waters_Count,
         TN_Perc = TN_Count/Waters_Count)
count_PM_TM_PN_TN$HUC4_COWARDIN <- paste(count_PM_TM_PN_TN$HUC4,"_",count_PM_TM_PN_TN$COWARDIN_FINAL)

count_PM_TM_PN_TN_2 <- count_PM_TM_PN_TN[c("HUC4_COWARDIN","Waters_Count","PM_Perc","TM_Perc","PN_Perc","TN_Perc")]

# Merge PM, TM, PN, TN percentages with ORM_Dataset_Final_NoCWRNoMit dataset
ORM_Dataset_Final_NoCWRNoMit <- merge(x=ORM_Dataset_Final_NoCWRNoMit,y=count_PM_TM_PN_TN_2,by="HUC4_COWARDIN",all.x=TRUE)
rm(list = c('count_PM_TM_PN_TN','count_PM_TM_PN_TN_2'))

# Convert LF (permanent & temporary) to acres, using average buffer width of 50 ft. (25 ft. on each side of the stream) 
ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_PERM_ACRE_CONV <- ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_PERM*50*0.0000229568
ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_TEMP_ACRE_CONV <- ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_TEMP*50*0.0000229568

# Calculate total permanent and temporary impacts (acres + LF converted to acres)
ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_ACRES_PERM <- ifelse(is.na(ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_ACRES_PERM),0,ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_ACRES_PERM)
ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_PERM_ACRE_CONV <- ifelse(is.na(ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_PERM_ACRE_CONV),0,ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_PERM_ACRE_CONV)
ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_PERM_TOTAL <- ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_ACRES_PERM+ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_PERM_ACRE_CONV

ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_ACRES_TEMP <- ifelse(is.na(ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_ACRES_TEMP),0,ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_ACRES_TEMP)
ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_TEMP_ACRE_CONV <- ifelse(is.na(ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_TEMP_ACRE_CONV),0,ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_TEMP_ACRE_CONV)
ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_TEMP_TOTAL <- ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_ACRES_TEMP+ORM_Dataset_Final_NoCWRNoMit$IMPACT_AUTH_LF_TEMP_ACRE_CONV

################################################################
#Calculate impacts not covered by a permit under Rapanos & NWPR
################################################################

# Remove observations with HIST1, missing COWARDIN_FINAL, or missing HUC12
ORM_Dataset_Final_NoCWRNoMit  <- ORM_Dataset_Final_NoCWRNoMit[ which(ORM_Dataset_Final_NoCWRNoMit$COWARDIN_FINAL != "HIST1"), ]
ORM_Dataset_Final_NoCWRNoMit  <- ORM_Dataset_Final_NoCWRNoMit[ which(ORM_Dataset_Final_NoCWRNoMit$COWARDIN_FINAL != ''), ]
ORM_Dataset_Final_NoCWRNoMit <- ORM_Dataset_Final_NoCWRNoMit[!(is.na(ORM_Dataset_Final_NoCWRNoMit$HUC12)), ]

# Multiply impacts by Cowardin Code-level Rapanos and NWPR adjustments (to account for impacts not covered under Rapanos vs. NWPR)
# Multiply Rapanos/NWPR-corrected impacts by HUC4/Cowardin Code mitigation percentages (PM_Perc & TM_Perc)
# Divide Rapanos/NWPR- & PM/TM-corrected impacts by proportion of 10-year period that state followed Rapanos to get annual average
# Use the same methodology to count impacted waters not covered under Rapanos vs. NWPR
ORM_Dataset_Final_NoCWRNoMit <- ORM_Dataset_Final_NoCWRNoMit %>%  
  mutate(IMPACT_AUTH_ACRES_PERM_PM_Rapanos_10Yr = (((IMPACT_AUTH_ACRES_PERM/Rapanos_Perc_Permitted)-IMPACT_AUTH_ACRES_PERM)*PM_Perc)/Prop_10YrPeriod_NonCWR,
         IMPACT_AUTH_LF_PERM_ACRE_CONV_PM_Rapanos_10Yr = (((IMPACT_AUTH_LF_PERM_ACRE_CONV/Rapanos_Perc_Permitted)-IMPACT_AUTH_LF_PERM_ACRE_CONV)*PM_Perc)/Prop_10YrPeriod_NonCWR,
         IMPACT_AUTH_PERM_TOTAL_PM_Rapanos_10Yr = (((IMPACT_AUTH_PERM_TOTAL/Rapanos_Perc_Permitted)-IMPACT_AUTH_PERM_TOTAL)*PM_Perc)/Prop_10YrPeriod_NonCWR,
         IMPACT_AUTH_ACRES_TEMP_TM_Rapanos_10Yr = (((IMPACT_AUTH_ACRES_TEMP/Rapanos_Perc_Permitted)-IMPACT_AUTH_ACRES_TEMP)*TM_Perc)/Prop_10YrPeriod_NonCWR,
         IMPACT_AUTH_LF_TEMP_ACRE_CONV_TM_Rapanos_10Yr = (((IMPACT_AUTH_LF_TEMP_ACRE_CONV/Rapanos_Perc_Permitted)-IMPACT_AUTH_LF_TEMP_ACRE_CONV)*TM_Perc)/Prop_10YrPeriod_NonCWR,
         IMPACT_AUTH_TEMP_TOTAL_TM_Rapanos_10Yr = (((IMPACT_AUTH_TEMP_TOTAL/Rapanos_Perc_Permitted)-IMPACT_AUTH_TEMP_TOTAL)*TM_Perc)/Prop_10YrPeriod_NonCWR,
         IMPACT_AUTH_ACRES_PERM_PM_NWPR_10Yr = (((IMPACT_AUTH_ACRES_PERM/Rapanos_Perc_Permitted)-(IMPACT_AUTH_ACRES_PERM/Rapanos_Perc_Permitted)*NWPR_Perc_Permitted)*PM_Perc)/Prop_10YrPeriod_NonCWR,
         IMPACT_AUTH_LF_PERM_ACRE_CONV_PM_NWPR_10Yr = (((IMPACT_AUTH_LF_PERM_ACRE_CONV/Rapanos_Perc_Permitted)-(IMPACT_AUTH_LF_PERM_ACRE_CONV/Rapanos_Perc_Permitted)*NWPR_Perc_Permitted)*PM_Perc)/Prop_10YrPeriod_NonCWR,
         IMPACT_AUTH_PERM_TOTAL_PM_NWPR_10Yr = (((IMPACT_AUTH_PERM_TOTAL/Rapanos_Perc_Permitted)-(IMPACT_AUTH_PERM_TOTAL/Rapanos_Perc_Permitted)*NWPR_Perc_Permitted)*PM_Perc)/Prop_10YrPeriod_NonCWR,
         IMPACT_AUTH_ACRES_TEMP_TM_NWPR_10Yr = (((IMPACT_AUTH_ACRES_TEMP/Rapanos_Perc_Permitted)-(IMPACT_AUTH_ACRES_TEMP/Rapanos_Perc_Permitted)*NWPR_Perc_Permitted)*TM_Perc)/Prop_10YrPeriod_NonCWR,
         IMPACT_AUTH_LF_TEMP_ACRE_CONV_TM_NWPR_10Yr = (((IMPACT_AUTH_LF_TEMP_ACRE_CONV/Rapanos_Perc_Permitted)-(IMPACT_AUTH_LF_TEMP_ACRE_CONV/Rapanos_Perc_Permitted)*NWPR_Perc_Permitted)*TM_Perc)/Prop_10YrPeriod_NonCWR,
         IMPACT_AUTH_TEMP_TOTAL_TM_NWPR_10Yr = (((IMPACT_AUTH_TEMP_TOTAL/Rapanos_Perc_Permitted)-(IMPACT_AUTH_TEMP_TOTAL/Rapanos_Perc_Permitted)*NWPR_Perc_Permitted)*TM_Perc)/Prop_10YrPeriod_NonCWR,
         PERM_Count_Rapanos_10Yr = ((Perm_Impact/Rapanos_Perc_Permitted)-Perm_Impact)/Prop_10YrPeriod_NonCWR,
         PERM_Count_NWPR_10Yr = ((Perm_Impact/Rapanos_Perc_Permitted)-(Perm_Impact/Rapanos_Perc_Permitted)*NWPR_Perc_Permitted)/Prop_10YrPeriod_NonCWR,
         PERM_Count_Change_10Yr = (Perm_Impact-(Perm_Impact*(NWPR_Perc_Permitted/Rapanos_Perc_Permitted)))/Prop_10YrPeriod_NonCWR,
         TEMP_Count_Rapanos_10Yr = ((Temp_Impact/Rapanos_Perc_Permitted)-Temp_Impact)/Prop_10YrPeriod_NonCWR,
         TEMP_Count_NWPR_10Yr = ((Temp_Impact/Rapanos_Perc_Permitted)-(Temp_Impact/Rapanos_Perc_Permitted)*NWPR_Perc_Permitted)/Prop_10YrPeriod_NonCWR,
         TEMP_Count_Change_10Yr = (Temp_Impact-(Temp_Impact*(NWPR_Perc_Permitted/Rapanos_Perc_Permitted)))/Prop_10YrPeriod_NonCWR,
         ImpactedWater_Count_Rapanos_10Yr = ((PermTemp_Impact/Rapanos_Perc_Permitted)-PermTemp_Impact)/Prop_10YrPeriod_NonCWR,
         ImpactedWater_Count_NWPR_10Yr = ((PermTemp_Impact/Rapanos_Perc_Permitted)-(PermTemp_Impact/Rapanos_Perc_Permitted)*NWPR_Perc_Permitted)/Prop_10YrPeriod_NonCWR,
         ImpactedWater_Count_Change_10Yr = (PermTemp_Impact-(PermTemp_Impact*(NWPR_Perc_Permitted/Rapanos_Perc_Permitted)))/Prop_10YrPeriod_NonCWR)

# Aggregate permanent & temporary impacts at the HUC12 level
Impacts_HUC12 <- ORM_Dataset_Final_NoCWRNoMit %>%  
  group_by(HUC12) %>%
  summarise(PERM_ACRES_Rapanos = sum(IMPACT_AUTH_ACRES_PERM_PM_Rapanos_10Yr),
            PERM_LF_ACRE_CONV_Rapanos = sum(IMPACT_AUTH_LF_PERM_ACRE_CONV_PM_Rapanos_10Yr),
            PERM_TOTAL_Rapanos = sum(IMPACT_AUTH_PERM_TOTAL_PM_Rapanos_10Yr),
            TEMP_ACRES_Rapanos = sum(IMPACT_AUTH_ACRES_TEMP_TM_Rapanos_10Yr),
            TEMP_LF_ACRE_CONV_Rapanos = sum(IMPACT_AUTH_LF_TEMP_ACRE_CONV_TM_Rapanos_10Yr),
            TEMP_TOTAL_Rapanos = sum(IMPACT_AUTH_TEMP_TOTAL_TM_Rapanos_10Yr),
            PERM_ACRES_NWPR = sum(IMPACT_AUTH_ACRES_PERM_PM_NWPR_10Yr),
            PERM_LF_ACRE_CONV_NWPR = sum(IMPACT_AUTH_LF_PERM_ACRE_CONV_PM_NWPR_10Yr),
            PERM_TOTAL_NWPR = sum(IMPACT_AUTH_PERM_TOTAL_PM_NWPR_10Yr),
            TEMP_ACRES_NWPR = sum(IMPACT_AUTH_ACRES_TEMP_TM_NWPR_10Yr),
            TEMP_LF_ACRE_CONV_NWPR = sum(IMPACT_AUTH_LF_TEMP_ACRE_CONV_TM_NWPR_10Yr),
            TEMP_TOTAL_NWPR = sum(IMPACT_AUTH_TEMP_TOTAL_TM_NWPR_10Yr)) %>%
  mutate(PERM_ACRES_CHANGE = PERM_ACRES_Rapanos - PERM_ACRES_NWPR,
         PERM_LF_ACRE_CONV_CHANGE = PERM_LF_ACRE_CONV_Rapanos - PERM_LF_ACRE_CONV_NWPR,
         PERM_TOTAL_CHANGE = PERM_TOTAL_Rapanos - PERM_TOTAL_NWPR,
         TEMP_ACRES_CHANGE = TEMP_ACRES_Rapanos - TEMP_ACRES_NWPR,
         TEMP_LF_ACRE_CONV_CHANGE = TEMP_LF_ACRE_CONV_Rapanos - TEMP_LF_ACRE_CONV_NWPR,
         TEMP_TOTAL_CHANGE = TEMP_TOTAL_Rapanos - TEMP_TOTAL_NWPR)
write.csv(Impacts_HUC12,paste(OUT_DIR,"Impacts_By_HUC12_Federalism.csv",sep="/"),row.names=FALSE)

# Aggregate permanent & temporary impacts at the HUC4 level
Impacts_HUC4 <- ORM_Dataset_Final_NoCWRNoMit %>%  
  group_by(HUC4) %>%
  summarise(PERM_ACRES_Rapanos = sum(IMPACT_AUTH_ACRES_PERM_PM_Rapanos_10Yr),
            PERM_LF_ACRE_CONV_Rapanos = sum(IMPACT_AUTH_LF_PERM_ACRE_CONV_PM_Rapanos_10Yr),
            PERM_TOTAL_Rapanos = sum(IMPACT_AUTH_PERM_TOTAL_PM_Rapanos_10Yr),
            TEMP_ACRES_Rapanos = sum(IMPACT_AUTH_ACRES_TEMP_TM_Rapanos_10Yr),
            TEMP_LF_ACRE_CONV_Rapanos = sum(IMPACT_AUTH_LF_TEMP_ACRE_CONV_TM_Rapanos_10Yr),
            TEMP_TOTAL_Rapanos = sum(IMPACT_AUTH_TEMP_TOTAL_TM_Rapanos_10Yr),
            PERM_ACRES_NWPR = sum(IMPACT_AUTH_ACRES_PERM_PM_NWPR_10Yr),
            PERM_LF_ACRE_CONV_NWPR = sum(IMPACT_AUTH_LF_PERM_ACRE_CONV_PM_NWPR_10Yr),
            PERM_TOTAL_NWPR = sum(IMPACT_AUTH_PERM_TOTAL_PM_NWPR_10Yr),
            TEMP_ACRES_NWPR = sum(IMPACT_AUTH_ACRES_TEMP_TM_NWPR_10Yr),
            TEMP_LF_ACRE_CONV_NWPR = sum(IMPACT_AUTH_LF_TEMP_ACRE_CONV_TM_NWPR_10Yr),
            TEMP_TOTAL_NWPR = sum(IMPACT_AUTH_TEMP_TOTAL_TM_NWPR_10Yr)) %>%
  mutate(PERM_ACRES_CHANGE = PERM_ACRES_Rapanos - PERM_ACRES_NWPR,
         PERM_LF_ACRE_CONV_CHANGE = PERM_LF_ACRE_CONV_Rapanos - PERM_LF_ACRE_CONV_NWPR,
         PERM_TOTAL_CHANGE = PERM_TOTAL_Rapanos - PERM_TOTAL_NWPR,
         TEMP_ACRES_CHANGE = TEMP_ACRES_Rapanos - TEMP_ACRES_NWPR,
         TEMP_LF_ACRE_CONV_CHANGE = TEMP_LF_ACRE_CONV_Rapanos - TEMP_LF_ACRE_CONV_NWPR,
         TEMP_TOTAL_CHANGE = TEMP_TOTAL_Rapanos - TEMP_TOTAL_NWPR)
write.csv(Impacts_HUC4,paste(OUT_DIR,"Impacts_By_HUC4_Federalism.csv",sep="/"),row.names=FALSE)

# Aggregate permanent & temporary impacts at the STATE_FINAL level
Impacts_STATE <- ORM_Dataset_Final_NoCWRNoMit %>%  
  group_by(STATE_FINAL) %>%
  summarise(PERM_ACRES_Rapanos = sum(IMPACT_AUTH_ACRES_PERM_PM_Rapanos_10Yr),
            PERM_LF_ACRE_CONV_Rapanos = sum(IMPACT_AUTH_LF_PERM_ACRE_CONV_PM_Rapanos_10Yr),
            PERM_TOTAL_Rapanos = sum(IMPACT_AUTH_PERM_TOTAL_PM_Rapanos_10Yr),
            TEMP_ACRES_Rapanos = sum(IMPACT_AUTH_ACRES_TEMP_TM_Rapanos_10Yr),
            TEMP_LF_ACRE_CONV_Rapanos = sum(IMPACT_AUTH_LF_TEMP_ACRE_CONV_TM_Rapanos_10Yr),
            TEMP_TOTAL_Rapanos = sum(IMPACT_AUTH_TEMP_TOTAL_TM_Rapanos_10Yr),
            PERM_ACRES_NWPR = sum(IMPACT_AUTH_ACRES_PERM_PM_NWPR_10Yr),
            PERM_LF_ACRE_CONV_NWPR = sum(IMPACT_AUTH_LF_PERM_ACRE_CONV_PM_NWPR_10Yr),
            PERM_TOTAL_NWPR = sum(IMPACT_AUTH_PERM_TOTAL_PM_NWPR_10Yr),
            TEMP_ACRES_NWPR = sum(IMPACT_AUTH_ACRES_TEMP_TM_NWPR_10Yr),
            TEMP_LF_ACRE_CONV_NWPR = sum(IMPACT_AUTH_LF_TEMP_ACRE_CONV_TM_NWPR_10Yr),
            TEMP_TOTAL_NWPR = sum(IMPACT_AUTH_TEMP_TOTAL_TM_NWPR_10Yr)) %>%
  mutate(PERM_ACRES_CHANGE = PERM_ACRES_Rapanos - PERM_ACRES_NWPR,
         PERM_LF_ACRE_CONV_CHANGE = PERM_LF_ACRE_CONV_Rapanos - PERM_LF_ACRE_CONV_NWPR,
         PERM_TOTAL_CHANGE = PERM_TOTAL_Rapanos - PERM_TOTAL_NWPR,
         TEMP_ACRES_CHANGE = TEMP_ACRES_Rapanos - TEMP_ACRES_NWPR,
         TEMP_LF_ACRE_CONV_CHANGE = TEMP_LF_ACRE_CONV_Rapanos - TEMP_LF_ACRE_CONV_NWPR,
         TEMP_TOTAL_CHANGE = TEMP_TOTAL_Rapanos - TEMP_TOTAL_NWPR)
write.csv(Impacts_STATE,paste(OUT_DIR,"Impacts_By_STATE_Federalism.csv",sep="/"),row.names=FALSE)

# Aggregate permanent & temporary impacts at the Cowardin Code level
Impacts_Cowardin <- ORM_Dataset_Final_NoCWRNoMit %>%  
  group_by(COWARDIN_FINAL) %>%
  summarise(PERM_ACRES_Rapanos = sum(IMPACT_AUTH_ACRES_PERM_PM_Rapanos_10Yr),
            PERM_LF_ACRE_CONV_Rapanos = sum(IMPACT_AUTH_LF_PERM_ACRE_CONV_PM_Rapanos_10Yr),
            PERM_TOTAL_Rapanos = sum(IMPACT_AUTH_PERM_TOTAL_PM_Rapanos_10Yr),
            TEMP_ACRES_Rapanos = sum(IMPACT_AUTH_ACRES_TEMP_TM_Rapanos_10Yr),
            TEMP_LF_ACRE_CONV_Rapanos = sum(IMPACT_AUTH_LF_TEMP_ACRE_CONV_TM_Rapanos_10Yr),
            TEMP_TOTAL_Rapanos = sum(IMPACT_AUTH_TEMP_TOTAL_TM_Rapanos_10Yr),
            PERM_ACRES_NWPR = sum(IMPACT_AUTH_ACRES_PERM_PM_NWPR_10Yr),
            PERM_LF_ACRE_CONV_NWPR = sum(IMPACT_AUTH_LF_PERM_ACRE_CONV_PM_NWPR_10Yr),
            PERM_TOTAL_NWPR = sum(IMPACT_AUTH_PERM_TOTAL_PM_NWPR_10Yr),
            TEMP_ACRES_NWPR = sum(IMPACT_AUTH_ACRES_TEMP_TM_NWPR_10Yr),
            TEMP_LF_ACRE_CONV_NWPR = sum(IMPACT_AUTH_LF_TEMP_ACRE_CONV_TM_NWPR_10Yr),
            TEMP_TOTAL_NWPR = sum(IMPACT_AUTH_TEMP_TOTAL_TM_NWPR_10Yr)) %>%
  mutate(PERM_ACRES_CHANGE = PERM_ACRES_Rapanos - PERM_ACRES_NWPR,
         PERM_LF_ACRE_CONV_CHANGE = PERM_LF_ACRE_CONV_Rapanos - PERM_LF_ACRE_CONV_NWPR,
         PERM_TOTAL_CHANGE = PERM_TOTAL_Rapanos - PERM_TOTAL_NWPR,
         TEMP_ACRES_CHANGE = TEMP_ACRES_Rapanos - TEMP_ACRES_NWPR,
         TEMP_LF_ACRE_CONV_CHANGE = TEMP_LF_ACRE_CONV_Rapanos - TEMP_LF_ACRE_CONV_NWPR,
         TEMP_TOTAL_CHANGE = TEMP_TOTAL_Rapanos - TEMP_TOTAL_NWPR)
write.csv(Impacts_Cowardin,paste(OUT_DIR,"Impacts_By_Cowardin_Federalism.csv",sep="/"),row.names=FALSE)

#################################################################################
#Calculate count of impacted waters not covered by a permit under Rapanos & NWPR
#################################################################################

# Aggregate permanent & temporary impacts at the HUC12 level
ImpactedWaters_HUC12 <- ORM_Dataset_Final_NoCWRNoMit %>%
  group_by(HUC12) %>%
  summarise(Waters_Count_PERM_Rapanos = sum(PERM_Count_Rapanos_10Yr),
            Waters_Count_PERM_NWPR = sum(PERM_Count_NWPR_10Yr),
            Waters_Count_PERM_Change = sum(PERM_Count_Change_10Yr),
            Waters_Count_TEMP_Rapanos = sum(TEMP_Count_Rapanos_10Yr),
            Waters_Count_TEMP_NWPR = sum(TEMP_Count_NWPR_10Yr),
            Waters_Count_TEMP_Change = sum(TEMP_Count_Change_10Yr),
            Waters_Count_Rapanos = sum(ImpactedWater_Count_Rapanos_10Yr),
            Waters_Count_NWPR = sum(ImpactedWater_Count_NWPR_10Yr),
            Waters_Count_Change = sum(ImpactedWater_Count_Change_10Yr))
write.csv(ImpactedWaters_HUC12,paste(OUT_DIR,"ImpactedWaters_By_HUC12_Federalism.csv",sep="/"),row.names=FALSE)

# Aggregate permanent & temporary impacts at the HUC4 level
ImpactedWaters_HUC4 <- ORM_Dataset_Final_NoCWRNoMit %>%
  group_by(HUC4) %>%
  summarise(Waters_Count_PERM_Rapanos = sum(PERM_Count_Rapanos_10Yr),
            Waters_Count_PERM_NWPR = sum(PERM_Count_NWPR_10Yr),
            Waters_Count_PERM_Change = sum(PERM_Count_Change_10Yr),
            Waters_Count_TEMP_Rapanos = sum(TEMP_Count_Rapanos_10Yr),
            Waters_Count_TEMP_NWPR = sum(TEMP_Count_NWPR_10Yr),
            Waters_Count_TEMP_Change = sum(TEMP_Count_Change_10Yr),
            Waters_Count_Rapanos = sum(ImpactedWater_Count_Rapanos_10Yr),
            Waters_Count_NWPR = sum(ImpactedWater_Count_NWPR_10Yr),
            Waters_Count_Change = sum(ImpactedWater_Count_Change_10Yr))
write.csv(ImpactedWaters_HUC4,paste(OUT_DIR,"ImpactedWaters_By_HUC4_Federalism.csv",sep="/"),row.names=FALSE)

# Aggregate permanent & temporary impacts at the STATE_FINAL level
ImpactedWaters_STATE <- ORM_Dataset_Final_NoCWRNoMit %>%
  group_by(STATE_FINAL) %>%
  summarise(Waters_Count_PERM_Rapanos = sum(PERM_Count_Rapanos_10Yr),
            Waters_Count_PERM_NWPR = sum(PERM_Count_NWPR_10Yr),
            Waters_Count_PERM_Change = sum(PERM_Count_Change_10Yr),
            Waters_Count_TEMP_Rapanos = sum(TEMP_Count_Rapanos_10Yr),
            Waters_Count_TEMP_NWPR = sum(TEMP_Count_NWPR_10Yr),
            Waters_Count_TEMP_Change = sum(TEMP_Count_Change_10Yr),
            Waters_Count_Rapanos = sum(ImpactedWater_Count_Rapanos_10Yr),
            Waters_Count_NWPR = sum(ImpactedWater_Count_NWPR_10Yr),
            Waters_Count_Change = sum(ImpactedWater_Count_Change_10Yr))
write.csv(ImpactedWaters_STATE,paste(OUT_DIR,"ImpactedWaters_By_STATE_Federalism.csv",sep="/"),row.names=FALSE)

# Aggregate permanent & temporary impacts at the Cowardin Code level
ImpactedWaters_Cowardin <- ORM_Dataset_Final_NoCWRNoMit %>%
  group_by(COWARDIN_FINAL) %>%
  summarise(Waters_Count_PERM_Rapanos = sum(PERM_Count_Rapanos_10Yr),
            Waters_Count_PERM_NWPR = sum(PERM_Count_NWPR_10Yr),
            Waters_Count_PERM_Change = sum(PERM_Count_Change_10Yr),
            Waters_Count_TEMP_Rapanos = sum(TEMP_Count_Rapanos_10Yr),
            Waters_Count_TEMP_NWPR = sum(TEMP_Count_NWPR_10Yr),
            Waters_Count_TEMP_Change = sum(TEMP_Count_Change_10Yr),
            Waters_Count_Rapanos = sum(ImpactedWater_Count_Rapanos_10Yr),
            Waters_Count_NWPR = sum(ImpactedWater_Count_NWPR_10Yr),
            Waters_Count_Change = sum(ImpactedWater_Count_Change_10Yr))
write.csv(ImpactedWaters_Cowardin,paste(OUT_DIR,"ImpactedWaters_By_Cowardin_Federalism.csv",sep="/"),row.names=FALSE)

####################################
#Calculate changes in permit counts
####################################

# Multiply Permit_Count by Cowardin Code-level (Rapanos_Perc - NWPR_Perc) adjustment (to calculate change in number of permits issued between the two scenarios)
# Multiply Permit_Count by SI_Perc and SG_Perc (calculate average number of Ind/Gen permits affecting one water at the HUC4/Cowardin Code level)
# Divide SI/SG- & Rapanos/NWPR-corrected Permit_Count by proportion of 10-year period that state followed Rapanos to get annual average
ORM_Dataset_Final_NoCWRNoMit <- ORM_Dataset_Final_NoCWRNoMit %>%  
  mutate(Avg_SI_Count = (Permit_Count*SI_Perc*(Rapanos_Perc_Permitted-NWPR_Perc_Permitted))/Prop_10YrPeriod_NonCWR,
         Avg_SG_Count = (Permit_Count*SG_Perc*(Rapanos_Perc_Permitted-NWPR_Perc_Permitted))/Prop_10YrPeriod_NonCWR)

write.csv(ORM_Dataset_Final_NoCWRNoMit,paste(QA_DIR,"ORM_Dataset_Final_Federalism.csv",sep="/"),row.names=FALSE)

# Subset dataset to retain permits affecting single waters only (permits affecting multiple waters will still need to be permitted)
Permit_Counts  <- ORM_Dataset_Final_NoCWRNoMit[ which(ORM_Dataset_Final_NoCWRNoMit$SingleWater==1), ]

# Subset dataset, count number of permits per HUC4/COWARDIN and HUC4/COWARDIN/STATE combo, merge with permit counts
Permit_Counts <- Permit_Counts[c("HUC4_COWARDIN","HUC4","STATE_FINAL","Avg_SI_Count","Avg_SG_Count")]
Permit_Counts$HUC4_COWARDIN_STATE <- paste(Permit_Counts$HUC4_COWARDIN,"_",Permit_Counts$STATE_FINAL)

Permit_Counts_HUC4_COWARDIN <- Permit_Counts %>%  
  group_by(HUC4_COWARDIN) %>%
  summarise(Permit_Count_HUC4_COWARDIN = n())

Permit_Counts_HUC4_COWARDIN_STATE <- Permit_Counts %>%  
  group_by(HUC4_COWARDIN_STATE) %>%
  summarise(Permit_Count_HUC4_COWARDIN_STATE = n())

Permit_Counts <- merge(x=Permit_Counts,y=Permit_Counts_HUC4_COWARDIN,by="HUC4_COWARDIN",all.x=TRUE)
Permit_Counts <- merge(x=Permit_Counts,y=Permit_Counts_HUC4_COWARDIN_STATE,by="HUC4_COWARDIN_STATE",all.x=TRUE)

# Weight adjust permit counts based on HUC4/COWARDIN/STATE divided by HUC4/COWARDIN
Permit_Counts <- Permit_Counts %>%  
  mutate(Avg_SI_Count_Adj = Avg_SI_Count * (Permit_Count_HUC4_COWARDIN_STATE/Permit_Count_HUC4_COWARDIN),
         Avg_SG_Count_Adj = Avg_SG_Count * (Permit_Count_HUC4_COWARDIN_STATE/Permit_Count_HUC4_COWARDIN))

# Subset by unique HUC4/COWARDIN/STATE
Permit_Counts2 <- Permit_Counts[!duplicated(Permit_Counts[,c("HUC4_COWARDIN_STATE")]),]
write.csv(Permit_Counts2,paste(QA_DIR,"PermitCounts_QA_Federalism.csv",sep="/"),row.names=FALSE)
rm(list = c('Permit_Counts','Permit_Counts_HUC4_COWARDIN','Permit_Counts_HUC4_COWARDIN_STATE'))

# Calculate change in average number of permits at the HUC4 level
Permits_HUC4 <- Permit_Counts2 %>%  
  group_by(HUC4) %>%
  summarise(Avg_SI_Count = sum(Avg_SI_Count_Adj),
            Avg_SG_Count = sum(Avg_SG_Count_Adj))
write.csv(Permits_HUC4,paste(OUT_DIR,"PermitCounts_By_HUC4_Federalism.csv",sep="/"),row.names=FALSE)

# Calculate change in average number of permits at the state level
Permits_State <- Permit_Counts2 %>%  
  group_by(STATE_FINAL) %>%
  summarise(Avg_SI_Count = sum(Avg_SI_Count_Adj),
            Avg_SG_Count = sum(Avg_SG_Count_Adj))
write.csv(Permits_State,paste(OUT_DIR,"PermitCounts_By_State_Federalism.csv",sep="/"),row.names=FALSE)