Lab Project #2

Section #1: Import necessary libraries

# Install packages
if (!require("tidycensus")) {install.packages("tidycensus")}

Loading required package: tidycensus

if (!require("tidyverse")) {install.packages("tidyverse")}

Loading required package: tidyverse

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

if (!require("readxl")) {install.packages("readxl")}

Loading required package: readxl

if (!require("elevatr")) {install.packages("elevatr")}

Loading required package: elevatr
elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'.  Use 
of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being 
deprecated; however, get_elev_raster continues to return a RasterLayer.  This 
will be dropped in future versions, so please plan accordingly.

if (!require("sf")) {install.packages("sf")}

Loading required package: sf
Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE

if (!require("ggplot2")) {install.packages("ggplot2")}
if (!require("raster")) {install.packages("raster")}

Loading required package: raster
Loading required package: sp

Attaching package: 'raster'

The following object is masked from 'package:dplyr':

    select

if (!require("terra")) {install.packages("terra")}

Loading required package: terra
terra 1.7.83

Attaching package: 'terra'

The following object is masked from 'package:tidyr':

    extract

if (!require("ggspatial")) {install.packages("ggspatial")}

Loading required package: ggspatial

if (!require("ggnewscale")) {install.packages("ggnewscale")}

Loading required package: ggnewscale

if (!require("tidyterra")) {install.packages("tidyterra")}

Loading required package: tidyterra

Attaching package: 'tidyterra'

The following object is masked from 'package:raster':

    select

The following object is masked from 'package:stats':

    filter

if (!require("topoDistance")) {install.packages("topoDistance")}

Loading required package: topoDistance

if (!require("progress")) {install.packages("progress")}

Loading required package: progress

if (!require("gdistance")) {install.packages("gdistance")}

Loading required package: gdistance
Loading required package: igraph

Attaching package: 'igraph'

The following object is masked from 'package:tidyterra':

    groups

The following objects are masked from 'package:terra':

    blocks, compare, union

The following object is masked from 'package:raster':

    union

The following objects are masked from 'package:lubridate':

    %--%, union

The following objects are masked from 'package:dplyr':

    as_data_frame, groups, union

The following objects are masked from 'package:purrr':

    compose, simplify

The following object is masked from 'package:tidyr':

    crossing

The following object is masked from 'package:tibble':

    as_data_frame

The following objects are masked from 'package:stats':

    decompose, spectrum

The following object is masked from 'package:base':

    union

Loading required package: Matrix

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack


Attaching package: 'gdistance'

The following object is masked from 'package:igraph':

    normalize

if (!require("sp")) {install.packages("sp")}
if (!require("leastcostpath")) {install.packages("leastcostpath")}

Loading required package: leastcostpath

if (!require("rgdal")) {install.packages("rgdal", dependencies = TRUE)}

Loading required package: rgdal

Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
logical.return = TRUE, : there is no package called 'rgdal'

Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.4'
(as 'lib' is unspecified)

Warning: package 'rgdal' is not available for this version of R

A version of this package for your version of R might be available elsewhere,
see the ideas at
https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages

if (!require("ggmap")) {install.packages("ggmap", dependencies = TRUE)}

Loading required package: ggmap
ℹ Google's Terms of Service: <https://mapsplatform.google.com>
  Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service/>
  OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles/>
ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
Attaching package: 'ggmap'

The following object is masked from 'package:terra':

    inset

if (!require("tmaptools")) {install.packages("tmaptools", dependencies = TRUE)}

Loading required package: tmaptools

# Load the libraries
library(tidycensus)
library(tidyverse)
library(readxl)
library(elevatr) 
library(sf) 
library(ggplot2)
library(raster)
library(terra)
library(ggspatial)
library(ggnewscale)
library(tidyterra)
library(topoDistance)
library(progress)
library(gdistance)
library(sp)
library(leastcostpath)
library(ggmap)
library(tmaptools)

Section #2: Import Voter Registration Data

# file path to data frame.
file_path <- "save/voter_registration.rds"

