Here I add my libraries and read in the various data files. I performs some cleaning as well.
library(tidyverse)
library(janitor)
library(sf)
library(mapview)
library(ggplot2)
library(ggmap)
wd <- getwd()
# read in the shapes
zip_path <- file.path(wd, "R-Spatial_II_Lab", "ZIP_CODE_040114", "ZIP_CODE_040114.shp")
zip_codes_sf <- st_read(zip_path)
nyc_food_path <- file.path(wd, "R-Spatial_II_Lab", "nycFoodStore.shp")
nyc_food_store_sf <- st_read(nyc_food_path)
census_path <- file.path(wd, "R-Spatial_II_Lab", "2010 Census Tracts", "geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a.shp")
census_tracts_2010_sf <- st_read(census_path)
# read in the csvs and clean names
health_facility_path <- file.path(wd, "R-Spatial_II_Lab", "NYS_Health_Facility.csv")
health_facility_df <- read_csv(health_facility_path, show_col_types = FALSE, lazy = FALSE) %>%
clean_names()
tests_path <- file.path(wd, "R-Spatial_II_Lab", "tests-by-zcta_2020_04_12.csv")
tests_by_zcta_2020_I_df <- read_csv(tests_path, show_col_types = FALSE, lazy = FALSE) %>%
clean_names()
tests_path <- file.path(wd, "R-Spatial_II_Lab", "tests-by-zcta_2020_04_19.csv")
tests_by_zcta_2020_II_df <- read_csv(tests_path, show_col_types = FALSE, lazy = FALSE) %>%
clean_names()
tests_path <- file.path(wd, "R-Spatial_II_Lab", "tests-by-zcta_2021_04_23.csv")
tests_by_zcta_2021_df <- read_csv(tests_path, show_col_types = FALSE, lazy = FALSE) %>%
clean_names()
acs_path <- file.path(wd, "R-Spatial_II_Lab", "ACSDP5Y2018.DP05_data_with_overlays_2020-04-22T132935.csv")
acs_demographics_housing_df <- read_csv(acs_path, skip = 1, show_col_types = FALSE, lazy = FALSE) %>%
clean_names()
# EDA
# glimpse(tests_by_zcta_2020_I_df)
# glimpse(tests_by_zcta_2020_II_df)
# glimpse(tests_by_zcta_2021_df)
#
# glimpse(zip_codes_sf)
#
# plot(st_geometry(zip_codes_sf))# join covid-19 data to zip code area data
# I used the 2021 data as it offered more attributes
# need to convert to chr first to make this work
tests_by_zcta_2021_df <- tests_by_zcta_2021_df %>%
mutate(modified_zcta = as.character(modified_zcta))
zip_covid_merged_sf <- zip_codes_sf %>%
left_join(tests_by_zcta_2021_df, by = c("ZIPCODE" = "modified_zcta"))
# Prepare the food
# I used the grepl function to select only records with any of the establishment codes:
# J, A, C, K. I included vehicles (K) but that might not be necessary if people social distance
nyc_food_store_prepared_sf <- nyc_food_store_sf %>%
filter(grepl("J", Estbl_T) |
grepl("A", Estbl_T) |
grepl("C", Estbl_T) |
grepl("K", Estbl_T))
# also have to transform it to NAD83
nyc_food_store_prepared_sf <- st_transform(nyc_food_store_prepared_sf, 2263)
# testing
stores_per_zipcode <- nyc_food_store_prepared_sf %>%
group_by(Zip_Cod) %>%
summarize(store_count = n()) %>%
ungroup()
print(stores_per_zipcode)
# now to aggregate with the zip codes
# first I create a sf with just summed stores in each zip
zip_food_query_sf <- nyc_food_store_prepared_sf %>%
st_join(zip_codes_sf, join = st_intersects) %>%
filter(!is.na(ZIPCODE)) %>%
group_by(ZIPCODE) %>%
summarize(store_count = n(), .groups = 'drop')
zip_food_query_df <- as.data.frame(zip_food_query_sf)
zip_food_query_df[3] <- NULL
# join to the zip codes by ZIPCODE and create the finished sf
# also remove some cruft
zip_fstore_agg_sf <- zip_codes_sf %>%
left_join(zip_food_query_df, by = 'ZIPCODE')
# join covid-19 data to zip code area data
# I used the 2021 data as it offered more attributes
# need to convert to chr first to make this work, also changed zcta to ZIPCODE for the join
tests_2021_df <- tests_by_zcta_2021_df %>%
mutate(modified_zcta = as.character(modified_zcta)) %>%
rename(ZIPCODE = modified_zcta)
zip_fstore_covid_sf <- left_join(zip_fstore_agg_sf, tests_2021_df, by = "ZIPCODE")
# aggregate the health facilities to zip codes
# I handle this similarly to the stores from above but because we start with df its a little more complicated
# I reuse some code from lab 1, this pipe gets only facilities in proper counties and filters for some anomalies
health_facility_sf <- health_facility_df %>%
filter(!is.na(facility_longitude) & !is.na(facility_latitude)) %>%
filter(facility_county %in% c("Bronx", "Kings", "New York", "Queens", "Richmond")) %>%
filter(!((facility_longitude == 0) & (facility_latitude == 0))) %>%
filter(facility_id != 5764) %>%
st_as_sf(coords = c("facility_longitude", "facility_latitude"), crs = 2263)
# check some stuff
# str(health_facility_sf)
# st_crs(health_facility_sf)
# I want to get the number and quantity of abbreviations
# abbr_counts <- health_facility_sf %>%
# count(short_description, name = "count") %>%
# arrange(desc(count))
# print(abbr_counts)
# Based on this I can select what I want to keep by filtering out certain short_description entries
filtered_health_facility_sf <- health_facility_sf %>%
filter(!short_description %in% c("DTC", "DTC-EC", "CHHA", "LTHHCP"))
# check using the previous code
abbr_counts_check <- filtered_health_facility_sf %>%
count(short_description, name = "count") %>%
arrange(desc(count))
print(abbr_counts_check)
# now we can do another aggregation via the same methodology
health_facility_query_sf <- zip_fstore_covid_sf %>%
st_join(filtered_health_facility_sf, join = st_within) %>%
group_by(ZIPCODE) %>%
summarize(health_facility_count = n(),
geometry = st_union(geometry))
health_facility_query_df <- as.data.frame(health_facility_query_sf)
health_facility_query_df[3] <- NULL
zip_covid_fstore_hfacility_agg_sf <- zip_fstore_covid_sf %>%
left_join(health_facility_query_df, by = 'ZIPCODE')
# now to join ACS to census
# complicated to join so we need to clean first
# I look for the strings of county names, converting where needed, and generate a new column; boro
acs_demographics_housing_clean_df <- acs_demographics_housing_df %>%
mutate(tract_number = str_extract(geographic_area_name, "\\d+")) %>%
mutate(
boro = case_when(
str_detect(`geographic_area_name`, "Bronx") ~ "Bronx",
str_detect(`geographic_area_name`, "Queens") ~ "Queens",
str_detect(`geographic_area_name`, "Kings") ~ "Brooklyn",
str_detect(`geographic_area_name`, "New York County") ~ "Manhattan",
str_detect(`geographic_area_name`, "Richmond") ~ "Staten Island"
)
) %>%
select(id, tract_number, boro, everything())
# now we can join
acs_demographics_housing_sf <- census_tracts_2010_sf %>%
left_join(acs_demographics_housing_clean_df, by = c("boro_name" = "boro", "ctlabel" = "tract_number"))Here I did the last few joins to create the final shapefile.
# almost there, just need to aggregate the census and zip code stuff, another big one
acs_demographics_housing_sf <- st_transform(acs_demographics_housing_sf, st_crs(zip_covid_fstore_hfacility_agg_sf))
joined_sf <- st_join(zip_covid_fstore_hfacility_agg_sf, acs_demographics_housing_sf, join = st_intersects)
# list_of_columns <- names(joined_sf)
# print(list_of_columns)
zip_fstores_hfacilities_demographics_sf <- joined_sf %>%
group_by(ZIPCODE) %>%
mutate(across(c(40:395), ~if(is.numeric(.x)) sum(.x, na.rm = TRUE) else NA_real_, .names = "{.col}_sum")) %>%
ungroup()
# Now to remove columns, first the unsummed columns and old columns from the census shape
zip_fstores_hfacilities_demographics_sf <- zip_fstores_hfacilities_demographics_sf %>%
select(1:26, 396:ncol(zip_fstores_hfacilities_demographics_sf))
# here I want to remove any variable (column) namse that include percent_estimate or margins of error as they are unneccessary
zip_fstores_hfacilities_demographics_sf <- zip_fstores_hfacilities_demographics_sf %>%
select(-matches("percent_estimate|margin"))
# I use this to determine what columns to keep or remove
list_of_columns <- names(zip_fstores_hfacilities_demographics_sf)
print(list_of_columns)