# ----- Code Purpose and Authors -----
# Lower Delaware River WQS Watershed Comparison
#
# Authors: Mary Sluder
# Purpose: To compare the Lower Delaware River watershed to watersheds used
# in the Netusil property values paper to inform benefit transfer.

# Date: 5/26/2023 | Last Updated: 6/28/2023

# Part 1: Set up --------------------------------------------------------------

# Clear global environment and memory
rm(list=ls())
gc()

# Install and load required packages

pacman::p_load(sf,dplyr,openxlsx,terra,rgdal,readr,rgeos,sp,tidycensus,tidyr,purrr)

# Directory Setup

var <- "59229"     ## Change var to user ID

InPath <-  paste0("C:/Users/",var,"/ICF/CWA Anniversary & Econ Support - TD3_Delaware River WQS/GIS/Data/HAWQS watersheds")

LandPath <- paste0("C:/Users/",var,"/OneDrive - ICF/Documents/LDR WQS Shapefiles")

OutPath <- paste0("C:/Users/",var,"/ICF/CWA Anniversary & Econ Support - TD3_Delaware River WQS/Property Values/Watershed Comparison/Outputs")

CensusPath <- paste0("C:/Users/",var,"/ICF/CWA Anniversary & Econ Support - TD3_Delaware River WQS/Property Values/Watershed Comparison/Inputs/ACS/")

# Part 2: Data Set Up ----------------------------------------------------------------

# Loading Data

# Land Cover (2011 NLCD)

NLCD <- rast(paste0(LandPath, "/NLCD 2011/nlcd_2011_land_cover_l48_20210604.img"))

# Watersheds

setwd(paste0(InPath, "/Burnt Bridge Creek, Washington/Watershed Shapefile from study authors/"))

BurntBridge <-  read_sf(dsn = ".", layer = "BBC whole")

BurntBridge_proj <- st_transform(BurntBridge,crs=st_crs(crs(NLCD)))

BurntBridge_proj$Area_Sqm_WS <- st_area(BurntBridge_proj)

BB_WSArea <- BurntBridge_proj %>%
  reframe(Area_Sqkm_WS = Area_Sqm_WS/1000000,
          Area_Sqkm_WS_sum = sum(Area_Sqkm_WS)) %>%
  distinct(Area_Sqkm_WS_sum) %>%
  mutate(Watershed = "BurntBridge")

setwd(paste0(InPath, "/Johnson Creek, Oregon/Watershed Shapefile from study authors/"))

JohnsonCreek <-  read_sf(dsn = ".", layer = "JohnsonCreek")

JohnsonCreek_proj <- st_transform(JohnsonCreek,crs=st_crs(crs(NLCD)))

JohnsonCreek_proj$Area_Sqm_WS <- st_area(JohnsonCreek_proj)

JC_WSArea <- JohnsonCreek_proj %>%
  reframe(Area_Sqkm_WS = Area_Sqm_WS/1000000,
          Area_Sqkm_WS_sum = sum(Area_Sqkm_WS)) %>%
  distinct(Area_Sqkm_WS_sum) %>%
  mutate(Watershed = "JohnsonCreek")

setwd(paste0(InPath, "/Lower Delaware River/"))

LowerDelaware <-  read_sf(dsn = ".", layer = "FMA_Intersecting_HUC12s")

LowerDelaware_proj <- st_transform(LowerDelaware,crs=st_crs(crs(NLCD)))

LowerDelaware_proj$Area_Sqm_WS <- st_area(LowerDelaware_proj)

LDR_WSArea <- LowerDelaware_proj %>%
  reframe(Area_Sqkm_WS = Area_Sqm_WS/1000000,
            Area_Sqkm_WS_sum = sum(Area_Sqkm_WS)) %>%
  distinct(Area_Sqkm_WS_sum) %>%
  mutate(Watershed = "LowerDelawareRiver")

LDR_WSArea_byHUC12 <- LowerDelaware_proj %>%
  group_by(HUC12) %>%
  summarise(Area_Sqkm_WS_byHUC12 = Area_Sqm_WS/1000000) %>%
  mutate(Watershed = "LowerDelawareRiver") %>%
  st_drop_geometry()

WS_area_list <- list(LDR_WSArea,BB_WSArea,JC_WSArea)

WS_area <- do.call(rbind,WS_area_list)

# Individual HUC12s

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020401/"))

LDR_HUC12_1 <- read_sf(dsn = ".", layer = "020402020401")

