Load and Prepare the Data

# Convert from NAD83 / New York Long Island (ftUS) to WGS84
acsCovidDat <- st_transform(acsCovidDat, 4326)

# Define color palettes
pal_covid <- brewer.pal(5, "YlOrRd")
pal_elderly <- brewer.pal(5, "YlGnBu")

STATIC MAP 1: COVID-19 CASES MAP (ggplot2)

# Create the COVID-19 map with labels
gg_covid <- ggplot(acsCovidDat) + 
  geom_sf(aes(fill=Positiv), color="black", size=0.2) +
  scale_fill_viridis_c(option = "plasma", name="COVID Cases") +
  coord_sf() +
  labs(title="COVID-19 Confirmed Cases by ZIP Code",
       x="Longitude", y="Latitude") +
  theme_minimal() +
  geom_text(data = acsCovidDat_centroids, 
            aes(x = x, y = y, label = ZIPCODE), 
            size = 1, color = "gray")

# Set legend position
gg_covid <- gg_covid + theme(legend.position = "bottom")

# Save the static map
ggsave("covid_static_map.png", gg_covid, width=8, height=6)

# Display the image in the rendered document
knitr::include_graphics("covid_static_map.png")

STATIC MAP 2: ELDERLY POPULATION MAP (ggplot2)

# Create the Elderly Population map with labels
gg_elderly <- ggplot(acsCovidDat) + 
  geom_sf(aes(fill=eldrlyP), color="black", size=0.2) +
  scale_fill_gradientn(colors=pal_elderly, name="Elderly Population (%)") +
  coord_sf() +
  labs(title="Elderly Population by ZIP Code",
       x="Longitude", y="Latitude") +
  theme_minimal() +
  geom_text(data = acsCovidDat_centroids, 
            aes(x = x, y = y, label = ZIPCODE), 
            size = 1, color = "black")  # Smaller, less distracting labels

# Set legend position
gg_elderly <- gg_elderly + theme(legend.position = "bottom")

# Save the static map
ggsave("elderly_population_map.png", gg_elderly, width=8, height=6)

# Display the image in the rendered document
knitr::include_graphics("elderly_population_map.png")

MULTI-MAP FIGURE: COVID CASES & ELDERLY POPULATION (ggpubr)

# -------------------------------------------------------------------------
# MULTI-MAP FIGURE: COVID CASES & ELDERLY POPULATION (ggpubr)
# -------------------------------------------------------------------------

# Combine maps with equal heights
multi_map <- ggarrange(gg_covid, gg_elderly, 
                       ncol=2, nrow=1, heights = c(1,1))

# Save multi-map figure
ggsave("multi_map_figure.png", multi_map, width=12, height=6)

# Display the image in the rendered document
knitr::include_graphics("multi_map_figure.png")

INTERACTIVE WEB MAP (Leaflet)

# -------------------------------------------------------------------------
# INTERACTIVE WEB MAP (Leaflet)
# -------------------------------------------------------------------------
# Define Leaflet color palettes
pal_fun <- colorQuantile("YlOrRd", NULL, n = 5)
pal_elderly <- colorQuantile("YlGnBu", domain = range(acsCovidDat$eldrlyP, na.rm=TRUE))

# Create popups for interactivity
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)

# Build interactive map
htmlMap <- leaflet(acsCovidDat) %>%
  addTiles(group = "OpenStreetMap") %>%
  addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
  addPolygons(
    fillColor = ~pal_elderly(eldrlyP),
    fillOpacity = 0.8, color = "grey", weight = 1,
    popup = covid_popup,
    group = "Elderly Population") %>%
  addPolygons(
    fillColor = ~pal_fun(Positiv),
    fillOpacity = 0.8, stroke = FALSE,
    popup = covid_popup,
    group = "COVID-19 Cases") %>%
  addLayersControl(baseGroups = c("OpenStreetMap", "Carto"),
                   overlayGroups = c("COVID-19 Cases", "Elderly Population"))

# Save the interactive map as an HTML file
htmlwidgets::saveWidget(htmlMap, "covid_interactive_map.html", selfcontained = FALSE)

htmlMap