# if the data frame exists, load it.
if (file.exists(file_path)) {
  voter_registration <- readRDS(file_path)
  message("Dataframe loaded from file.")
} else {
    # Assuming you have a Excel file with voter registration data
  voter_registration <- read_excel("/cloud/project/current voterregstatsbycongressionaldistricts.xlsx")
  
  i <- c(1, 3:9)
  voter_registration[ , i] <- apply(voter_registration[ , i], 2, function(x) as.numeric(as.character(x)))
  
  # Calculate proportion of Republican voters
  voter_registration <- voter_registration %>%
    mutate(RepublicanProportion = Republican / Total)
  
  # Calculate proportion of Democratic voters
  voter_registration <- voter_registration %>%
    mutate(DemocraticProportion = Democratic / Total)
  
  # Calculate proportion of Libertarian voters
  voter_registration <- voter_registration %>%
    mutate(LibertarianProportion = Libertarian / Total)
  
  # Calculate proportion of Green voters
  voter_registration <- voter_registration %>%
    mutate(GreenProportion = Green / Total)
  
  # Calculate proportion of No Affiliation voters
  voter_registration <- voter_registration %>%
    mutate(NoAffiliationProportion = `No Affiliation` / Total)
  
  # Calculate proportion of Other voters
  voter_registration <- voter_registration %>%
    mutate(OtherProportion = Other / Total)
  
  saveRDS(voter_registration, file_path) 
  message("Dataframe generated and saved to file.")
}

Dataframe loaded from file.

Section #3: Import ACS variable table for later use.

# Define the file path
file_path <- "save/acs_vars.rda"

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data frame from the file
  acs_variables <- readRDS(file_path)
  message("Dataframe loaded from file.")
  
} else {
  
  # Generate the data frame
  acs_variables <- load_variables(dataset = "acs1", year = 2021)
  
  # Save the data frame to the file
  saveRDS(acs_variables, file_path)
  message("Dataframe generated and saved to file.")
}

Dataframe loaded from file.

Section #4: Import PA state geometry.