LDR_HUC12_1_proj <- st_transform(LDR_HUC12_1,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020402/"))

LDR_HUC12_2 <- read_sf(dsn = ".", layer = "020402020402")

LDR_HUC12_2_proj <- st_transform(LDR_HUC12_2,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020403/"))

LDR_HUC12_3 <- read_sf(dsn = ".", layer = "020402020403")

LDR_HUC12_3_proj <- st_transform(LDR_HUC12_3,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020404/"))

LDR_HUC12_4 <- read_sf(dsn = ".", layer = "020402020404")

LDR_HUC12_4_proj <- st_transform(LDR_HUC12_4,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020405/"))

LDR_HUC12_5 <- read_sf(dsn = ".", layer = "020402020405")

LDR_HUC12_5_proj <- st_transform(LDR_HUC12_5,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020503/"))

LDR_HUC12_6 <- read_sf(dsn = ".", layer = "020402020503")

LDR_HUC12_6_proj <- st_transform(LDR_HUC12_6,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020505/"))

LDR_HUC12_7 <- read_sf(dsn = ".", layer = "020402020505")

LDR_HUC12_7_proj <- st_transform(LDR_HUC12_7,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020506/"))

LDR_HUC12_8 <- read_sf(dsn = ".", layer = "020402020506")

LDR_HUC12_8_proj <- st_transform(LDR_HUC12_8,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020507/"))

LDR_HUC12_9 <- read_sf(dsn = ".", layer = "020402020507")

LDR_HUC12_9_proj <- st_transform(LDR_HUC12_9,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020606/"))

LDR_HUC12_10 <- read_sf(dsn = ".", layer = "020402020606")

LDR_HUC12_10_proj <- st_transform(LDR_HUC12_10,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020607/"))

LDR_HUC12_11 <- read_sf(dsn = ".", layer = "020402020607")

LDR_HUC12_11_proj <- st_transform(LDR_HUC12_11,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402020608/"))

LDR_HUC12_12 <- read_sf(dsn = ".", layer = "020402020608")

LDR_HUC12_12_proj <- st_transform(LDR_HUC12_12,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402031008/"))

LDR_HUC12_13 <- read_sf(dsn = ".", layer = "020402031008")

LDR_HUC12_13_proj <- st_transform(LDR_HUC12_13,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402040000/"))

LDR_HUC12_14 <- read_sf(dsn = ".", layer = "020402040000")

LDR_HUC12_14_proj <- st_transform(LDR_HUC12_14,crs=st_crs(crs(NLCD)))

setwd(paste0(InPath, "/Lower Delaware River/Individual HUC12s/020402050505/"))

LDR_HUC12_15 <- read_sf(dsn = ".", layer = "020402050505")

LDR_HUC12_15_proj <- st_transform(LDR_HUC12_15,crs=st_crs(crs(NLCD)))

# Census Data Prep

# Activate API Key

census_api_key("0d3711e1f7f03d54ea326708428a15bcec0c0621") # Key registered to Sam Ennett

census_year <- c(2000)

ACS_2010_var <- c("B01001_001", "B19013_001E", "DP04_0001E", "DP04_0088E","DP04_0002E","DP04_0007E","DP04_0008E","DP04_0009E","DP04_0010E","DP04_0011E",
             "DP04_0012E","DP04_0013E","DP04_0014E","DP04_0015E")
# Population, Median Household Income, Number of Households, Median Housing Value,Number occupied housing units,
# 1 unit-detached, 1 unit-attached, 2 units, 3 or 4 units, 5-9 units, 10-29 units, 20 or more units, mobile home, boat/RV/van

ACS_2021_var <- c("B01001_001", "B19013_001E", "DP04_0001E", "DP04_0089E","DP04_0002E","DP04_0007E","DP04_0008E","DP04_0009E","DP04_0010E","DP04_0011E",
                  "DP04_0012E","DP04_0013E","DP04_0014E","DP04_0015E")
# Population, Median Household Income, Number of Households, Median Housing Value,Number occupied housing units,
# 1 unit-detached, 1 unit-attached, 2 units, 3 or 4 units, 5-9 units, 10-29 units, 20 or more units, mobile home, boat/RV/van

Census_var4 <- c("DP4_C2","DP4_C4","DP4_C6","DP4_C8","DP4_C10","DP4_C12","DP4_C14","DP4_C16","DP4_C18","DP4_C125","DP1_C0")
# Number 1 unit detached, Number 1 unit attached, Number 2 units, Number 3 or 4 units,
# Number 5-9 units, Number 10 to 19, Number 20 or more, Number mobile home, number boat, rv, van, etc., Median owner-occupied housing value, total population
Census_var3 <- c("DP3_C112","DP4_C55","DP4_C0") # median household income, Number of occupied housing units, Number total households

# Note DP1 and DP3 cannot mix; throws error

# Pull 2000 census tracts

Census_2000_df4 <- get_decennial(geography = "tract",
                  year = census_year,
                  sumfile = c("sf4profile"),
                  variables = Census_var4,
                  state = c("PA", "NJ", "DE", "WA", "OR"),
                  geometry = T)

Census_2000_df4 <- Census_2000_df4 %>%
  st_transform(crs=st_crs(crs(NLCD)))

Census_2000_df3 <- get_decennial(geography = "tract",
                         year = census_year,
                         sumfile = c("sf3profile"),
                         variables = Census_var3,
                         state = c("PA", "NJ", "DE", "WA", "OR"),
                         geometry = T)

Census_2000_df3 <- Census_2000_df3 %>%
  st_transform(crs=st_crs(crs(NLCD)))

Census_2000_df <- rbind(Census_2000_df3,Census_2000_df4) %>%
  pivot_wider(names_from = "variable", values_from = "value") %>%
  dplyr::select(-NAME)  %>%
  rename("Med_Inc" = "DP3_C112",
         "Occupied_Housing" = "DP4_C55",
         "Med_Value" = "DP4_C125",
         "Population" = "DP1_C0",
         "Households" = "DP4_C0",
         "One_Detached" = "DP4_C2",
         "One_Attached" = "DP4_C4",
         "Two_units" = "DP4_C6",
         "Three_or_Four" = "DP4_C8",
         "Five_to_Nine" = "DP4_C10",
         "Ten_to_Nineteen" = "DP4_C12",
         "Twenty_or_More" = "DP4_C14",
         "Mobile_home" = "DP4_C16",
         "Boat_RV" = "DP4_C18")

Census_2000_df_centroid <- Census_2000_df %>%
  st_centroid(Census_2000_df) %>%
  st_zm()

# Pull 2007 ACS tracts

# ACS_2010_year <- c(2010)
#
# ACS_2010_df <- get_acs(geography = "census tract",
#                   year = ACS_2010_year,
#                   survey = "acs5",
#                   state = c("PA", "NJ", "DE", "WA", "OR"),
#                   variables = ACS_var,
#                   geometry = T) %>%
#   pivot_wider(names_from = "variable", values_from = "estimate")
#
# ACS_2021_year <- c(2021)
#
# ACS_2021_df <- get_acs(geography = "census tract",
#                        year = ACS_2021_year,
#                        survey = "acs5",
#                        state = c("PA", "NJ", "DE", "WA", "OR"),
#                        variables = ACS_var,
#                        geometry = T) %>%
#   pivot_wider(names_from = "variable", values_from = "estimate")

# Since API for ACS is not working; pulling in shapefiles & datasets manually
# @MS Note: Need "tract" not "census tract"

setwd(paste0(CensusPath,"Shapefiles/Oregon/2010/"))

OR_2010 <- read_sf(dsn = ".", layer = "tl_2010_41_tract10")

OR_2010_proj <- st_transform(OR_2010,crs=st_crs(crs(NLCD)))

setwd(paste0(CensusPath,"Shapefiles/Oregon/2021/"))

OR_2021 <- read_sf(dsn = ".", layer = "tl_2021_41_tract")

OR_2021_proj <- st_transform(OR_2021,crs=st_crs(crs(NLCD)))

setwd(paste0(CensusPath,"Shapefiles/Washington/2021/"))

WA_2021 <- read_sf(dsn = ".", layer = "tl_2021_53_tract")

WA_2021_proj <- st_transform(WA_2021,crs=st_crs(crs(NLCD)))

setwd(paste0(CensusPath,"Shapefiles/Washington/2010/"))

WA_2010 <- read_sf(dsn = ".", layer = "tl_2010_53_tract10")

WA_2010_proj <- st_transform(WA_2010,crs=st_crs(crs(NLCD)))

setwd(paste0(CensusPath,"Shapefiles/Pennsylvania/2021/"))

PA_2021 <- read_sf(dsn = ".", layer = "tl_2021_42_tract")

PA_2021_proj <- st_transform(PA_2021,crs=st_crs(crs(NLCD)))

setwd(paste0(CensusPath,"Shapefiles/Pennsylvania/2010/"))

PA_2010 <- read_sf(dsn = ".", layer = "tl_2010_42_tract10")

PA_2010_proj <- st_transform(PA_2010,crs=st_crs(crs(NLCD)))

setwd(paste0(CensusPath,"Shapefiles/Delaware/2021/"))

DE_2021 <- read_sf(dsn = ".", layer = "tl_2021_10_tract")

DE_2021_proj <- st_transform(DE_2021,crs=st_crs(crs(NLCD)))

setwd(paste0(CensusPath,"Shapefiles/Delaware/2010/"))

DE_2010 <- read_sf(dsn = ".", layer = "tl_2010_10_tract10")

DE_2010_proj <- st_transform(DE_2010,crs=st_crs(crs(NLCD)))

setwd(paste0(CensusPath,"Shapefiles/New Jersey/2021/"))

NJ_2021 <- read_sf(dsn = ".", layer = "tl_2021_34_tract")

NJ_2021_proj <- st_transform(NJ_2021,crs=st_crs(crs(NLCD)))

setwd(paste0(CensusPath,"Shapefiles/New Jersey/2010/"))

NJ_2010 <- read_sf(dsn = ".", layer = "tl_2010_34_tract10")

NJ_2010_proj <- st_transform(NJ_2010,crs=st_crs(crs(NLCD)))

House_Characteristics_2010 <- read.csv(paste0(CensusPath, "Household Characteristics/2010.csv")) %>%
  select(GEO_ID,DP04_0001E,DP04_0088E,DP04_0002E,DP04_0007E,DP04_0008E,DP04_0009E,DP04_0010E,DP04_0011E,DP04_0012E,DP04_0013E,DP04_0014E,DP04_0015E)

House_Characteristics_2021 <- read.csv(paste0(CensusPath, "Household Characteristics/2021.csv")) %>%
  select(GEO_ID,DP04_0001E,DP04_0089E,DP04_0002E,DP04_0007E,DP04_0008E,DP04_0009E,DP04_0010E,DP04_0011E,DP04_0012E,DP04_0013E,DP04_0014E,DP04_0015E)

House_Inc_2010 <- read.csv(paste0(CensusPath, "Household Income/2010.csv")) %>%
  select(GEO_ID,B19013_001E)

House_Inc_2021 <- read.csv(paste0(CensusPath, "Household Income/2021.csv")) %>%
  select(GEO_ID,B19013_001E)

Population_2010 <- read.csv(paste0(CensusPath, "Population/2010.csv")) %>%
  select(GEO_ID,B01001_001E)

Population_2021 <- read.csv(paste0(CensusPath, "Population/2021.csv")) %>%
  select(GEO_ID,B01001_001E)

Shapefile_list_2010 <- list(PA_2010_proj,WA_2010_proj,NJ_2010_proj,DE_2010_proj,OR_2010_proj)

All_States_2010 <- do.call(rbind, Shapefile_list_2010) %>%
  select(GEOID10,geometry)

Shapefile_list_2021 <- list(PA_2021_proj,WA_2021_proj,NJ_2021_proj,DE_2021_proj,OR_2021_proj)

All_States_2021 <- do.call(rbind, Shapefile_list_2021) %>%
  select(GEOID,geometry)

# Part 3: Land Cover Analysis -------------------------------------------------------------

# All HUCs

# Intersect NLCD dataset with watersheds

BB_landcover <- terra::extract(NLCD,BurntBridge_proj, df = T)

JC_landcover <- terra::extract(NLCD,JohnsonCreek_proj, df = T)

LDR_landcover <- terra::extract(NLCD,LowerDelaware_proj, df = T)

# Rename

names(BB_landcover) <- c("ID","NLCD_Land_Cover_Class")

names(JC_landcover) <- c("ID","NLCD_Land_Cover_Class")

names(LDR_landcover) <- c("ID","NLCD_Land_Cover_Class")

# Find percent land coverage

watersheds <- list(BB_landcover,JC_landcover,LDR_landcover)

names(watersheds) <- c("BurntBridge","JohnsonCreek","LowerDelawareRiver")

for (i in seq_along(watersheds)) {

  LU_watersheds <- watersheds[[i]]

LU_tot <- LU_watersheds %>%
  summarise(
    LU_tot = n()
  )

LU_tot_NoWater <- LU_watersheds %>%
  filter(!c(NLCD_Land_Cover_Class == "Open Water" | NLCD_Land_Cover_Class == "Perennial Ice/Snow" |
           NLCD_Land_Cover_Class == "Woody Wetlands" | NLCD_Land_Cover_Class == "Emergent Herbaceous Wetlands")) %>%
  summarise(
    LU_tot_NoWater = n()
  )

# Water

LU_Water <- LU_watersheds %>%
  filter(NLCD_Land_Cover_Class == "Open Water" | NLCD_Land_Cover_Class == "Perennial Ice/Snow") %>%
  summarise(
    Water = n()
  )

Landuse_Water <- merge(LU_Water,LU_tot,all=T)

Landuse_Water[is.na(Landuse_Water)] = 0

# Developed

LU_Developed <- LU_watersheds %>%
  filter(NLCD_Land_Cover_Class == "Developed, Open Space" | NLCD_Land_Cover_Class == "Developed, Low Intensity"
         | NLCD_Land_Cover_Class == "Developed, Medium Intensity" | NLCD_Land_Cover_Class == "Developed, High Intensity") %>%
  summarise(
    Developed = n()
  )

Landuse_Developed <- merge(LU_Developed,LU_tot,all=T)

Landuse_Developed[is.na(Landuse_Developed)] = 0

Landuse_Developed_NoWater <- merge(LU_Developed,LU_tot_NoWater,all=T)

Landuse_Developed_NoWater[is.na(Landuse_Developed_NoWater)] = 0

# Barren

LU_Barren <- LU_watersheds %>%
  filter(NLCD_Land_Cover_Class == "Barren Land (Rock/Sand/Clay)") %>%
  summarise(
    Barren = n()
  )

Landuse_Barren <- merge(LU_Barren,LU_tot,all=T)

Landuse_Barren[is.na(Landuse_Barren)] = 0

Landuse_Barren_NoWater <- merge(LU_Barren,LU_tot_NoWater,all=T)

Landuse_Barren_NoWater[is.na(Landuse_Barren_NoWater)] = 0

# Forest

LU_Forest <- LU_watersheds %>%
  filter(NLCD_Land_Cover_Class == "Deciduous Forest" | NLCD_Land_Cover_Class == "Evergreen Forest"
         | NLCD_Land_Cover_Class == "Mixed Forest") %>%
  summarise(
    Forest = n()
  )

Landuse_Forest <- merge(LU_Forest,LU_tot,all=T)

Landuse_Forest[is.na(Landuse_Forest)] = 0

Landuse_Forest_NoWater <- merge(LU_Forest,LU_tot_NoWater,all=T)

Landuse_Forest_NoWater[is.na(Landuse_Forest_NoWater)] = 0

# Shrubland

LU_Shrubland <- LU_watersheds %>%
  filter(NLCD_Land_Cover_Class == "Dwarf Scrub" | NLCD_Land_Cover_Class == "Shrub/Scrub") %>%
  summarise(
    Shrubland = n()
  )

Landuse_Shrubland <- merge(LU_Shrubland,LU_tot,all=T)

Landuse_Shrubland[is.na(Landuse_Shrubland)] = 0

Landuse_Shrubland_NoWater <- merge(LU_Shrubland,LU_tot_NoWater,all=T)

Landuse_Shrubland_NoWater[is.na(Landuse_Shrubland_NoWater)] = 0

# Herbaceous

LU_Herbaceous <- LU_watersheds %>%
  filter(NLCD_Land_Cover_Class == "Grassland/Herbaceous" | NLCD_Land_Cover_Class == "Sedge/Herbaceous"
         | NLCD_Land_Cover_Class == "Lichens" | NLCD_Land_Cover_Class == "Moss") %>%
  summarise(
    Herbaceous = n()
  )

Landuse_Herbaceous <- merge(LU_Herbaceous,LU_tot,all=T)

Landuse_Herbaceous[is.na(Landuse_Herbaceous)] = 0

Landuse_Herbaceous_NoWater <- merge(LU_Herbaceous,LU_tot_NoWater,all=T)

Landuse_Herbaceous_NoWater[is.na(Landuse_Herbaceous_NoWater)] = 0

# Agriculture

LU_Agriculture <- LU_watersheds %>%
  filter(NLCD_Land_Cover_Class == "Hay/Pasture" | NLCD_Land_Cover_Class == "Cultivated Crops") %>%
  summarise(
    Agriculture = n()
  )

Landuse_Agriculture <- merge(LU_Agriculture,LU_tot,all=T)

Landuse_Agriculture[is.na(Landuse_Agriculture)] = 0

Landuse_Agriculture_NoWater <- merge(LU_Agriculture,LU_tot_NoWater,all=T)

Landuse_Agriculture_NoWater[is.na(Landuse_Agriculture_NoWater)] = 0

# Wetlands

LU_Wetlands <- LU_watersheds %>%
  filter(NLCD_Land_Cover_Class == "Woody Wetlands" | NLCD_Land_Cover_Class == "Emergent Herbaceous Wetlands") %>%
  summarise(
    Wetlands = n()
  )

Landuse_Wetlands <- merge(LU_Wetlands,LU_tot,all=T)

Landuse_Wetlands[is.na(Landuse_Wetlands)] = 0

# Percent coverage

LU_Water_Percent <- Landuse_Water %>%
  mutate(LU_Percent_Water = Water / LU_tot)

LU_Developed_Percent <- Landuse_Developed %>%
  mutate(LU_Percent_Developed = Developed / LU_tot)

LU_Barren_Percent <- Landuse_Barren %>%
  mutate(LU_Percent_Barren = Barren / LU_tot)

LU_Forest_Percent <- Landuse_Forest %>%
  mutate(LU_Percent_Forest = Forest / LU_tot)

LU_Shrubland_Percent <- Landuse_Shrubland %>%
  mutate(LU_Percent_Shrubland = Shrubland / LU_tot)

LU_Herbaceous_Percent <- Landuse_Herbaceous %>%
  mutate(LU_Percent_Herbaceous = Herbaceous / LU_tot)

LU_Agriculture_Percent <- Landuse_Agriculture %>%
  mutate(LU_Percent_Agriculture = Agriculture / LU_tot)

LU_Wetlands_Percent <- Landuse_Wetlands %>%
  mutate(LU_Percent_Wetlands = Wetlands / LU_tot)

# Percent coverage -  no water /wetlands

LU_Developed_Percent_NoWater <- Landuse_Developed_NoWater %>%
  mutate(LU_Percent_Developed_NoWater = Developed / LU_tot_NoWater)

LU_Barren_Percent_NoWater <- Landuse_Barren_NoWater %>%
  mutate(LU_Percent_Barren_NoWater = Barren / LU_tot_NoWater)

LU_Forest_Percent_NoWater <- Landuse_Forest_NoWater %>%
  mutate(LU_Percent_Forest_NoWater = Forest / LU_tot_NoWater)

LU_Shrubland_Percent_NoWater <- Landuse_Shrubland_NoWater %>%
  mutate(LU_Percent_Shrubland_NoWater = Shrubland / LU_tot_NoWater)

LU_Herbaceous_Percent_NoWater <- Landuse_Herbaceous_NoWater %>%
  mutate(LU_Percent_Herbaceous_NoWater = Herbaceous / LU_tot_NoWater)

LU_Agriculture_Percent_NoWater <- Landuse_Agriculture_NoWater %>%
  mutate(LU_Percent_Agriculture_NoWater = Agriculture / LU_tot_NoWater)

# Bind together

if (exists("LU_Water_Percent_bind")) {
  LU_Water_Percent_bind <- rbind(LU_Water_Percent_bind, LU_Water_Percent)
} else {
  LU_Water_Percent_bind <- LU_Water_Percent
}

if (exists("LU_Developed_Percent_bind")) {
  LU_Developed_Percent_bind <- rbind(LU_Developed_Percent_bind, LU_Developed_Percent)
} else {
  LU_Developed_Percent_bind <- LU_Developed_Percent
}

if (exists("LU_Barren_Percent_bind")) {
  LU_Barren_Percent_bind <- rbind(LU_Barren_Percent_bind, LU_Barren_Percent)
} else {
  LU_Barren_Percent_bind <- LU_Barren_Percent
}

if (exists("LU_Forest_Percent_bind")) {
  LU_Forest_Percent_bind <- rbind(LU_Forest_Percent_bind, LU_Forest_Percent)
} else {
  LU_Forest_Percent_bind <- LU_Forest_Percent
}

if (exists("LU_Shrubland_Percent_bind")) {
  LU_Shrubland_Percent_bind <- rbind(LU_Shrubland_Percent_bind, LU_Shrubland_Percent)
} else {
  LU_Shrubland_Percent_bind <- LU_Shrubland_Percent
}

if (exists("LU_Herbaceous_Percent_bind")) {
  LU_Herbaceous_Percent_bind <- rbind(LU_Herbaceous_Percent_bind, LU_Herbaceous_Percent)
} else {
  LU_Herbaceous_Percent_bind <- LU_Herbaceous_Percent
}

if (exists("LU_Agriculture_Percent_bind")) {
  LU_Agriculture_Percent_bind <- rbind(LU_Agriculture_Percent_bind, LU_Agriculture_Percent)
} else {
  LU_Agriculture_Percent_bind <- LU_Agriculture_Percent
}

if (exists("LU_Wetlands_Percent_bind")) {
  LU_Wetlands_Percent_bind <- rbind(LU_Wetlands_Percent_bind, LU_Wetlands_Percent)
} else {
  LU_Wetlands_Percent_bind <- LU_Wetlands_Percent
}

# No water/wetlands

if (exists("LU_Developed_Percent_NoWater_bind")) {
  LU_Developed_Percent_NoWater_bind <- rbind(LU_Developed_Percent_NoWater_bind, LU_Developed_Percent_NoWater)
} else {
  LU_Developed_Percent_NoWater_bind <- LU_Developed_Percent_NoWater
}

if (exists("LU_Barren_Percent_NoWater_bind")) {
  LU_Barren_Percent_NoWater_bind <- rbind(LU_Barren_Percent_NoWater_bind, LU_Barren_Percent_NoWater)
} else {
  LU_Barren_Percent_NoWater_bind <- LU_Barren_Percent_NoWater
}

if (exists("LU_Forest_Percent_NoWater_bind")) {
  LU_Forest_Percent_NoWater_bind <- rbind(LU_Forest_Percent_NoWater_bind, LU_Forest_Percent_NoWater)
} else {
  LU_Forest_Percent_NoWater_bind <- LU_Forest_Percent_NoWater
}

if (exists("LU_Shrubland_Percent_NoWater_bind")) {
  LU_Shrubland_Percent_NoWater_bind <- rbind(LU_Shrubland_Percent_NoWater_bind, LU_Shrubland_Percent_NoWater)
} else {
  LU_Shrubland_Percent_NoWater_bind <- LU_Shrubland_Percent_NoWater
}

if (exists("LU_Herbaceous_Percent_NoWater_bind")) {
  LU_Herbaceous_Percent_NoWater_bind <- rbind(LU_Herbaceous_Percent_NoWater_bind, LU_Herbaceous_Percent_NoWater)
} else {
  LU_Herbaceous_Percent_NoWater_bind <- LU_Herbaceous_Percent_NoWater
}

if (exists("LU_Agriculture_Percent_NoWater_bind")) {
  LU_Agriculture_Percent_NoWater_bind <- rbind(LU_Agriculture_Percent_NoWater_bind, LU_Agriculture_Percent_NoWater)
} else {
  LU_Agriculture_Percent_NoWater_bind <- LU_Agriculture_Percent_NoWater
}

}

# No water

Results_list_NoWater <- list(LU_Developed_Percent_NoWater_bind,LU_Barren_Percent_NoWater_bind,LU_Forest_Percent_NoWater_bind,
                     LU_Shrubland_Percent_NoWater_bind,LU_Herbaceous_Percent_NoWater_bind,LU_Agriculture_Percent_NoWater_bind)

Results_bind_NoWater <- Reduce(function(dtf1, dtf2) merge(dtf1, dtf2, all = T),
                             Results_list_NoWater) %>%
  select(LU_Percent_Developed_NoWater,LU_Percent_Barren_NoWater,LU_Percent_Forest_NoWater,
         LU_Percent_Shrubland_NoWater,LU_Percent_Herbaceous_NoWater,LU_Percent_Agriculture_NoWater)

Results_bind_NoWater$Watershed <- names(watersheds)

Results_subfinal_NoWater <- Results_bind_NoWater %>%
  select(Watershed,LU_Percent_Developed_NoWater,LU_Percent_Barren_NoWater,LU_Percent_Forest_NoWater,
         LU_Percent_Shrubland_NoWater,LU_Percent_Herbaceous_NoWater,LU_Percent_Agriculture_NoWater)

Results_final_NoWater <- merge(Results_subfinal_NoWater,WS_area,all=T)

# write.csv(Results_final_NoWater,paste0(OutPath,"/LandCoverComparison_NoWater-2023-06-26.csv"), row.names = F)

# With water

Results_list <- list(LU_Water_Percent_bind,LU_Developed_Percent_bind,LU_Barren_Percent_bind,LU_Forest_Percent_bind,
                     LU_Shrubland_Percent_bind,LU_Herbaceous_Percent_bind,LU_Agriculture_Percent_bind,LU_Wetlands_Percent_bind)

Results_bind <- do.call(cbind, Results_list)

Results_bind$Watershed <- names(watersheds)

Results_subfinal <- Results_bind %>%
  select(Watershed,LU_Percent_Water,LU_Percent_Developed,LU_Percent_Barren,LU_Percent_Forest,
         LU_Percent_Shrubland,LU_Percent_Herbaceous,LU_Percent_Agriculture,LU_Percent_Wetlands)

Results_final <- merge(Results_subfinal,WS_area,all=T)

# write.csv(Results_final,paste0(OutPath,"/LandCoverComparison-2023-06-26.csv"), row.names = F)

# By HUC12, with water

LDR_landcover$HUC12 <- LowerDelaware_proj$HUC12[LDR_landcover$ID]

LU_tot_byHUC <- LDR_landcover %>%
  group_by(HUC12) %>%
  summarise(
    LU_tot_byHUC = n()
  )

# Water

LU_Water_byHUC <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Open Water" | NLCD_Land_Cover_Class == "Perennial Ice/Snow") %>%
  group_by(HUC12) %>%
  summarise(
    Water_byHUC = n()
  )

Landuse_Water_byHUC <- merge(LU_Water_byHUC,LU_tot_byHUC,all=T)

Landuse_Water_byHUC[is.na(Landuse_Water_byHUC)] = 0

# Developed

LU_Developed_byHUC <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Developed, Open Space" | NLCD_Land_Cover_Class == "Developed, Low Intensity"
         | NLCD_Land_Cover_Class == "Developed, Medium Intensity" | NLCD_Land_Cover_Class == "Developed, High Intensity") %>%
  group_by(HUC12) %>%
  summarise(
    Developed_byHUC = n()
  )

Landuse_Developed_byHUC <- merge(LU_Developed_byHUC,LU_tot_byHUC,all=T)

Landuse_Developed_byHUC[is.na(Landuse_Developed_byHUC)] = 0

# Barren

LU_Barren_byHUC <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Barren Land (Rock/Sand/Clay)") %>%
  group_by(HUC12) %>%
  summarise(
    Barren_byHUC = n()
  )

Landuse_Barren_byHUC <- merge(LU_Barren_byHUC,LU_tot_byHUC,all=T)

Landuse_Barren_byHUC[is.na(Landuse_Barren_byHUC)] = 0

# Forest

LU_Forest_byHUC <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Deciduous Forest" | NLCD_Land_Cover_Class == "Evergreen Forest"
         | NLCD_Land_Cover_Class == "Mixed Forest") %>%
  group_by(HUC12) %>%
  summarise(
    Forest_byHUC = n()
  )

Landuse_Forest_byHUC <- merge(LU_Forest_byHUC,LU_tot_byHUC,all=T)

Landuse_Forest_byHUC[is.na(Landuse_Forest_byHUC)] = 0

# Shrubland

LU_Shrubland_byHUC <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Dwarf Scrub" | NLCD_Land_Cover_Class == "Shrub/Scrub") %>%
  group_by(HUC12) %>%
  summarise(
    Shrubland_byHUC = n()
  )

Landuse_Shrubland_byHUC <- merge(LU_Shrubland_byHUC,LU_tot_byHUC,all=T)

Landuse_Shrubland_byHUC[is.na(Landuse_Shrubland_byHUC)] = 0

# Herbaceous

LU_Herbaceous_byHUC <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Grassland/Herbaceous" | NLCD_Land_Cover_Class == "Sedge/Herbaceous"
         | NLCD_Land_Cover_Class == "Lichens" | NLCD_Land_Cover_Class == "Moss") %>%
  group_by(HUC12) %>%
  summarise(
    Herbaceous_byHUC = n()
  )

Landuse_Herbaceous_byHUC <- merge(LU_Herbaceous_byHUC,LU_tot_byHUC,all=T)

Landuse_Herbaceous_byHUC[is.na(Landuse_Herbaceous_byHUC)] = 0

# Agriculture

LU_Agriculture_byHUC <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Hay/Pasture" | NLCD_Land_Cover_Class == "Cultivated Crops") %>%
  group_by(HUC12) %>%
  summarise(
    Agriculture_byHUC = n()
  )

Landuse_Agriculture_byHUC <- merge(LU_Agriculture_byHUC,LU_tot_byHUC,all=T)

Landuse_Agriculture_byHUC[is.na(Landuse_Agriculture_byHUC)] = 0

# Wetlands

LU_Wetlands_byHUC <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Woody Wetlands" | NLCD_Land_Cover_Class == "Emergent Herbaceous Wetlands") %>%
  group_by(HUC12) %>%
  summarise(
    Wetlands_byHUC = n()
  )

Landuse_Wetlands_byHUC <- merge(LU_Wetlands_byHUC,LU_tot_byHUC,all=T)

Landuse_Wetlands_byHUC[is.na(Landuse_Wetlands_byHUC)] = 0

# Percent coverage

LU_Water_Percent_byHUC <- Landuse_Water_byHUC %>%
  mutate(LU_Percent_Water_byHUC = Water_byHUC / LU_tot_byHUC)

LU_Developed_Percent_byHUC <- Landuse_Developed_byHUC %>%
  mutate(LU_Percent_Developed_byHUC = Developed_byHUC / LU_tot_byHUC)

LU_Barren_Percent_byHUC <- Landuse_Barren_byHUC %>%
  mutate(LU_Percent_Barren_byHUC = Barren_byHUC / LU_tot_byHUC)

LU_Forest_Percent_byHUC <- Landuse_Forest_byHUC %>%
  mutate(LU_Percent_Forest_byHUC = Forest_byHUC / LU_tot_byHUC)

LU_Shrubland_Percent_byHUC <- Landuse_Shrubland_byHUC %>%
  mutate(LU_Percent_Shrubland_byHUC = Shrubland_byHUC / LU_tot_byHUC)

LU_Herbaceous_Percent_byHUC <- Landuse_Herbaceous_byHUC %>%
  mutate(LU_Percent_Herbaceous_byHUC = Herbaceous_byHUC / LU_tot_byHUC)

LU_Agriculture_Percent_byHUC <- Landuse_Agriculture_byHUC %>%
  mutate(LU_Percent_Agriculture_byHUC = Agriculture_byHUC / LU_tot_byHUC)

LU_Wetlands_Percent_byHUC <- Landuse_Wetlands_byHUC %>%
  mutate(LU_Percent_Wetlands_byHUC = Wetlands_byHUC / LU_tot_byHUC)

Results_list_byHUC <- list(LU_Water_Percent_byHUC,LU_Developed_Percent_byHUC,LU_Barren_Percent_byHUC,LU_Forest_Percent_byHUC,
                     LU_Shrubland_Percent_byHUC,LU_Herbaceous_Percent_byHUC,LU_Agriculture_Percent_byHUC,LU_Wetlands_Percent_byHUC)

Results_byHUC_bind <- Reduce(function(dtf1, dtf2) merge(dtf1, dtf2, all = T),
       Results_list_byHUC) %>%
  select(HUC12,LU_Percent_Water_byHUC,LU_Percent_Developed_byHUC,LU_Percent_Barren_byHUC,LU_Percent_Forest_byHUC,
         LU_Percent_Shrubland_byHUC,LU_Percent_Herbaceous_byHUC,LU_Percent_Agriculture_byHUC,LU_Percent_Wetlands_byHUC)

Results_byHUC_final <- merge(Results_byHUC_bind,LDR_WSArea_byHUC12,all=T)

# By HUC12, no water

LU_tot_byHUC_NoWater <- LDR_landcover %>%
  filter(!c(NLCD_Land_Cover_Class == "Open Water" | NLCD_Land_Cover_Class == "Perennial Ice/Snow" |
              NLCD_Land_Cover_Class == "Woody Wetlands" | NLCD_Land_Cover_Class == "Emergent Herbaceous Wetlands")) %>%
  group_by(HUC12) %>%
  summarise(
    LU_tot_byHUC_NoWater = n()
  )

# Developed

LU_Developed_byHUC_NoWater <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Developed, Open Space" | NLCD_Land_Cover_Class == "Developed, Low Intensity"
         | NLCD_Land_Cover_Class == "Developed, Medium Intensity" | NLCD_Land_Cover_Class == "Developed, High Intensity") %>%
  group_by(HUC12) %>%
  summarise(
    Developed_byHUC_NoWater = n()
  )

Landuse_Developed_byHUC_NoWater <- merge(LU_Developed_byHUC_NoWater,LU_tot_byHUC_NoWater,all=T)

Landuse_Developed_byHUC_NoWater[is.na(Landuse_Developed_byHUC_NoWater)] = 0

# Barren

LU_Barren_byHUC_NoWater <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Barren Land (Rock/Sand/Clay)") %>%
  group_by(HUC12) %>%
  summarise(
    Barren_byHUC_NoWater = n()
  )

Landuse_Barren_byHUC_NoWater <- merge(LU_Barren_byHUC_NoWater,LU_tot_byHUC_NoWater,all=T)

Landuse_Barren_byHUC_NoWater[is.na(Landuse_Barren_byHUC_NoWater)] = 0

# Forest

LU_Forest_byHUC_NoWater <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Deciduous Forest" | NLCD_Land_Cover_Class == "Evergreen Forest"
         | NLCD_Land_Cover_Class == "Mixed Forest") %>%
  group_by(HUC12) %>%
  summarise(
    Forest_byHUC_NoWater = n()
  )

Landuse_Forest_byHUC_NoWater <- merge(LU_Forest_byHUC_NoWater,LU_tot_byHUC_NoWater,all=T)

Landuse_Forest_byHUC_NoWater[is.na(Landuse_Forest_byHUC_NoWater)] = 0

# Shrubland

LU_Shrubland_byHUC_NoWater <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Dwarf Scrub" | NLCD_Land_Cover_Class == "Shrub/Scrub") %>%
  group_by(HUC12) %>%
  summarise(
    Shrubland_byHUC_NoWater = n()
  )

Landuse_Shrubland_byHUC_NoWater <- merge(LU_Shrubland_byHUC_NoWater,LU_tot_byHUC_NoWater,all=T)

Landuse_Shrubland_byHUC_NoWater[is.na(Landuse_Shrubland_byHUC_NoWater)] = 0

# Herbaceous

LU_Herbaceous_byHUC_NoWater <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Grassland/Herbaceous" | NLCD_Land_Cover_Class == "Sedge/Herbaceous"
         | NLCD_Land_Cover_Class == "Lichens" | NLCD_Land_Cover_Class == "Moss") %>%
  group_by(HUC12) %>%
  summarise(
    Herbaceous_byHUC_NoWater = n()
  )

Landuse_Herbaceous_byHUC_NoWater <- merge(LU_Herbaceous_byHUC_NoWater,LU_tot_byHUC_NoWater,all=T)

Landuse_Herbaceous_byHUC_NoWater[is.na(Landuse_Herbaceous_byHUC_NoWater)] = 0

# Agriculture

LU_Agriculture_byHUC_NoWater <- LDR_landcover %>%
  filter(NLCD_Land_Cover_Class == "Hay/Pasture" | NLCD_Land_Cover_Class == "Cultivated Crops") %>%
  group_by(HUC12) %>%
  summarise(
    Agriculture_byHUC_NoWater = n()
  )

Landuse_Agriculture_byHUC_NoWater <- merge(LU_Agriculture_byHUC_NoWater,LU_tot_byHUC_NoWater,all=T)

Landuse_Agriculture_byHUC_NoWater[is.na(Landuse_Agriculture_byHUC_NoWater)] = 0

# Percent coverage

LU_Developed_Percent_byHUC_NoWater <- Landuse_Developed_byHUC_NoWater %>%
  mutate(LU_Percent_Developed_byHUC_NoWater = Developed_byHUC_NoWater / LU_tot_byHUC_NoWater)

LU_Barren_Percent_byHUC_NoWater <- Landuse_Barren_byHUC_NoWater %>%
  mutate(LU_Percent_Barren_byHUC_NoWater = Barren_byHUC_NoWater / LU_tot_byHUC_NoWater)

LU_Forest_Percent_byHUC_NoWater <- Landuse_Forest_byHUC_NoWater %>%
  mutate(LU_Percent_Forest_byHUC_NoWater = Forest_byHUC_NoWater / LU_tot_byHUC_NoWater)

LU_Shrubland_Percent_byHUC_NoWater <- Landuse_Shrubland_byHUC_NoWater %>%
  mutate(LU_Percent_Shrubland_byHUC_NoWater = Shrubland_byHUC_NoWater / LU_tot_byHUC_NoWater)

LU_Herbaceous_Percent_byHUC_NoWater <- Landuse_Herbaceous_byHUC_NoWater %>%
  mutate(LU_Percent_Herbaceous_byHUC_NoWater = Herbaceous_byHUC_NoWater / LU_tot_byHUC_NoWater)

LU_Agriculture_Percent_byHUC_NoWater <- Landuse_Agriculture_byHUC_NoWater %>%
  mutate(LU_Percent_Agriculture_byHUC_NoWater = Agriculture_byHUC_NoWater / LU_tot_byHUC_NoWater)

Results_list_byHUC_NoWater <- list(LU_Developed_Percent_byHUC_NoWater,LU_Barren_Percent_byHUC_NoWater,LU_Forest_Percent_byHUC_NoWater,
                           LU_Shrubland_Percent_byHUC_NoWater,LU_Herbaceous_Percent_byHUC_NoWater,LU_Agriculture_Percent_byHUC_NoWater)

Results_byHUC_bind_NoWater <- Reduce(function(dtf1, dtf2) merge(dtf1, dtf2, all = T),
                             Results_list_byHUC_NoWater) %>%
  select(HUC12,LU_Percent_Developed_byHUC_NoWater,LU_Percent_Barren_byHUC_NoWater,LU_Percent_Forest_byHUC_NoWater,
         LU_Percent_Shrubland_byHUC_NoWater,LU_Percent_Herbaceous_byHUC_NoWater,LU_Percent_Agriculture_byHUC_NoWater)

Results_byHUC_NoWater_final <- merge(Results_byHUC_bind_NoWater,LDR_WSArea_byHUC12,all=T)

# Part 4: Housing Market Analysis -------------------------------------------------------------

# All HUCs

# Subset to watersheds

Census_2000_df_JC <- st_filter(Census_2000_df_centroid, y = JohnsonCreek_proj, .predicate = st_intersects)

temp_Census_2000_df_JC <- Census_2000_df_JC %>%
  dplyr::select(GEOID) %>%
  distinct(GEOID)

Census_2000_df_JC_Area <- Census_2000_df %>%
  mutate(Area_Sqm = st_area(Census_2000_df)) %>%
  merge(temp_Census_2000_df_JC, all.y = TRUE) %>%
  st_drop_geometry()

Census_2000_df_BB <- st_filter(Census_2000_df_centroid, y = BurntBridge_proj, .predicate = st_intersects)

temp_Census_2000_df_BB <- Census_2000_df_BB %>%
  dplyr::select(GEOID) %>%
  distinct(GEOID)

Census_2000_df_BB_Area <- Census_2000_df %>%
  mutate(Area_Sqm = st_area(Census_2000_df)) %>%
  merge(temp_Census_2000_df_BB, all.y = TRUE) %>%
  st_drop_geometry()

Census_2000_df_LDR <- st_filter(Census_2000_df_centroid, y = LowerDelaware_proj, .predicate = st_intersects)

temp_Census_2000_df_LDR <- Census_2000_df_LDR %>%
  dplyr::select(GEOID) %>%
  distinct(GEOID)

Census_2000_df_LDR_Area <- Census_2000_df %>%
  mutate(Area_Sqm = st_area(Census_2000_df)) %>%
  merge(temp_Census_2000_df_LDR, all.y = TRUE) %>%
  st_drop_geometry()

All_States_2010_centroid <- All_States_2010 %>%
  st_centroid(All_States_2010) %>%
  st_zm()

All_States_2021_centroid <- All_States_2021 %>%
  st_centroid(All_States_2021) %>%
  st_zm()

ACS_2010_df_JC <- st_filter(All_States_2010_centroid, y = JohnsonCreek_proj, .predicate = st_intersects) %>%
  select(GEOID10,geometry)

temp_ACS_2010_df_JC <- ACS_2010_df_JC %>%
  dplyr::select(GEOID10) %>%
  distinct(GEOID10)

ACS_2010_df_JC_Area <- All_States_2010 %>%
  mutate(Area_Sqm = st_area(All_States_2010)) %>%
  merge(temp_ACS_2010_df_JC, all.y = TRUE) %>%
  st_drop_geometry()

ACS_2021_df_JC <- st_filter(All_States_2021_centroid, y = JohnsonCreek_proj, .predicate = st_intersects) %>%
  select(GEOID,geometry)

temp_ACS_2021_df_JC <- ACS_2021_df_JC %>%
  dplyr::select(GEOID) %>%
  distinct(GEOID)

ACS_2021_df_JC_Area <- All_States_2021 %>%
  mutate(Area_Sqm = st_area(All_States_2021)) %>%
  merge(temp_ACS_2021_df_JC, all.y = TRUE) %>%
  st_drop_geometry()

ACS_2010_df_BB <- st_filter(All_States_2010_centroid, y = BurntBridge_proj, .predicate = st_intersects) %>%
  select(GEOID10,geometry)

temp_ACS_2010_df_BB <- ACS_2010_df_BB %>%
  dplyr::select(GEOID10) %>%
  distinct(GEOID10)

ACS_2010_df_BB_Area <- All_States_2010 %>%
  mutate(Area_Sqm = st_area(All_States_2010)) %>%
  merge(temp_ACS_2010_df_BB, all.y = TRUE) %>%
  st_drop_geometry()

ACS_2021_df_BB <- st_filter(All_States_2021_centroid, y = BurntBridge_proj, .predicate = st_intersects) %>%
  select(GEOID,geometry)

temp_ACS_2021_df_BB <- ACS_2021_df_BB %>%
  dplyr::select(GEOID) %>%
  distinct(GEOID)

ACS_2021_df_BB_Area <- All_States_2021 %>%
  mutate(Area_Sqm = st_area(All_States_2021)) %>%
  merge(temp_ACS_2021_df_BB, all.y = TRUE) %>%
  st_drop_geometry()

ACS_2010_df_LDR <- st_filter(All_States_2010_centroid, y = LowerDelaware_proj, .predicate = st_intersects) %>%
  select(GEOID10,geometry)

temp_ACS_2010_df_LDR <- ACS_2010_df_LDR %>%
  dplyr::select(GEOID10) %>%
  distinct(GEOID10)

ACS_2010_df_LDR_Area <- All_States_2010 %>%
  mutate(Area_Sqm = st_area(All_States_2010)) %>%
  merge(temp_ACS_2010_df_LDR, all.y = TRUE) %>%
  st_drop_geometry()

ACS_2021_df_LDR <- st_filter(All_States_2021_centroid, y = LowerDelaware_proj, .predicate = st_intersects) %>%
  select(GEOID,geometry)

temp_ACS_2021_df_LDR <- ACS_2021_df_LDR %>%
  dplyr::select(GEOID) %>%
  distinct(GEOID)

ACS_2021_df_LDR_Area <- All_States_2021 %>%
  mutate(Area_Sqm = st_area(All_States_2021)) %>%
  merge(temp_ACS_2021_df_LDR, all.y = TRUE) %>%
  st_drop_geometry()

# Merge with demographics

Demo_merge_2010 <- merge(House_Inc_2010, Population_2010, all = T)
Demographics_2010 <- merge(Demo_merge_2010,House_Characteristics_2010, all = T) %>%
  rename("GEOID10" = "GEO_ID",
         "Population" = "B01001_001E",
         "Med_Value" = "DP04_0088E",
         "Med_Inc" = "B19013_001E",
         "Households" = "DP04_0001E",
         "Occupied_Housing" = "DP04_0002E",
         "One_Detached" = "DP04_0007E",
         "One_Attached" = "DP04_0008E",
         "Two_units" = "DP04_0009E",
         "Three_or_Four" = "DP04_0010E",
         "Five_to_Nine" = "DP04_0011E",
         "Ten_to_Nineteen" = "DP04_0012E",
         "Twenty_or_More" = "DP04_0013E",
         "Mobile_home" = "DP04_0014E",
         "Boat_RV" = "DP04_0015E")

Demographics_2010$GEOID10 <- gsub("1400000US","",Demographics_2010$GEOID10)

Demo_merge_2021 <- merge(House_Inc_2021, Population_2021, all = T)
Demographics_2021 <- merge(Demo_merge_2021,House_Characteristics_2021, all = T) %>%
  rename("GEOID" = "GEO_ID",
         "Population" = "B01001_001E",
         "Med_Value" = "DP04_0089E",
         "Med_Inc" = "B19013_001E",
         "Households" = "DP04_0001E",
         "Occupied_Housing" = "DP04_0002E",
         "One_Detached" = "DP04_0007E",
         "One_Attached" = "DP04_0008E",
         "Two_units" = "DP04_0009E",
         "Three_or_Four" = "DP04_0010E",
         "Five_to_Nine" = "DP04_0011E",
         "Ten_to_Nineteen" = "DP04_0012E",
         "Twenty_or_More" = "DP04_0013E",
         "Mobile_home" = "DP04_0014E",
         "Boat_RV" = "DP04_0015E")

Demographics_2021$GEOID <- gsub("1400000US","",Demographics_2021$GEOID)

JC_2010 <- merge(ACS_2010_df_JC_Area,Demographics_2010,all.x=T)
JC_2021 <- merge(ACS_2021_df_JC_Area,Demographics_2021,all.x=T)
BB_2010 <- merge(ACS_2010_df_BB_Area,Demographics_2010,all.x=T)
BB_2021 <- merge(ACS_2021_df_BB_Area,Demographics_2021,all.x=T)
LDR_2010 <- merge(ACS_2010_df_LDR_Area,Demographics_2010,all.x=T)
LDR_2021 <- merge(ACS_2021_df_LDR_Area,Demographics_2021,all.x=T)

column_list <- c("Med_Inc", "Population", "Households","Med_Value","Occupied_Housing","One_Detached","One_Attached","Two_units",
                 "Three_or_Four","Five_to_Nine","Ten_to_Nineteen","Twenty_or_More","Mobile_home","Boat_RV")

ACS_2010_list <- list(JC_2010,BB_2010,LDR_2010)
ACS_2021_list <- list(JC_2021,BB_2021,LDR_2021)
names(ACS_2010_list) <- c("JohnsonCreek","BurntBridge","LowerDelawareRiver")
names(ACS_2021_list) <- c("JohnsonCreek","BurntBridge","LowerDelawareRiver")

Census_2000_list <- list(Census_2000_df_JC_Area,Census_2000_df_LDR_Area,Census_2000_df_BB_Area)
names(Census_2000_list) <- c("JohnsonCreek","LowerDelawareRiver","BurntBridge")

# Summary Statistics - 2010 ACS

for (i in seq_along(ACS_2010_list)) {

  ACS_2010 <- ACS_2010_list[[i]]

  ACS_2010[column_list] <- lapply(ACS_2010[column_list], gsub, pattern = "-", replacement = NA)

  ACS_2010_summary <- ACS_2010 %>%
    mutate(Area_Sqkm = Area_Sqm/1000000,
           Pop_Density = as.numeric(Population)/Area_Sqkm,
           Occupied_Percent = as.numeric(Occupied_Housing)/as.numeric(Households),
           One_Total = (as.numeric(One_Detached) + as.numeric(One_Attached))*Occupied_Percent,
           Two_Total = (as.numeric(Two_units)*Occupied_Percent),
           Three_Four_Total = (as.numeric(Three_or_Four)*Occupied_Percent),
           Five_or_More = (as.numeric(Five_to_Nine) + as.numeric(Ten_to_Nineteen,Twenty_or_More))*Occupied_Percent,
           Other = (as.numeric(Mobile_home) + as.numeric(Boat_RV))*Occupied_Percent) %>%
    summarise(
      Median_Income = median(as.numeric(Med_Inc), na.rm = T),
      Households_Sum = sum(as.numeric(Households), na.rm = T),
      Occupied_Sum = sum(as.numeric(Occupied_Housing),na.rm=T),
      Median_Value = median(as.numeric(Med_Value), na.rm = T),
      Population_Sum = sum(as.numeric(Population), na.rm = T),
      Avg_Pop_Density_Sqkm = mean(Pop_Density, na.rm = T),
      Area_Sqkm = sum(Area_Sqkm,na.rm = T),
      One_Unit = sum(One_Total,na.rm = T),
      Two_Units = sum(Two_Total,na.rm=T),
      Three_Four_Units = sum(Three_Four_Total,na.rm=T),
      Five_orMore_Units = sum(Five_or_More, na.rm=T),
      Other_Unit = sum(Other,na.rm=T)
    )

  if (exists("ACS_2010_summary_bind")) {
    ACS_2010_summary_bind <- rbind(ACS_2010_summary_bind, ACS_2010_summary)
  } else {
    ACS_2010_summary_bind <- ACS_2010_summary
  }
}

ACS_2010_summary_bind$Watershed <- names(ACS_2010_list)

ACS_2010_summary_bind$Survey <- c("ACS_2010")

ACS_2010_final <- ACS_2010_summary_bind %>%
  select(Watershed,Survey,Area_Sqkm,Population_Sum,Avg_Pop_Density_Sqkm,Households_Sum,Occupied_Sum,One_Unit,Two_Units,Three_Four_Units,Five_orMore_Units,Other_Unit,Median_Value,Median_Income)

# Summary Statistics - 2021 ACS

for (i in seq_along(ACS_2021_list)) {

  ACS_2021 <- ACS_2021_list[[i]]

  ACS_2021[column_list] <- lapply(ACS_2021[column_list], gsub, pattern = "-", replacement = NA)

  ACS_2021_summary <- ACS_2021 %>%
    mutate(Area_Sqkm = Area_Sqm/1000000,
           Pop_Density = as.numeric(Population)/Area_Sqkm,
           Occupied_Percent = as.numeric(Occupied_Housing)/as.numeric(Households),
           One_Total = (as.numeric(One_Detached) + as.numeric(One_Attached))*Occupied_Percent,
           Two_Total = (as.numeric(Two_units)*Occupied_Percent),
           Three_Four_Total = (as.numeric(Three_or_Four)*Occupied_Percent),
           Five_or_More = (as.numeric(Five_to_Nine) + as.numeric(Ten_to_Nineteen,Twenty_or_More))*Occupied_Percent,
           Other = (as.numeric(Mobile_home) + as.numeric(Boat_RV))*Occupied_Percent) %>%
    summarise(
      Median_Income = median(as.numeric(Med_Inc), na.rm = T),
      Households_Sum = sum(as.numeric(Households), na.rm = T),
      Occupied_Sum = sum(as.numeric(Occupied_Housing),na.rm=T),
      Median_Value = median(as.numeric(Med_Value), na.rm = T),
      Population_Sum = sum(as.numeric(Population), na.rm = T),
      Avg_Pop_Density_Sqkm = mean(Pop_Density, na.rm = T),
      Area_Sqkm = sum(Area_Sqkm,na.rm = T),
      One_Unit = sum(One_Total,na.rm = T),
      Two_Units = sum(Two_Total,na.rm=T),
      Three_Four_Units = sum(Three_Four_Total,na.rm=T),
      Five_orMore_Units = sum(Five_or_More, na.rm=T),
      Other_Unit = sum(Other,na.rm=T)
    )

  if (exists("ACS_2021_summary_bind")) {
    ACS_2021_summary_bind <- rbind(ACS_2021_summary_bind, ACS_2021_summary)
  } else {
    ACS_2021_summary_bind <- ACS_2021_summary
  }
}

ACS_2021_summary_bind$Watershed <- names(ACS_2021_list)

ACS_2021_summary_bind$Survey <- c("ACS_2021")

ACS_2021_final <- ACS_2021_summary_bind %>%
  select(Watershed,Survey,Area_Sqkm,Population_Sum,Avg_Pop_Density_Sqkm,Households_Sum,Occupied_Sum,One_Unit,Two_Units,Three_Four_Units,Five_orMore_Units,Other_Unit,Median_Value,Median_Income)

# Summary statistics - 2000 Census

for (i in seq_along(Census_2000_list)) {

  Census_2000 <- Census_2000_list[[i]]

  Census_2000_summary <- Census_2000 %>%
    mutate(Area_Sqkm = Area_Sqm/1000000,
           Pop_Density = as.numeric(Population)/Area_Sqkm,
           Occupied_Percent = as.numeric(Occupied_Housing)/as.numeric(Households),
           One_Total = (as.numeric(One_Detached) + as.numeric(One_Attached))*Occupied_Percent,
           Two_Total = (as.numeric(Two_units)*Occupied_Percent),
           Three_Four_Total = (as.numeric(Three_or_Four)*Occupied_Percent),
           Five_or_More = (as.numeric(Five_to_Nine) + as.numeric(Ten_to_Nineteen,Twenty_or_More))*Occupied_Percent,
           Other = (as.numeric(Mobile_home) + as.numeric(Boat_RV))*Occupied_Percent) %>%
    summarise(
      Median_Income = median(Med_Inc, na.rm = T),
      Households_Sum = sum(Households, na.rm = T),
      Occupied_Sum = sum(as.numeric(Occupied_Housing),na.rm=T),
      Median_Value = median(Med_Value, na.rm = T),
      Population_Sum = sum(Population, na.rm = T),
      Avg_Pop_Density_Sqkm = mean(Pop_Density, na.rm = T),
      Area_Sqkm = sum(Area_Sqkm,na.rm = T),
      One_Unit = sum(One_Total,na.rm = T),
      Two_Units = sum(Two_Total,na.rm=T),
      Three_Four_Units = sum(Three_Four_Total,na.rm=T),
      Five_orMore_Units = sum(Five_or_More, na.rm=T),
      Other_Unit = sum(Other,na.rm=T)
    )

  if (exists("Census_2000_summary_bind")) {
    Census_2000_summary_bind <- rbind(Census_2000_summary_bind, Census_2000_summary)
  } else {
    Census_2000_summary_bind <- Census_2000_summary
  }
}

Census_2000_summary_bind$Watershed <- names(Census_2000_list)

Census_2000_summary_bind$Survey <- c("Census_2000")

Census_2000_final <- Census_2000_summary_bind %>%
  select(Watershed,Survey,Area_Sqkm,Population_Sum,Avg_Pop_Density_Sqkm,Households_Sum,Occupied_Sum,One_Unit,Two_Units,Three_Four_Units,Five_orMore_Units,Other_Unit,Median_Value,Median_Income)

Final_housing <- list(Census_2000_final,ACS_2010_final,ACS_2021_final)

Results_Housing <- do.call(rbind,Final_housing)

Results_Final_Housing <- merge(Results_Housing,WS_area,all=T)

write.csv(Results_Final_Housing,paste0(OutPath,"/HousingComparison_AllHUCs-2023-06-26.csv"), row.names = F)

# Intersect individual HUC12s

LDR_HUC12_list <- list(LDR_HUC12_1_proj,LDR_HUC12_2_proj,LDR_HUC12_3_proj,LDR_HUC12_4_proj,LDR_HUC12_5_proj,LDR_HUC12_6_proj,
                       LDR_HUC12_7_proj,LDR_HUC12_8_proj,LDR_HUC12_9_proj,LDR_HUC12_10_proj,LDR_HUC12_11_proj,LDR_HUC12_12_proj,
                       LDR_HUC12_13_proj,LDR_HUC12_14_proj,LDR_HUC12_15_proj)

Census_2000_LDR_byHUC12 <- list()
ACS_2010_LDR_byHUC12 <- list()
ACS_2021_LDR_byHUC12 <- list()

for (i in seq_along(LDR_HUC12_list)) {
  LDR_HUC12_proj <- LDR_HUC12_list[[i]]

  Census_2000_LDR_byHUC12[[i]] <- st_filter(Census_2000_df_centroid, y = LDR_HUC12_proj, .predicate = st_intersects)

  ACS_2010_LDR_byHUC12[[i]] <- st_filter(All_States_2010_centroid, y = LDR_HUC12_proj, .predicate = st_intersects)

  ACS_2021_LDR_byHUC12[[i]] <- st_filter(All_States_2021_centroid, y = LDR_HUC12_proj, .predicate = st_intersects)

  Census_2000_LDR_byHUC12[[i]]$HUC12 <- LDR_HUC12_proj$HUC12

  ACS_2010_LDR_byHUC12[[i]]$HUC12 <- LDR_HUC12_proj$HUC12

  ACS_2021_LDR_byHUC12[[i]]$HUC12 <- LDR_HUC12_proj$HUC12

  if (exists("Census_2000_LDR_byHUC12_bind")) {
    Census_2000_LDR_byHUC12_bind <- rbind(Census_2000_LDR_byHUC12_bind, Census_2000_LDR_byHUC12[[i]])
  } else {
    Census_2000_LDR_byHUC12_bind <- Census_2000_LDR_byHUC12[[i]]
  }

  if (exists("ACS_2010_LDR_byHUC12_bind")) {
    ACS_2010_LDR_byHUC12_bind <- rbind(ACS_2010_LDR_byHUC12_bind, ACS_2010_LDR_byHUC12[[i]])
  } else {
    ACS_2010_LDR_byHUC12_bind <- ACS_2010_LDR_byHUC12[[i]]
  }

  if (exists("ACS_2021_LDR_byHUC12_bind")) {
    ACS_2021_LDR_byHUC12_bind <- rbind(ACS_2021_LDR_byHUC12_bind, ACS_2021_LDR_byHUC12[[i]])
  } else {
    ACS_2021_LDR_byHUC12_bind <- ACS_2021_LDR_byHUC12[[i]]
  }
}

# Add area

Census_2000_LDR_byHUC12_bind_Area <- merge(Census_2000_LDR_byHUC12_bind,LDR_WSArea_byHUC12) %>%
  st_drop_geometry()

Census_2000_LDR_byHUC12_Area <- merge(Census_2000_LDR_byHUC12_bind_Area,Census_2000_df_LDR_Area)

ACS_2010_LDR_byHUC12_bind_Area <- merge(ACS_2010_LDR_byHUC12_bind,LDR_WSArea_byHUC12) %>%
  st_drop_geometry()

ACS_2010_LDR_byHUC12_Area <- merge(ACS_2010_LDR_byHUC12_bind_Area,ACS_2010_df_LDR_Area)

ACS_2021_LDR_byHUC12_bind_Area <- merge(ACS_2021_LDR_byHUC12_bind,LDR_WSArea_byHUC12) %>%
  st_drop_geometry()

ACS_2021_LDR_byHUC12_Area <- merge(ACS_2021_LDR_byHUC12_bind_Area,ACS_2021_df_LDR_Area)

ACS_LDR_2010 <- merge(ACS_2010_LDR_byHUC12_Area,Demographics_2010,all.x=T)

ACS_LDR_2010[column_list] <- lapply(ACS_LDR_2010[column_list], gsub, pattern = "-", replacement = NA)

ACS_LDR_2021 <- merge(ACS_2021_LDR_byHUC12_Area,Demographics_2021,all.x=T)

ACS_LDR_2021[column_list] <- lapply(ACS_LDR_2021[column_list], gsub, pattern = "-", replacement = NA)

# Summarize

Census_2000_summary_byHUC12 <- Census_2000_LDR_byHUC12_Area %>%
  mutate(Area_Sqkm = Area_Sqm/1000000,
         Pop_Density = as.numeric(Population)/Area_Sqkm,
         Occupied_Percent = as.numeric(Occupied_Housing)/as.numeric(Households),
         One_Total = (as.numeric(One_Detached) + as.numeric(One_Attached))*Occupied_Percent,
         Two_Total = (as.numeric(Two_units)*Occupied_Percent),
         Three_Four_Total = (as.numeric(Three_or_Four)*Occupied_Percent),
         Five_or_More = (as.numeric(Five_to_Nine) + as.numeric(Ten_to_Nineteen,Twenty_or_More))*Occupied_Percent,
         Other = (as.numeric(Mobile_home) + as.numeric(Boat_RV))*Occupied_Percent) %>%
  group_by(HUC12) %>%
  summarise(
    Median_Income = median(as.numeric(Med_Inc), na.rm = T),
    Households_Sum = sum(as.numeric(Households), na.rm = T),
    Occupied_Sum = sum(as.numeric(Occupied_Housing),na.rm=T),
    Median_Value = median(as.numeric(Med_Value), na.rm = T),
    Population_Sum = sum(as.numeric(Population), na.rm = T),
    Avg_Pop_Density_Sqkm = mean(Pop_Density, na.rm = T),
    Area_Sqkm = sum(Area_Sqkm,na.rm = T),
    One_Unit = sum(One_Total,na.rm = T),
    Two_Units = sum(Two_Total,na.rm=T),
    Three_Four_Units = sum(Three_Four_Total,na.rm=T),
    Five_orMore_Units = sum(Five_or_More, na.rm=T),
    Other_Unit = sum(Other,na.rm=T)
  )

ACS_2010_summary_byHUC12 <- ACS_LDR_2010 %>%
  mutate(Area_Sqkm = Area_Sqm/1000000,
         Pop_Density = as.numeric(Population)/Area_Sqkm,
         Occupied_Percent = as.numeric(Occupied_Housing)/as.numeric(Households),
         One_Total = (as.numeric(One_Detached) + as.numeric(One_Attached))*Occupied_Percent,
         Two_Total = (as.numeric(Two_units)*Occupied_Percent),
         Three_Four_Total = (as.numeric(Three_or_Four)*Occupied_Percent),
         Five_or_More = (as.numeric(Five_to_Nine) + as.numeric(Ten_to_Nineteen,Twenty_or_More))*Occupied_Percent,
         Other = (as.numeric(Mobile_home) + as.numeric(Boat_RV))*Occupied_Percent) %>%
  group_by(HUC12) %>%
  summarise(
    Median_Income = median(as.numeric(Med_Inc), na.rm = T),
    Households_Sum = sum(as.numeric(Households), na.rm = T),
    Occupied_Sum = sum(as.numeric(Occupied_Housing),na.rm=T),
    Median_Value = median(as.numeric(Med_Value), na.rm = T),
    Population_Sum = sum(as.numeric(Population), na.rm = T),
    Avg_Pop_Density_Sqkm = mean(Pop_Density, na.rm = T),
    Area_Sqkm = sum(Area_Sqkm,na.rm = T),
    One_Unit = sum(One_Total,na.rm = T),
    Two_Units = sum(Two_Total,na.rm=T),
    Three_Four_Units = sum(Three_Four_Total,na.rm=T),
    Five_orMore_Units = sum(Five_or_More, na.rm=T),
    Other_Unit = sum(Other,na.rm=T)
  )

ACS_2021_summary_byHUC12 <- ACS_LDR_2021 %>%
  mutate(Area_Sqkm = Area_Sqm/1000000,
         Pop_Density = as.numeric(Population)/Area_Sqkm,
         Occupied_Percent = as.numeric(Occupied_Housing)/as.numeric(Households),
         One_Total = (as.numeric(One_Detached) + as.numeric(One_Attached))*Occupied_Percent,
         Two_Total = (as.numeric(Two_units)*Occupied_Percent),
         Three_Four_Total = (as.numeric(Three_or_Four)*Occupied_Percent),
         Five_or_More = (as.numeric(Five_to_Nine) + as.numeric(Ten_to_Nineteen,Twenty_or_More))*Occupied_Percent,
         Other = (as.numeric(Mobile_home) + as.numeric(Boat_RV))*Occupied_Percent) %>%
  group_by(HUC12) %>%
  summarise(
    Median_Income = median(as.numeric(Med_Inc), na.rm = T),
    Households_Sum = sum(as.numeric(Households), na.rm = T),
    Occupied_Sum = sum(as.numeric(Occupied_Housing),na.rm=T),
    Median_Value = median(as.numeric(Med_Value), na.rm = T),
    Population_Sum = sum(as.numeric(Population), na.rm = T),
    Avg_Pop_Density_Sqkm = mean(Pop_Density, na.rm = T),
    Area_Sqkm = sum(Area_Sqkm,na.rm = T),
    One_Unit = sum(One_Total,na.rm = T),
    Two_Units = sum(Two_Total,na.rm=T),
    Three_Four_Units = sum(Three_Four_Total,na.rm=T),
    Five_orMore_Units = sum(Five_or_More, na.rm=T),
    Other_Unit = sum(Other,na.rm=T)
  )

Census_2000_summary_byHUC12$Watershed <- c("LowerDelawareRiver")

Census_2000_summary_byHUC12$Survey <- c("Census_2000")

Census_2000_final_byHUC12 <- Census_2000_summary_byHUC12 %>%
  select(Watershed,HUC12,Survey,Area_Sqkm,Population_Sum,Avg_Pop_Density_Sqkm,Households_Sum,Occupied_Sum,One_Unit,Two_Units,Three_Four_Units,Five_orMore_Units,Other_Unit,Median_Value,Median_Income)

ACS_2010_summary_byHUC12$Watershed <- c("LowerDelawareRiver")

ACS_2010_summary_byHUC12$Survey <- c("ACS_2010")

ACS_2010_final_byHUC12 <- ACS_2010_summary_byHUC12 %>%
  select(Watershed,HUC12,Survey,Area_Sqkm,Population_Sum,Avg_Pop_Density_Sqkm,Households_Sum,Occupied_Sum,One_Unit,Two_Units,Three_Four_Units,Five_orMore_Units,Other_Unit,Median_Value,Median_Income)

ACS_2021_summary_byHUC12$Watershed <- c("LowerDelawareRiver")

ACS_2021_summary_byHUC12$Survey <- c("ACS_2021")

ACS_2021_final_byHUC12 <- ACS_2021_summary_byHUC12 %>%
  select(Watershed,HUC12,Survey,Area_Sqkm,Population_Sum,Avg_Pop_Density_Sqkm,Households_Sum,Occupied_Sum,One_Unit,Two_Units,Three_Four_Units,Five_orMore_Units,Other_Unit,Median_Value,Median_Income)

Final_housing_byHUC12 <- list(Census_2000_final_byHUC12,ACS_2010_final_byHUC12,ACS_2021_final_byHUC12)

Results_Housing_byHUC12 <- do.call(rbind,Final_housing_byHUC12)

Results_Housing_byHUC12_final <- merge(Results_Housing_byHUC12,LDR_WSArea_byHUC12)

write.csv(Results_Housing_byHUC12_final,paste0(OutPath,"/HousingComparison_byHUC12-2023-06-26.csv"), row.names = F)

# QA

Results_QA <- Results_final %>%
  group_by(Watershed) %>%
  summarise(Sum_Percent = sum(LU_Percent_Agriculture,LU_Percent_Water,LU_Percent_Barren,LU_Percent_Developed,LU_Percent_Forest,LU_Percent_Herbaceous,LU_Percent_Shrubland,LU_Percent_Wetlands))

Results_QA_NoWater <- Results_final_NoWater %>%
  group_by(Watershed) %>%
  summarise(Sum_Percent = sum(LU_Percent_Agriculture_NoWater,LU_Percent_Barren_NoWater,LU_Percent_Developed_NoWater,LU_Percent_Forest_NoWater,LU_Percent_Herbaceous_NoWater,LU_Percent_Shrubland_NoWater))

Results_QA_byHUC <- Results_byHUC_final %>%
  group_by(HUC12) %>%
  summarise(Sum_Percent = sum(LU_Percent_Agriculture_byHUC,LU_Percent_Water_byHUC,LU_Percent_Barren_byHUC,LU_Percent_Developed_byHUC,LU_Percent_Forest_byHUC,LU_Percent_Herbaceous_byHUC,LU_Percent_Shrubland_byHUC,LU_Percent_Wetlands_byHUC))

Results_QA_byHUC_NoWater <- Results_byHUC_NoWater_final %>%
  group_by(HUC12) %>%
  summarise(Sum_Percent = sum(LU_Percent_Agriculture_byHUC_NoWater,LU_Percent_Barren_byHUC_NoWater,LU_Percent_Developed_byHUC_NoWater,LU_Percent_Forest_byHUC_NoWater,LU_Percent_Herbaceous_byHUC_NoWater,LU_Percent_Shrubland_byHUC_NoWater))