acsCovidDat <- acsCovidDat %>%
  mutate(covid_rate = (Positiv / totPop) * 10000)

acsCovidDat <- st_transform(acsCovidDat, 4326)
breaks_covid <- classIntervals(c(min(acsCovidDat$covid_rate, na.rm=TRUE) - 0.00001,
                                  acsCovidDat$covid_rate),
                                n = 6, style = "quantile", na.rm = TRUE)

acsCovidDat <- acsCovidDat %>%
  mutate(covid_rate_cat = cut(covid_rate, breaks_covid$brks, dig.lab = 4))

breaks_eld <- classIntervals(c(min(acsCovidDat$eldrlyP, na.rm=TRUE) - 0.00001,
                                acsCovidDat$eldrlyP),
                              n = 6, style = "quantile", na.rm = TRUE)

acsCovidDat <- acsCovidDat %>%
  mutate(elderly_cat = cut(eldrlyP, breaks_eld$brks,
                           labels = c("0–2,876", "2,876–4,265", "4,265–6,160",
                                      "6,160–7,737", "7,737–10,040", "10,040–16,570")))
mapA <- ggplot(acsCovidDat) +
  geom_sf(aes(fill = covid_rate_cat), color = "white", linewidth = 0.2) +
  scale_fill_brewer(palette = "YlOrRd", name = "Cases per\n10,000 people",
                    na.value = "grey80") +
  coord_sf(crs = st_crs(4326)) +
  labs(
    title = "NYC COVID-19 Confirmed Case Rate by ZIP Code",
    caption = "Source: NYC Health, US Census ACS",
    x = "Longitude",
    y = "Latitude"
  ) +
  annotation_scale(location = "bl", width_hint = 0.3) +
  annotation_north_arrow(location = "tr", which_north = "true",
                         style = north_arrow_fancy_orienteering()) +
  theme_bw() +
  theme(
    panel.border = element_rect(color = "black", linewidth = 1),
    legend.position = "right"
  )

mapA

mapB <- ggplot(acsCovidDat) +
  geom_sf(aes(fill = elderly_cat), color = "white", linewidth = 0.2) +
  scale_fill_brewer(palette = "YlGnBu", name = "Elderly\nPopulation",
                    na.value = "grey80") +
  coord_sf(crs = st_crs(4326)) +
  labs(
    title = "NYC Elderly Population by ZIP Code",
    caption = "Source: US Census ACS",
    x = "Longitude",
    y = "Latitude"
  ) +
  annotation_scale(location = "bl", width_hint = 0.3) +
  annotation_north_arrow(location = "tr", which_north = "true",
                         style = north_arrow_fancy_orienteering()) +
  theme_bw() +
  theme(
    panel.border = element_rect(color = "black", linewidth = 1),
    legend.position = "right"
  )

mapB

library(patchwork)

mapA_side <- mapA + 
  labs(title = "COVID-19 Case Rate") +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 7),
        legend.title = element_text(size = 8),
        plot.title = element_text(size = 14, hjust = 0.5))

mapB_side <- mapB + 
  labs(title = "Elderly Population") +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 7),
        legend.title = element_text(size = 8),
        plot.title = element_text(size = 14, hjust = 0.5))

mapA_side + plot_spacer() + mapB_side + plot_layout(ncol = 3, widths = c(10, 1, 10))

pal_covid <- colorQuantile("YlOrRd", domain = acsCovidDat$covid_rate, n = 5)
pal_elderly <- colorQuantile("YlGnBu", domain = acsCovidDat$eldrlyP, n = 5)
pal_hispanic <- colorQuantile("PuRd", domain = acsCovidDat$hspncPp, n = 5)

covid_popup <- paste0("<strong>ZIP Code: </strong>", acsCovidDat$ZIPCODE,
                      "<br/><strong>Neighborhood: </strong>", acsCovidDat$PO_NAME,
                      "<br/><strong>COVID Rate: </strong>", round(acsCovidDat$covid_rate, 1),
                      " per 10,000")


elderly_popup <- paste0("<strong>ZIP Code: </strong>", acsCovidDat$ZIPCODE,
                        "<br/><strong>Elderly Population: </strong>", acsCovidDat$eldrlyP)

hispanic_popup <- paste0("<strong>ZIP Code: </strong>", acsCovidDat$ZIPCODE,
                         "<br/><strong>Hispanic Population: </strong>", acsCovidDat$hspncPp)

htmlMap <- leaflet(acsCovidDat %>% sf::st_transform(4326)) %>%
  addPolygons(fillColor = ~pal_covid(covid_rate),
              fillOpacity = 0.8, stroke = TRUE, color = "white", weight = 0.5,
              label = ~paste0(PO_NAME, ": ", round(covid_rate, 1), " per 10,000"),
              popup = covid_popup,
              highlightOptions = highlightOptions(opacity = 1.0, fillColor = "black"),
              group = "COVID Rate") %>%
  addPolygons(fillColor = ~pal_elderly(eldrlyP),
              fillOpacity = 0.8, stroke = TRUE, color = "white", weight = 0.5,
              label = ~paste0(PO_NAME, ": ", eldrlyP, " elderly residents"),
              popup = elderly_popup,
              highlightOptions = highlightOptions(opacity = 1.0, fillColor = "black"),
              group = "Elderly Population") %>%
  addPolygons(fillColor = ~pal_hispanic(hspncPp),
              fillOpacity = 0.8, stroke = TRUE, color = "white", weight = 0.5,
              label = ~paste0(PO_NAME, ": ", hspncPp, " Hispanic residents"),
              popup = hispanic_popup,
              highlightOptions = highlightOptions(opacity = 1.0, fillColor = "black"),
              group = "Hispanic Population") %>%
  addTiles(group = "OSM") %>%
  addProviderTiles("CartoDB.Positron", group = "Positron") %>%
  addProviderTiles("CartoDB.DarkMatter", group = "Dark") %>%
  addLegend("bottomright",
          colors = RColorBrewer::brewer.pal(5, "YlOrRd"),
          labels = paste0("up to ", round(quantile(acsCovidDat$covid_rate, 
                          c(0.2, 0.4, 0.6, 0.8, 1.0), na.rm = TRUE), 1)),
          title = "COVID Rate<br>per 10,000",
          group = "COVID Rate",
          opacity = 0.8)%>%
  addLayersControl(baseGroups = c("Positron", "OSM", "Dark"),
                   overlayGroups = c("COVID Rate", "Elderly Population", "Hispanic Population"),
                   options = layersControlOptions(collapsed = FALSE)) %>%
  hideGroup(c("Elderly Population", "Hispanic Population"))

htmlMap
htmlwidgets::saveWidget(htmlMap, 'lab3_interactive.html')

Extra features added: multiple basemap providers (Positron, OSM, Dark Matter), polygon highlight on hover via highlightOptions, and layer toggle control with addLayersControl.