1. High quality static maps

tract_2_zip <- st_transform(tract_2_zip,4326)
pal <- brewer.pal(5,"OrRd")
pal_1 <- brewer.pal(5,"Greens")
options(scipen = 999)
pal <- brewer.pal(5,"OrRd")
pal_1 <- brewer.pal(5,"Greens")

plot(tract_2_zip['total_tests'],
     main = "Total Covid-19 Tests Administered in NYC",
     breaks = "jenks", nbreaks = 5,
     pal = pal,
     graticule = st_crs(4326),
     axes = TRUE)

plot(tract_2_zip['adultPop'],
     main = "Adult Population in NYC",
     breaks = "jenks", nbreaks=5,
     pal = pal_1,
     graticule = st_crs(4326),
     axes = TRUE)

2. Create a multimap figure

breaks_qt <- classIntervals(c(min(tract_2_zip$case_count) -.00001,
                             tract_2_zip$case_count),n = 7, style = "jenks")
breaks_qt_1 <- classIntervals(c(min(tract_2_zip$elderlyPop) -.00001,
                              tract_2_zip$elderlyPop),n = 7, style = "jenks")
tract_2_zip <- mutate(tract_2_zip,case_count_cat = cut(case_count,breaks_qt$brks))
tract_2_zip <- mutate(tract_2_zip,elderlyPop_cat = cut(elderlyPop,breaks_qt_1$brks))
p1 <- ggplot(tract_2_zip) +
  geom_sf(aes(fill = case_count)) + 
  geom_sf_label(data = tract_2_zip %>% filter(case_count > 10000),
                aes(label = case_count),
                label.size = .09,
                size = 3) +
  scale_fill_gradient(low = "orange", high = "red", name = "Rate",labels = comma) +
  coord_sf(default_crs = sf::st_crs(4326)) +
  labs(x='Longitude',y='Latitude',title = "Case Count") + 
  theme(plot.title = element_text(hjust = 0.5),axis.text.x = element_text(size = 7.5),
        legend.position = "bottom",
        legend.key.size = unit(0.3, 'cm'), 
        legend.key.height = unit(0.3, 'cm'), 
        legend.key.width = unit(0.5, 'cm'), 
        legend.title = element_text(size=7), 
        legend.text = element_text(size=5),
        panel.background = element_rect(fill = "white",
        colour = "white", size = 0.5, linetype = "solid"),
        panel.grid.major = element_line(size = 0.5, linetype = 'solid',colour = "lightgrey"), 
  panel.grid.minor = element_line(size = 0.25, linetype = 'solid',colour = "lightgrey"))

p2 <- ggplot(tract_2_zip) + geom_sf(aes(fill = elderlyPop)) + 
  geom_sf_label(data = tract_2_zip %>% filter(elderlyPop > 15000),
                aes(label = elderlyPop),
                label.size = .08,
                size = 3) +
  scale_fill_gradient(low = "lightblue", high = "blue", name = "Populaiton",labels = comma) +
  coord_sf(default_crs = sf::st_crs(4326)) + 
  labs(x = 'Longitude', y = 'Latitude', title = " Elderly Population") +
  theme(plot.title = element_text(hjust=0.5),axis.text.x = element_text(size = 7.5),
        legend.position = "bottom",
        legend.key.size = unit(0.3, 'cm'), 
        legend.key.height = unit(0.3, 'cm'), 
        legend.key.width = unit(0.5, 'cm'), 
        legend.title = element_text(size=7), 
        legend.text = element_text(size=5),
        panel.background = element_rect(fill = "white",
        colour = "white", size = 0.5, linetype = "solid"),
        panel.grid.major = element_line(size = 0.5, linetype = 'solid',colour = "lightgrey"), 
  panel.grid.minor = element_line(size = 0.25, linetype = 'solid',colour = "lightgrey")) 
grid.arrange(p1,p2,ncol = 2,
  top=textGrob("Covid-19 Cases and Elderly Population in NYC",gp=gpar(fontsize=15,font=8)))

3. Create a web based interactive map

pal_fun <- colorQuantile("YlOrRd",NULL,n=5)
p_tip <-paste0("GEOID: ", tract_2_zip$ZIPCODE);

p_popup <- paste0("<strong>Covid-19 Case Rate: </strong>",
                  tract_2_zip$case_rate%>%round(3)%>%format(nsmall=3),
                  "<br/>",
                  "<strong>Number of Covid-19 Cases: </strong>",
                  tract_2_zip$case_count,
                  sep = "")
polyHighlightOption <- leaflet::highlightOptions(opacity = 1.0, fillColor = 'black')
polyLabelOption <- leaflet::labelOptions(opacity = 0.6)

c19_html_map <- leaflet(tract_2_zip) %>%
  addPolygons(
    stroke = FALSE,
    fillColor = ~pal_fun(case_rate),
    fillOpacity = 0.8, smoothFactor =0.5,
    popup = p_popup,
    label = p_tip,
    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 = "Covid-19 <br> Case Rate") %>% 
  addLayersControl(baseGroups = c("OSM","Carto"),
                   overlayGroups = c("nyc"))
c19_html_map
htmlwidgets::saveWidget(c19_html_map, 'nyc_CovidCaseRate_leaflet_widget.html')