Load packages

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

#Using annotation to enhance map view
#Map for COVID-19 cases
map_covid <- ggplot() +
  geom_sf(data = acsCovidDat, aes(fill = Positiv), color = "black", size = 0.2) +
  scale_fill_viridis_c(name = "COVID-19 Cases") +
  theme_minimal() +
  ggtitle("COVID-19 Positive Cases by ZIP Code") +
  theme(legend.position = "bottom") +
  annotation_scale(location = "br") + 
  annotation_north_arrow(location = "tl", which_north = "true", style = north_arrow_fancy_orienteering())

#Map for Elderly population
map_elderly <- ggplot() +
  geom_sf(data = acsCovidDat, aes(fill = eldrlyP), color = "black", size = 0.2) +
  scale_fill_viridis_c(name = "Elderly Population") +
  theme_minimal() +
  ggtitle("Elderly Population by ZIP Code") +
  theme(legend.position = "bottom") +
  annotation_scale(location = "br") +  
  annotation_north_arrow(location = "tl", which_north = "true", style = north_arrow_fancy_orienteering()) +
  coord_sf(datum = NA)

#Combine the two maps side by side using patchwork Reference -> https://stackoverflow.com/questions/1249548/side-by-side-plots-with-ggplot2 
final_plot <- map_covid + map_elderly + plot_layout(ncol = 2)

#View plot
print(final_plot)

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

#Color palettes for COVID data and elderly population
pal_fun <- colorQuantile("YlOrRd", NULL, n = 5)
pal_elderly <- colorQuantile("YlGnBu", 
                        domain = c(min(acsCovidDat$eldrlyP, na.rm = TRUE), max(acsCovidDat$eldrlyP, na.rm = TRUE)),
                        alpha = TRUE)

#Popup with information about each ZIP code
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="")

#Highlight options when hovered over them
polyHighlightOption <- leaflet::highlightOptions(opacity = 1.0, fillColor = 'black')
polyLabelOption <- leaflet::labelOptions(opacity = 1.0, textsize = '11px')

#Create the interactive map
htmlMap <- leaflet(acsCovidDat %>% sf::st_transform(4326)) %>%
        addPolygons(
                stroke = TRUE, 
                fillOpacity = 0.8,
                fillColor = ~pal_elderly(eldrlyP),
                color = 'grey', 
                weight = 1, 
                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") %>%
        addLayersControl(baseGroups = c("OSM", "Carto"), 
                         overlayGroups = c("NYC-COVID", "Elderly Population"))  

#View the map
htmlMap
#Save the map as an HTML file
htmlwidgets::saveWidget(htmlMap, "NYC_COVID_Map.html")