R-Spatial Week 08 Assignment
Ella Bayer
2025-03-29
Week 08 Assignment Task Overview
- Join the COVID-19 data to the NYC zip code area data (sf or sp polygons).
- Aggregate the NYC food retails store data (points) to the zip code data, so that we know how many retail stores are 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.
- Aggregate the NYC health facilities (points) to the zip code data. Similarly, choose appropriate subtypes such as nursing homes from the facilities.
- Join the Census ACS population, race, and age data to the NYC Planning Census Tract Data.
- Aggregate the ACS census data to zip code area data.
In the end, we should have the confirmed and tested cases of covid-19, numbers of specific types of food stores, numbers of specific types of health facilities, and population (total population, elderly, by race, etc.) at the zip code level. We should also have boroughs, names, etc. for each zip code area.
Task 1:
Join the COVID-19 data to the NYC zip code area data.
#COVID-19 data
covid19 <- readr::read_csv("/Users/ellabayer/Documents/Hunter College/DataVisR/R-spatial/R-Spatial_II_Lab/tests-by-zcta_2021_04_23.csv")## Rows: 177 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): NEIGHBORHOOD_NAME, BOROUGH_GROUP, label
## dbl (10): MODIFIED_ZCTA, lat, lon, COVID_CASE_COUNT, COVID_CASE_RATE, POP_DE...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#NYC Zip Codes
nycZipCodes_sf <- sf::st_read("/Users/ellabayer/Documents/Hunter College/DataVisR/R-spatial/R-Spatial_I_Lab/ZIP_CODE_040114/ZIP_CODE_040114.shp") %>%
st_transform(crs = 4326)## Reading layer `ZIP_CODE_040114' from data source
## `/Users/ellabayer/Documents/Hunter College/DataVisR/R-spatial/R-Spatial_I_Lab/ZIP_CODE_040114/ZIP_CODE_040114.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 263 features and 12 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 913129 ymin: 120020.9 xmax: 1067494 ymax: 272710.9
## Projected CRS: NAD83 / New York Long Island (ftUS)#Merged
covid_to_zip <- base::merge(nycZipCodes_sf, covid19, by.x="ZIPCODE", by.y="MODIFIED_ZCTA")Task 2:
Aggregate the NYC food retails store data (points) to the zip code data, so that we know how many retail stores are in each zip code area.
nyc_retail <- st_read("/Users/ellabayer/Documents/Hunter College/DataVisR/R-spatial/R-Spatial_II_Lab/nycFoodStore.shp")## Reading layer `nycFoodStore' from data source
## `/Users/ellabayer/Documents/Hunter College/DataVisR/R-spatial/R-Spatial_II_Lab/nycFoodStore.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 11300 features and 16 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -74.2484 ymin: 40.50782 xmax: -73.67061 ymax: 40.91008
## Geodetic CRS: WGS 84nyc_food <- nyc_retail %>%
dplyr::filter(Estbl_T != "A") %>% #any Establishment type OTHER THAN "A" is a food retail store
sf::st_join(nycZipCodes_sf, ., join= st_contains) %>%
group_by(ZIPCODE) %>%
summarise(number_food = n())Visual Verification of Data:
plot(nyc_food["number_food"], breaks = "jenks", main = "Number of Food Retail 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.
nys_health <- read.csv("/Users/ellabayer/Documents/Hunter College/DataVisR/R-spatial/R-Spatial_II_Lab/NYS_Health_Facility.csv")
nyc_NursHomes <- nys_health %>%
dplyr::filter(Facility.County == c('Bronx', 'Kings', 'New York', 'Queens', 'Richmond')) %>% #filter for NYC counties
dplyr::filter(Short.Description == 'NH') %>% #filter for Nursing Homes
sf::st_as_sf(coords = c("Facility.Longitude", "Facility.Latitude"), crs = 4326) %>% #convert to sf object
sf::st_join(nycZipCodes_sf, ., join= st_contains) #join to zip code data
#Number of nursing homes by zip code
nyc_NursHomes_COUNT <- nyc_NursHomes %>%
group_by(ZIPCODE) %>%
summarise(num_NH = n())Visual Verification of Data:
plot(nyc_NursHomes_COUNT["num_NH"], breaks = "jenks", main = "Number of Nursing Homes by Zip Code")
Task 4:
Join the Census ACS population, race, and age data to the NYC Planning Census Tract Data.
#Import Census Tracts
nycCensus <- sf::st_read("/Users/ellabayer/Documents/Hunter College/DataVisR/R-spatial/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 `/Users/ellabayer/Documents/Hunter College/DataVisR/R-spatial/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)nycCensus %<>% 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='')
)#Import ACS Data
acsData<-readLines("/Users/ellabayer/Documents/Hunter College/DataVisR/R-spatial/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, # >= 65
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))Check first 10 rows to ensure it looks as expected:
acsData %>%
magrittr::extract(1:10,)## GEO_ID totPop elderlyPop malePop femalePop whitePop blackPop
## 1 1400000US36005000100 7080 51 6503 577 1773 4239
## 2 1400000US36005000200 4542 950 2264 2278 2165 1279
## 3 1400000US36005000400 5634 710 2807 2827 2623 1699
## 4 1400000US36005001600 5917 989 2365 3552 2406 2434
## 5 1400000US36005001900 2765 76 1363 1402 585 1041
## 6 1400000US36005002000 9409 977 4119 5290 3185 4487
## 7 1400000US36005002300 4600 648 2175 2425 479 2122
## 8 1400000US36005002400 172 0 121 51 69 89
## 9 1400000US36005002500 5887 548 2958 2929 903 1344
## 10 1400000US36005002701 2868 243 1259 1609 243 987
## asianPop hispanicPop adultPop citizenAdult censusCode
## 1 130 2329 6909 6100 005000100
## 2 119 3367 3582 2952 005000200
## 3 226 3873 4507 4214 005000400
## 4 68 3603 4416 3851 005001600
## 5 130 1413 2008 1787 005001900
## 6 29 5905 6851 6170 005002000
## 7 27 2674 3498 3056 005002300
## 8 14 0 131 42 005002400
## 9 68 4562 4237 2722 005002500
## 10 0 1985 1848 1412 005002701# merge Census and ACS data
popData_merged <- merge(nycCensus, acsData, by.x ='tractFIPS', by.y = 'censusCode')Task 5:
Aggregate the ACS census data to zip code area data.
#Set Census/ACS CRS to same as Zip Code CRS
popData_nyc <- sf::st_transform(popData_merged, st_crs(covid_to_zip))# Aggregate ACS/Census data to zip code areas
popData_nyc_to_zip <- sf::st_join(covid_to_zip,
popData_nyc %>% sf::st_centroid(),
join = st_contains) %>%
group_by(ZIPCODE, PO_NAME, POPULATION, COUNTY, COVID_CASE_COUNT, TOTAL_COVID_TESTS) %>%
summarise(totPop = sum(totPop),
malePctg = sum(malePop)/totPop*100,
asianPop = sum(asianPop),
blackPop = sum(blackPop),
hispanicPop = sum(hispanicPop),
whitePop = sum(whitePop))## Warning: st_centroid assumes attributes are constant over geometries## `summarise()` has grouped output by 'ZIPCODE', 'PO_NAME', 'POPULATION',
## 'COUNTY', 'COVID_CASE_COUNT'. You can override using the `.groups` argument.Visual Verification of Data (Covid):
plot(popData_nyc_to_zip["COVID_CASE_COUNT"], breaks='jenks', main = "Covid Cases by Zip Code")
Visual Verification of Data (Population):
plot(popData_nyc_to_zip["totPop"], breaks='jenks', main = "Total Population by Zip Code")