3.7 Lab Assignment

Task 1

Plot at least two high-quality static maps at the zip code level, with one using the COVID-19 testing data and one using a related factor such as the number of nursing homes or the elderly population. You can use either plot method for sf or ggplot method.

acsCovidDat <- sf::st_read('./data/acsPopByZip.gpkg', layer="NYC_COVID_Pop_by_Zip")
## Reading layer `NYC_COVID_Pop_by_Zip' from data source 
##   `/Users/staiceyandreauy/Desktop/GTECH 38520/Assignments/R-Spatial/data/acsPopByZip.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 180 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 913129 ymin: 120020.9 xmax: 1067113 ymax: 272710.9
## Projected CRS: NAD83 / New York Long Island (ftUS)
breaks_qt <- classIntervals(c(min(acsCovidDat$Positiv) - .00001,
                              acsCovidDat$Positiv), n = 7, style = "quantile")

breaks_qt
## style: quantile
##   one of 28,530,983,404 possible partitions of this variable into 7 classes
##      [13.99999,256)      [256,377.8571) [377.8571,490.8571) [490.8571,709.4286) 
##                  25                  27                  26                  25 
## [709.4286,999.7143) [999.7143,1302.143)     [1302.143,2817] 
##                  26                  26                  26
acsCovidDat <- mutate(acsCovidDat, 
                            acsCovidDat_cat = cut(Positiv, breaks_qt$brks, 
                                                 dig.lab = 4, digits=1)) 

ggplot(acsCovidDat) + 
  geom_sf(aes(fill=acsCovidDat_cat)) +
  scale_fill_brewer(palette = "OrRd", name='COVID Cases') +
  coord_sf(default_crs = sf::st_crs(4326) ) +
  labs(x='Longitude', y='Latitude', 
       title='NYC Covid-19 Cases by Zipcode',
       caption = 'Data Source: Census ACS Population Data') +
  theme(legend.position = "right")

breaks_qt_elder <- classIntervals(
  c(min(acsCovidDat$eldrlyP) - 0.00001, acsCovidDat$eldrlyP),
  n = 7,
  style = "quantile"
)
## Warning in classIntervals(c(min(acsCovidDat$eldrlyP) - 1e-05,
## acsCovidDat$eldrlyP), : var has missing values, omitted in finding classes
acsCovidDat <- acsCovidDat %>%
  mutate(eldrly_cat = cut(eldrlyP, breaks_qt_elder$brks, dig.lab = 4, digits = 1))

ggplot(acsCovidDat) +
  geom_sf(aes(fill = eldrly_cat)) +
  scale_fill_brewer(palette = "Blues", name = "Percentage of Elderly Population") +
  coord_sf(default_crs = sf::st_crs(4326)) +
  labs(
    x = "Longitude",
    y = "Latitude",
    title = "NYC Elderly Population by ZIP Code",
    caption = "Data Source: Census ACS Population Data"
  ) +
  theme(legend.position = "right")

Task 2

Use ggplot2 and other ggplot-compatible packages to create a multi-map figure illustrating the possible relationship between COVID-19 confirmed cases or rate and another factor (e.g., the number of nursing homes, neighborhood racial composition, elderly population, etc.).

require(classInt)
require(patchwork)
## Loading required package: patchwork
g_1 <- ggplot(acsCovidDat) +
  geom_sf(aes(fill=acsCovidDat_cat)) +
  scale_fill_brewer(palette="OrRd", name="NYC COVID-19 Cases") +
  coord_sf(default_crs=sf::st_crs(4326)) +
  labs(x="Longitude", y="Latitude",
       title="NYC COVID-19 Cases",
       caption="Data Source: Census ACS Population Data") +
  theme(legend.position="left")

g_2 <- ggplot(acsCovidDat) +
  geom_sf(aes(fill=eldrly_cat)) +
  scale_fill_brewer(palette="Blues", name="% of Elderly Pop") +
  coord_sf(default_crs=sf::st_crs(4326)) +
  labs(x="Longitude", y="Latitude",
       title="NYC Elderly Population",
       caption="Data Source: Census ACS Population Data") +
  theme(legend.position="bottom")

g_1 + g_2 + plot_layout(ncol=2)

Task 3

Create a web-based interactive map for COIVD-19 data using tmap, mapview, or leaflet package and save it as a HTML file.

library(leaflet) 

st_crs(acsCovidDat)
## Coordinate Reference System:
##   User input: NAD83 / New York Long Island (ftUS) 
##   wkt:
## PROJCRS["NAD83 / New York Long Island (ftUS)",
##     BASEGEOGCRS["NAD83",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["SPCS83 New York Long Island zone (US survey foot)",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",40.1666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-74,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",41.0333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",40.6666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",984250,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["United States (USA) - New York - counties of Bronx; Kings; Nassau; New York; Queens; Richmond; Suffolk."],
##         BBOX[40.47,-74.26,41.3,-71.8]],
##     ID["EPSG",2263]]
pal_fun <- colorQuantile("YlOrRd", NULL, n = 5)

pal_elderly <- colorQuantile("YlGnBu", 
                        domain = c(min(acsCovidDat$eldrlyP, na.rm = T), max(acsCovidDat$eldrlyP, na.rm = T)),
                        alpha = TRUE)

covid_popup <- paste0("<strong>Zip Code: </strong>", 
                  acsCovidDat$ZIPCODE,
                   " <br/>",
                  "<strong> Place Name: </strong>",
                  acsCovidDat$PO_NAME,
                  " <br/>",
                  "<strong> Confirmed COVID-19 Cases: </strong>",
                  acsCovidDat$Positiv,
                  sep="")

polyHighlightOption <- leaflet::highlightOptions(opacity = 1.0, fillColor = 'black')
polyLabelOption <- leaflet::labelOptions(opacity = 1.0, textsize = '11px')

htmlMap <- leaflet(acsCovidDat %>% sf::st_transform(4326)) %>%
        addPolygons(
                stroke = TRUE, # show the polygon boundary
                fillOpacity = 0.8,
                fillColor = ~pal_elderly(eldrlyP),
                color = 'grey', # color of the boundary
                weight = 1, # width of the boundary
                highlightOptions = polyHighlightOption,
                label = paste(acsCovidDat$eldrlyP),
                group = "Elderly Population") %>%
        addPolygons(
                stroke = FALSE,
                fillColor = ~pal_fun(Positiv),
                fillOpacity = 0.8, smoothFactor = 0.5,
                popup = covid_popup,
                group = "NYC-COVID",
                label = paste0(acsCovidDat$Positiv, " at ", acsCovidDat$PO_NAME, ", ", acsCovidDat$ZIPCODE),
                highlightOptions = polyHighlightOption,
                labelOptions = polyLabelOption) %>%
        addTiles(group = "OSM") %>%
        addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
  addLegend("bottomright", 
            group = "NYC",
            colors = brewer.pal(7, "YlOrRd"), 
            labels = paste0("up to ", format(breaks_qt$brks[-1], digits = 2)),
            title = "NYC <br> COVID Cases") %>%
        addLayersControl(baseGroups = c("OSM", "Carto"), 
                         overlayGroups = c("NYC-COVID", "Elderly Population"))  

htmlMap
htmlwidgets::saveWidget(htmlMap, 'NYC_Covid_&_ElderlyPop.html', 
                        selfcontained = FALSE, libdir = 'widget-lib')