Plot at least two high-quality static maps with one using the
COVID-19 data and one using a related factor. You can use either plot
method for sf or ggplot method.
Initial Data Preparation (adapted from Week 08 Assignment):
#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)
#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.
#Merged
covidBYzip <- base::merge(nycZipCodes_sf, covid19, by.x="ZIPCODE", by.y="MODIFIED_ZCTA")
#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) %>%
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=''))
## 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)
#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))
popData_merged <- merge(nycCensus, acsData, by.x ='tractFIPS', by.y = 'censusCode') %>%
sf::st_transform(st_crs(covidBYzip))
# Aggregate ACS/Census data to zip code areas and COVID-19 data
popData_nyc <- sf::st_join(covidBYzip,
popData_merged %>% sf::st_centroid(),
join = st_contains) %>%
group_by(ZIPCODE, PO_NAME, POPULATION, COUNTY, COVID_CASE_COUNT, TOTAL_COVID_TESTS) %>%
summarise(
totPop = sum(totPop, na.rm = TRUE),
malePctg = sum(malePop, na.rm = TRUE)/totPop*100,
asianPop = sum(asianPop, na.rm = TRUE),
blackPop = sum(blackPop, na.rm = TRUE),
hispanicPop = sum(hispanicPop, na.rm = TRUE),
whitePop = sum(whitePop, na.rm = TRUE),
elderlyPop = sum(elderlyPop, na.rm = TRUE),
COVID_CASE_COUNT = sum(COVID_CASE_COUNT, na.rm = TRUE))
## 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.
Map #1:
plot(covidBYzip["TOTAL_COVID_TESTS"],
main = "Total COVID-19 Tests; April 23, 2020",
breaks = "jenks",
graticule = st_crs(4326),
axes=TRUE,
reset=FALSE)

Map #2:
plot(covidBYzip["PERCENT_POSITIVE"],
main = "Percentage of COVID-19 Tests Reporting Positive; April 23, 2020",
breaks = "jenks",
graticule = st_crs(4326),
axes=TRUE,
reset=FALSE)

Map #3:
plot(popData_nyc["elderlyPop"],
breaks='jenks',
main = "Elderly Population (age 65+) by Zip Code")
