R Spaital Lab Assignment # 2

task 1: Join the COVID-19 data to the NYC zip code area data (sf or sp polygons).

# load in RData file from assignment #1 with processed datasets as sf objects
load("R-Spatial_II_Lab/assignment1_data.RData")

# read in covid test data
covid_tests_04_12_2020 <- readr::read_csv("R-Spatial_II_Lab/tests-by-zcta_2020_04_12.csv", lazy = FALSE, show_col_types = FALSE)
colnames(covid_tests_04_12_2020) <- c("zipcode", "posTests_04_12_2020", "totTests__04_12_2020", "pcPosTests_04_12_2020")

# join covid data with zip code sf
sf_nyc_zip_merged <- base::merge(nyc_zip_sf, covid_tests_04_12_2020, by.x = "ZIPCODE", by.y = "zipcode")

# map merged sf to view and confirm
mapview::mapview(sf_nyc_zip_merged, zcol = 'posTests_04_12_2020', layer.name = "Positive COVID Tests 04-12-20 by Zip Code")

task 2: Aggregate the NYC food retails store data (points) to the zip code data, so that we know how many retail stores in each zip code area. Note that not all locations are for food retail. And we need to choose the specific types according to the data.

# filter food store type and count the number contained per zip code
nyc_tot_foodStores_zc <- nyc_rfs %>% dplyr::filter(stringr::str_detect(establishment_type, '[AJD]')) %>%
  sf::st_join(sf_nyc_zip_merged, ., join = st_contains) %>%
  group_by(ZIPCODE) %>%
  summarise(tot_foodStores = n()) %>%
  st_drop_geometry()

# join total number of food stores by zip code to zip code sf
sf_nyc_zip_merged <- merge(sf_nyc_zip_merged, nyc_tot_foodStores_zc, by.x = "ZIPCODE", by.y = "ZIPCODE")

# map total food stores by zip code
mapview::mapview(sf_nyc_zip_merged, zcol = 'tot_foodStores', layer.name = "Total Food Stores by Zip Code")

task 3: Aggregate the NYC health facilities (points) to the zip code data. Similarly, choose appropriate subtypes such as nursing homes from the facilities.

nyc_tot_nursingHomes_zc <- nys_hf %>% 
  dplyr::filter(short_description == 'NH') %>%
  sf::st_join(sf_nyc_zip_merged, ., join = st_contains) %>%
  group_by(ZIPCODE) %>%
  summarise(tot_nursingHomes = n()) %>%
  st_drop_geometry()

# join total number of nursing homes by zip code to zip code sf
sf_nyc_zip_merged <- merge(sf_nyc_zip_merged, nyc_tot_nursingHomes_zc, by.x = "ZIPCODE", by.y = "ZIPCODE")

# map total food stores by zip code
mapview::mapview(sf_nyc_zip_merged, zcol = 'tot_nursingHomes', layer.name = "Total Nursing Homes by Zip Code")

task 4: Join the Census ACS population, race, and age data to the NYC Planning Census Tract Data.

# read in census tracts and make sure its in the correct crs
nyc_tracts <- sf::st_read('R-Spatial_II_Lab/2010 Census Tracts/geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a.shp', stringsAsFactors = FALSE)
## Reading layer `geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a' from data source `C:\Users\jdhoc\Desktop\Hunter\SPRING 2024\GEOG 39300 R\R-spatial\Session_8\R-Spatial_II_Lab\2010 Census Tracts\geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2165 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91553
## Geodetic CRS:  WGS84(DD)
nyc_tracts <- st_transform(nyc_tracts, 4326)

# create new columns for county and tracts FIPS codes based on boroughs
nyc_tracts %<>% dplyr::mutate(cntyFIPS = case_when(
  boro_name == 'Bronx' ~ '005',
  boro_name == 'Brooklyn' ~ '047',
  boro_name == 'Manhattan' ~ '061',
  boro_name == 'Queens' ~ '081',
  boro_name == 'Staten Island' ~ '085'),
  tractFIPS = paste(cntyFIPS, ct2010, sep='')
)

# read in acs data, exclude acs column names, and filter for desired population stats
acs <- readLines("./R-Spatial_II_Lab/ACSDP5Y2018.DP05_data_with_overlays_2020-04-22T132935.csv") %>%
  magrittr::extract(-2) %>% 
  textConnection() %>%
  read.csv(header=TRUE, quote= "\"") %>%
  dplyr::select(GEO_ID, 
                totPop = DP05_0001E, elderlyPop = DP05_0024E,
                malePop = DP05_0002E, femalePop = DP05_0003E,  
                whitePop = DP05_0037E, blackPop = DP05_0038E,
                asianPop = DP05_0067E, hispanicPop = DP05_0071E,
                adultPop = DP05_0021E, citizenAdult = DP05_0087E) %>%
  dplyr::mutate(censusCode = stringr::str_sub(GEO_ID, -9,-1))

# join tracts and acs data
sf_nyc_tracts_merged <- merge(nyc_tracts, acs, by.x ='tractFIPS', by.y = 'censusCode')

# map the tracts and an acs pop stat to view and test
mapview::mapview(sf_nyc_tracts_merged, zcol = 'elderlyPop', layer.name = "Elderly (65+) Population by Census Tract")

task 5: Aggregate the ACS census data to zip code area data.

# aggregate acs stats by zip code
acs_zip <- sf::st_join(sf_nyc_zip_merged, 
                    sf_nyc_tracts_merged %>% sf::st_centroid(),
                    join = st_contains) %>% 
  group_by(ZIPCODE) %>%
  summarise(totPop = sum(totPop), elderlyPop = sum(elderlyPop),
            malePop = sum(malePop), femalePop = sum(femalePop), 
            asianPop = sum(asianPop), blackPop = sum(blackPop), 
            hispanicPop = sum(hispanicPop), whitePop = sum(whitePop), 
            adultPop = sum(adultPop), adultPctg = sum(adultPop)/totPop*100,
            citizenAdult = sum(citizenAdult))
## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x
# merge aggregated stats to zip code sf
acs_zip <- st_drop_geometry(acs_zip)
sf_nyc_zip_merged <- merge(sf_nyc_zip_merged, acs_zip, by.x = "ZIPCODE", by.y = "ZIPCODE")

# map to view and test
mapview::mapview(sf_nyc_zip_merged, zcol = 'adultPctg', layer.name = "Adult Percent of Total Population by Zip Code")