Introduction

This lab focuses on visualizing spatial data at the ZIP code level in New York City. The goal of the assignment is to create both static and interactive maps that clearly communicate patterns in COVID-19 cases and a related demographic factor.

For this analysis, I used a GeoPackage dataset that contains ZIP code boundaries along with COVID-19 case counts and demographic variables. I chose to compare COVID-19 confirmed cases with the elderly population because there is a clear relationship between vulnerability and age. This allows for a meaningful spatial comparison across neighborhoods.

The assignment is divided into three main tasks. First, I create two high-quality static choropleth maps. Second, I combine those maps into a single comparison figure. Finally, I build an interactive map using leaflet with multiple layers and user interaction features.

Setup

library(sf)
library(ggplot2)
library(classInt)
library(dplyr)
library(RColorBrewer)
library(ggpubr)
library(grid)
library(leaflet)
library(htmlwidgets)

Task 1: Static Choropleth Maps

COVID-19 Confirmed Cases by ZIP Code (NYC)

breaks_qt <- classIntervals(acs_data$Positiv, n = 5, style = "quantile")
breaks_clean <- round(breaks_qt$brks, 0)

labels_clean <- paste0(
  head(breaks_clean, -1), " – ", tail(breaks_clean, -1)
)

acs_data$Positiv_cat <- cut(
  acs_data$Positiv,
  breaks = breaks_clean,
  labels = labels_clean,
  include.lowest = TRUE
)

covid_map <- ggplot(data = acs_data) +
  geom_sf(aes(fill = Positiv_cat), color = "grey40", linewidth = 0.1) +
  scale_fill_brewer(
    palette = "Blues",
    name = "COVID Cases",
    na.value = "lightgrey",
    na.translate = FALSE
  ) +
  labs(
    title = "COVID-19 Confirmed Cases",
    x = "Longitude",
    y = "Latitude",
    caption = "Data Source: NYC Open Data, US Census"
  ) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray80"),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    plot.title = element_text(face = "bold", size = 16),
    plot.caption = element_text(size = 9, color = "gray40")
  )

covid_map

Elderly Population by ZIP Code (NYC)

breaks_elderly <- classIntervals(
  acs_data$eldrlyP[!is.na(acs_data$eldrlyP)],
  n = 5,
  style = "quantile"
)
breaks_elderly_clean <- round(breaks_elderly$brks, 0)

labels_elderly <- paste0(
  head(breaks_elderly_clean, -1), " – ", tail(breaks_elderly_clean, -1)
)

acs_data$eldrly_cat <- cut(
  acs_data$eldrlyP,
  breaks = breaks_elderly_clean,
  labels = labels_elderly,
  include.lowest = TRUE
)

elderly_map <- ggplot(data = acs_data) +
  geom_sf(aes(fill = eldrly_cat), color = "grey40", linewidth = 0.1) +
  scale_fill_brewer(
    palette = "Reds",
    name = "Elderly Population",
    na.value = "lightgrey",
    na.translate = FALSE
  ) +
  labs(
    title = "Elderly Population",
    x = "Longitude",
    y = "Latitude",
    caption = "Data Source: NYC Open Data, US Census"
  ) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray80"),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    plot.title = element_text(face = "bold", size = 16),
    plot.caption = element_text(size = 9, color = "gray40")
  )

elderly_map

Task 2: Side-by-Side Comparison

combined_map <- ggarrange(
  covid_map,
  elderly_map,
  ncol = 2,
  nrow = 1,
  align = "hv"
)

combined_map <- annotate_figure(
  combined_map,
  top = text_grob(
    "Comparison of COVID-19 Cases and Elderly Population by ZIP Code (NYC)",
    face = "bold",
    size = 18
  )
)

combined_map

Task 3: Interactive Map

acs_leaflet <- st_transform(acs_data, 4326)

pal_covid <- colorQuantile("YlOrRd", domain = acs_leaflet$Positiv, n = 5)
pal_elderly <- colorQuantile("Blues", domain = acs_leaflet$eldrlyP, n = 5)

popup_text <- paste0(
  "<strong>ZIP Code: </strong>", acs_leaflet$ZIPCODE, "<br/>",
  "<strong>COVID Cases: </strong>", acs_leaflet$Positiv, "<br/>",
  "<strong>Elderly Population: </strong>", acs_leaflet$eldrlyP
)
map <- leaflet(acs_leaflet) %>%
  
  addProviderTiles("CartoDB.Positron", group = "Light Map") %>%
  addProviderTiles("CartoDB.DarkMatter", group = "Dark Map") %>%
  
  addPolygons(
    fillColor = ~pal_covid(Positiv),
    fillOpacity = 0.7,
    color = "white",
    weight = 1,
    label = ~paste0("ZIP: ", ZIPCODE),
    popup = popup_text,
    group = "COVID Cases"
  ) %>%
  
  addPolygons(
    fillColor = ~ifelse(
      is.na(eldrlyP),
      "#f0f0f0",
      pal_elderly(eldrlyP)
    ),
    fillOpacity = 0.7,
    color = "white",
    weight = 1,
    label = ~paste0("ZIP: ", ZIPCODE),
    popup = popup_text,
    group = "Elderly Population"
  ) %>%
  
  addLegend(
    pal = pal_covid,
    values = acs_leaflet$Positiv,
    title = "COVID Cases",
    position = "bottomright",
    group = "COVID Cases"
  ) %>%
  
  addLegend(
    pal = pal_elderly,
    values = acs_leaflet$eldrlyP[!is.na(acs_leaflet$eldrlyP)],
    title = "Elderly Population",
    position = "bottomleft",
    group = "Elderly Population"
  ) %>%
  
  addLayersControl(
    baseGroups = c("Light Map", "Dark Map"),
    overlayGroups = c("COVID Cases", "Elderly Population"),
    options = layersControlOptions(collapsed = FALSE)
  )

map
saveWidget(map, "interactive_map.html", selfcontained = TRUE)

Conclusion

In this lab, I created both static and interactive maps to explore spatial patterns of COVID-19 cases and elderly population across New York City ZIP codes. The static maps provide a clear visual comparison, while the interactive map allows for deeper exploration through layers, popups, and user interaction. Together, these approaches highlight how different mapping techniques can be used to better understand geographic patterns in real-world data.