# Define the file path
file_path <- "save/pa_geom.rda"

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data from the file
  pa_shape <- readRDS(file_path)
  message("Data loaded from file.")
  
} else {
  # Your code to generate the data
  pa_shape <- get_acs(geography = "state",
                      variables = c("B25077_001"),  # Random Variable to ignore.
                      state = "PA",
                      geometry = TRUE) %>%
    subset(select = -c(variable, estimate, moe))
  
  # Save the data  to the file
  saveRDS(pa_shape, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Section #5: Import Population Density Data for PA

# Define the file path
file_path <- "save/pop_raster.rds"

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data from the file
  pop_raster <- readRDS(file_path)
  message("Data loaded from file.")
  
} else {
  
  # Generate the data
  pop_raster <- raster("/cloud/project/pden2010_block/pden2010_60m.tif")
  
  # Transform the CRS of the multi-polygon to match the raster CRS
  pa_multipolygon <- st_transform(pa_shape, crs(pop_raster))
  
  # Crop the raster by the multi-polygon
  cropped_raster <- crop(pop_raster, extent(pa_multipolygon))
  plot(log(cropped_raster+1))
  
  # Mask the raster to set values outside the multi-polygon to NA
  masked_raster <- mask(cropped_raster, pa_multipolygon)
  plot(log(masked_raster+1))
  
  # Log transform of the population.
  pop_raster <- (log(masked_raster+1))
    
  # Save the data to the file
  saveRDS(pop_raster, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Section #6: Import Housing Cost Data for PA

# Define the file path
file_path <- "save/pa_housing_data.rds"

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data from the file
  pa_dists <- readRDS(file_path)
  message("Data loaded from file.")
  
} else {
  
  # Query housing price data for Pennsylvania's congressional districts
  pa_housing_data <- get_acs(geography = "congressional district",
                             variables = c("B25077_001"),  # Median home value, Median monthly housing costs
                             state = "PA",
                             geometry = TRUE)
  
  # Format the data layout.
  colnames(voter_registration)[1] <- "GEOID"
  pa_housing_data$GEOID <- as.numeric(pa_housing_data$GEOID)
  
  # Merge the ACS data with voter registration data
  merged_data <- pa_housing_data %>%
    left_join(voter_registration, by = "GEOID")
  
  # Transform the district data
  pa_dists <- st_transform(merged_data, crs(pop_raster))
  
  # Save the data frame to the file
  saveRDS(pa_dists, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Section #7: Visualize the proportion of registered republican voters for PA congressional districts

# Create the ggplot
ggplot(data = pa_dists) +
  geom_sf(aes(fill = RepublicanProportion), color = "darkred", size = 0.5) +
  scale_fill_gradient2(low = "white", mid = "grey", high = "red", midpoint = 0.5) +
  coord_sf() +  # Coordinate system for spatial data
  theme_minimal() +
  labs(title = "Proportion of Republican Voters in Pennsylvania Congressional Districts",
       fill = "Republican Proportion")

Section #8: Visualize the proportion of registered democrat voters for PA congressional districts

# Create the ggplot
ggplot(data = pa_dists) +
  geom_sf(aes(fill = DemocraticProportion), color = "darkblue", size = 0.5) +
  scale_fill_gradient2(low = "white", mid = "grey", high = "blue", midpoint = 0.5) +
  coord_sf() +  # Coordinate system for spatial data
  theme_minimal() +
  labs(title = "Proportion of Democratic Voters in Pennsylvania Congressional Districts",
       fill = "Democratic Proportion")

Section #9: Visualize the estimate median housing costs for PA congressional districts

# Create the ggplot
ggplot(data = pa_dists) +
  geom_sf(aes(fill = estimate), color = "darkblue", size = 0.5) +
  scale_fill_gradient2(low = "white", mid = "grey", high = "darkgreen", midpoint = 0.5) +
  coord_sf() +  # Coordinate system for spatial data
  theme_minimal() +
  labs(title = "Median Home Value in Pennsylvania Congressional Districts",
       fill = "Median Home Value")

Section #10: Visualize the population density for PA (Log Scale)

ggplot() +
  geom_spatraster(data=rast(pop_raster), show.legend = TRUE) +
  geom_sf(data = pa_dists, alpha=0, color='black', size = 0.5) +
  coord_sf() +  # Coordinate system for spatial data
  theme_minimal() +
  labs(title = '')

<SpatRaster> resampled to 500720 cells.

Section #11: Import the DEM for PA

# Define the file path
file_path <- "save/pa_dem.rds"

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data from the file
  pa_dem <- readRDS(file_path)
  message("Data loaded from file.")
  
} else {
  
  # Get the bounding box of the geometry 
  bounding_box <- st_bbox(pa_shape)
  pa_sf <- st_as_sfc(pa_shape)
  
  # Get the DEM data
  pa_dem <- get_elev_raster(locations = pa_dists, z = 7, src = "aws", prj = crs(pop_raster))
  
  # Create Spat Raster
  pa_dem_spatrast <- rast(pa_dem)
  
  # Transform the CRS of the multi-polygon to match the raster CRS
  pa_multipolygon <- st_transform(pa_shape, crs(pa_dem_spatrast))

  # Crop the raster by the multi-polygon
  cropped_dem <- crop(pa_dem_spatrast, extent(pa_multipolygon))
  plot(log(cropped_dem+1))

  # Mask the raster to set values outside the multi-polygon to NA
  pa_dem <- mask(cropped_dem, pa_multipolygon)
  plot(log(pa_dem+1))
  
  # Save the data to the file
  saveRDS(pa_dem, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Section #12: Generater the slope from PA DEM

# Define the file path
file_path <- "save/pa_slope.rds"

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data from the file
  pa_slope <- readRDS(file_path)
  message("Data loaded from file.")
  
} else {
  
  # Generate the data
  pa_slope <- terra::terrain(pa_dem, v = "slope",
                                unit = "radians", neighbors = 8)
  
  # Save the data to the file
  saveRDS(pa_slope, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Section #13: Generater the aspect from PA DEM

# Define the file path
file_path <- "save/pa_aspect.rds"

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data from the file
  pa_aspect <- readRDS(file_path)
  message("Data loaded from file.")
  
} else {
  
  # Generate the data
  pa_aspect <- terra::terrain(pa_dem, v = "aspect",
                                 unit = "radians", neighbors = 8)
  
  # Save the data to the file
  saveRDS(pa_aspect, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Section #14: Generater the shade from PA DEM

# Define the file path
file_path <- "save/pa_shade.rds"

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data from the file
  pa_shade <- readRDS(file_path)
  message("Data loaded from file.")
  
} else {
  
  # Generate the data
  pa_shade <- terra::shade(slope = pa_slope,
                              aspect = pa_aspect,
                              angle = 30, direction = 315,
                              normalize = TRUE)
  
  # Save the data to the file
  saveRDS(pa_shade, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Section #15: Visualize PA Elevation

ggplot() +
  geom_spatraster(data=pa_shade, show.legend = FALSE) +
  scale_fill_distiller(palette = "Greys") +
  new_scale_fill() + 
  geom_spatraster(data=pa_dem, alpha = 0.7) + 
  #alpha is the proportion opaque for the layer, so we're saying to make 
  #this 30% transparent
  scale_fill_hypso_tint_c() +
  labs(fill = "Elevation (masl)")

<SpatRaster> resampled to 500720 cells.
<SpatRaster> resampled to 500720 cells.

Section _#: Import the coordinates for Hospitals in PA

# Define the file path
file_path <- "save/hospital_coords.rds"

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data from the file
  hospital_coords <- readRDS(file_path)
  message("Data loaded from file.")
  
} else {
  
  # Read Hospital Data
  hospital_data <- read_csv("Hospitals.csv")
  coords <- select(hospital_data, LONGITUDE, LATITUDE)
  coords <- st_as_sf(x = coords, 
                     coords = c("LONGITUDE", "LATITUDE"),
                     crs="EPSG:4258")
  hospital_coords <- st_transform(coords, '+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs')
  
  # Save the data to the file
  saveRDS(hospital_coords, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Section _#: Generate Cost Distance for Hospitals in PA

# Define the file path
file_path <- "save/hospital_cost_dist.rds"

# Define a cost function based on slope 
cost_function <- function(x) { 
  slope <- abs(diff(x)) # Calculate the slope 
  return(slope + 1) # Add 1 to avoid zero cost 
}

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data from the file
  hospitals_cost_distance <- readRDS(file_path)
  message("Data loaded from file.")
  
} else {
  
  # Create the transition layer 
  trans <- transition(as(pa_dem, "Raster"), 
                      transitionFunction = cost_function, 
                      directions = 8)
  conductance <- geoCorrection(trans)
  
  # Calculate the cost distance raster 
  hospitals_cost_distance <- accCost(conductance, fromCoords=as(hospital_coords, "Spatial"))
  plot(hospitals_cost_distance)
  
  # Save the data to the file
  saveRDS(hospitals_cost_distance, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Section _#: Import the coordinates for Schools in PA

# replace "api_key" with your API key
register_google(key = 'AIzaSyABGVOEA6qrU5PQNOH33QrjJgHzR04ivB4')

# file path for the data frame
file_path <- "save/schls_ggmap.Rda"

# if the data frame exists, load it.
if (file.exists(file_path)) {
  
  # Load the data from the file
  schls_ggmap <- readRDS(file_path)
  message("Data loaded from file.")
  
} else { 
  
  # read/extract data for schools.
  schools <- read_excel("ExtractPublicSchools.xlsx")
  schools <- subset(schools, 
                    GradeCtgy == 'SECONDARY')
  schools <- select(schools, 
                    SchoolName, 
                    LEALocationCity, 
                    LEALocationState)
  schools$Concat <- with(data = schools, 
                         paste(SchoolName, 
                               ',', 
                               LEALocationCity,
                               ',', 
                               LEALocationState))
  schools <- schools[!apply(schools, 1, function(x) any(x=="")),] 
  
  # run the geocode function from ggmap package
  schls_ggmap <- geocode(location = schools$Concat,
                         output = "more", 
                         source = "google")
  schls_ggmap <- cbind(schools, 
                       schls_ggmap)
  
  # print the results
  schls_ggmap <- schls_ggmap[complete.cases(schls_ggmap), ]
  
  # Save the data to the file
  saveRDS(schls_ggmap, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Section _#: Clean and Mutate the School Coords

# Define the file path
file_path <- "save/school_coords.rds"

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data from the file
  school_coords <- readRDS(file_path)
  message("Data loaded from file.")
  
} else {
  
  # Read Hospital Data
  coords <- select(schls_ggmap, lon, lat)
  coords <- st_as_sf(x = coords, 
                     coords = c("lon", "lat"),
                     crs="EPSG:4258")
  school_coords <- st_transform(coords, '+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs')
  
  # Save the data to the file
  saveRDS(school_coords, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Section _#: Generate Cost Distance for School Availability

# Define the file path
file_path <- "save/school_cost_dist.rds"

# Define a cost function based on slope 
cost_function <- function(x) { 
  slope <- abs(diff(x)) # Calculate the slope 
  return(slope + 1) # Add 1 to avoid zero cost 
}

# Check if the file exists
if (file.exists(file_path)) {
  
  # Load the data from the file
  school_cost_dist <- readRDS(file_path)
  message("Data loaded from file.")
  
} else {
  
  # Create the transition layer 
  trans <- transition(as(pa_dem, "Raster"), 
                      transitionFunction = cost_function, 
                      directions = 8)
  conductance <- geoCorrection(trans)
  
  # Calculate the cost distance raster 
  school_cost_dist <- accCost(conductance, fromCoords=as(school_coords, "Spatial"))
  plot(school_cost_dist)
  
  # Save the data to the file
  saveRDS(school_cost_dist, file_path)
  message("Data generated and saved to file.")
}

Data loaded from file.

Visualize the Cost Distance Rasters

ggplot() +
  geom_spatraster(data=rast(school_cost_dist), show.legend = TRUE) +
  geom_sf(data = pa_dists, alpha=0, color='black', size = 0.5) +
  coord_sf() +  # Coordinate system for spatial data
  theme_minimal() +
  labs(title = '')

<SpatRaster> resampled to 500720 cells.

ggplot() +
  geom_spatraster(data=rast(hospitals_cost_distance), show.legend = TRUE) +
  geom_sf(data = pa_dists, alpha=0, color='black', size = 0.5) +
  coord_sf() +  # Coordinate system for spatial data
  theme_minimal() +
  labs(title = '')

<SpatRaster> resampled to 500720 cells.

Cont.