I ended up reading the outputs from the prior lab [see lab8] into a new file to be moved into the data directory for this visualization as well, as it already grouped the tract geometry plus demographic attributes.

covidPopZipNYC <- st_read(
  "data/covid_acs.gpkg",
  layer = "final_zip_covid_acs",
  quiet = TRUE
)

Task1:

For this first task, I used the final ZIP-level sf object created in the previous lab, since that dataset already contained the ZIP polygon geometry together with the COVID-19 variables and the ACS demographic variables aggregated to the ZIP code level. I then produced two separate static choropleth maps. The first map shows COVID-19 case rate by ZIP code, and the second shows the percentage of elderly residents by ZIP code as a related demographic factor. As for the second variable, the elderlyPct field was created by dividing by the total ZIP-level population. . Because both variables are continuous and somewhat unevenly distributed, I classified them into quantile-based categories before mapping them, following the cholropleth workflow shown in 3.3 (using the Brewer color palette).

# Use the final ZIP-level sf object from the previous lab
zip_maps <- covidPopZipNYC %>%
  mutate(
    elderlyPct = if_else(totPop > 0, elderlyPop / totPop * 100, NA_real_)
  )
# Extent box so that both maps have the same shape and size bb <- st_bbox(zip_maps)
bb <- st_bbox(zip_maps)


# quantile breaks for COVID case rate
covid_breaks <- classInt::classIntervals(
  c(min(zip_maps$COVID_CASE_RATE, na.rm = TRUE) - 0.00001,
    zip_maps$COVID_CASE_RATE),
  n = 7,
  style = "quantile"
)


# quantile breaks for elderly percentage
elderly_breaks <- classInt::classIntervals(
  c(min(zip_maps$elderlyPct, na.rm = TRUE) - 0.00001,
    zip_maps$elderlyPct),
  n = 7,
  style = "quantile"
)

# continuous variables --> choropleth categories
zip_maps <- zip_maps %>%
  mutate(
    covid_case_cat = cut(COVID_CASE_RATE, covid_breaks$brks, dig.lab = 4),
    elderly_cat = cut(elderlyPct, elderly_breaks$brks, dig.lab = 4)
  )

# Map 1: COVID-19 case rate by ZIP code
g_covid <- ggplot(zip_maps) +
  geom_sf(aes(fill = covid_case_cat), color = "white", linewidth = 0.15) +
  scale_fill_brewer(palette = "OrRd", name = "Case rate") +
  coord_sf(
    xlim = c(bb["xmin"], bb["xmax"]),
    ylim = c(bb["ymin"], bb["ymax"]),
    default_crs = st_crs(4326),
    expand = FALSE
  ) +
  labs(
    title = "COVID-19 Case Rate by ZIP Code",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "grey80", linewidth = 0.3),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
    plot.title = element_text(face = "bold", hjust = 0.5)
  )

# Map 2: elderly population percentage by ZIP code
g_elderly <- ggplot(zip_maps) +
  geom_sf(aes(fill = elderly_cat), color = "white", linewidth = 0.15) +
  scale_fill_brewer(palette = "Blues", name = "Percent age 65+") +
  coord_sf(
    xlim = c(bb["xmin"], bb["xmax"]),
    ylim = c(bb["ymin"], bb["ymax"]),
    default_crs = st_crs(4326),
    expand = FALSE
  ) +
  labs(
    title = "Elderly Population Percentage by ZIP Code",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "grey80", linewidth = 0.3),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
    plot.title = element_text(face = "bold", hjust = 0.5)
  )

#showing
par(mfrow = c(1, 2))

print(g_covid)

print(g_elderly)

par(mfrow = c(1, 1))

# had to explicitly instantiate because it wasn't showing up as a ggplot object otherwise for some reason...

Task 2 – Combined map of potential relatoinship between COVID-19 rate and elderly population share

Same as before; I used the same ZIP-level sf object and derived variables from the previous step to create a single comparative figure using plot_layout, which should make two maps the same size regardless of legends or other elements.

library(patchwork)
# Put the two maps side by side on one page
g_covid <- g_covid +
  theme(legend.position = "bottom") +
  labs(title = "COVID-19 Case Rate")

g_elderly <- g_elderly +
  theme(legend.position = "bottom") +
  labs(title = "Elderly Population Share")

g_covid + g_elderly + plot_layout(ncol = 2, nrow = 1)

Task 3 – Creating an Interactive Web map with Leaflet

To wrap up, First, I packaged a new share of % elderly residents grouped by ZIP code (leaflet_data) in WGS84 (just to make things a bit cleaner). After this, I added layer controls [*] which let the user switch basemaps, toggle overlay layers (light, dark, and OSM), and open one main layer visible with the other two hidden. There are three main layers: [1], The COVID Case rate; [2] The ZIP-level positivity rate; and [3], The Elderly (65+) share of population per ZIP. The popups and hover labels aren’t R-native, so I had to use HTML-formatted strings, which leaflet reads in the browser. lapply() is used to just tell Leaflet to interpret these strings as HTML rather than literal text. At this moment, I’ve just been using the add_legend fnuctionality, which unfortunately does not allow for ramp agglomoration, so I am rather wastefully trying to show 3 ramps at once (while, preferentially, I would just use one custom HTML with all three colors shown per % to save space).

library(leaflet) # using leaflet!
library(htmlwidgets) # to save to html file

#packaging elderly % per zip 
leaflet_data <- covidPopZipNYC %>%
  mutate(
    elderlyPct = if_else(totPop > 0, elderlyPop / totPop * 100, NA_real_)
  ) %>%
  st_transform(4326)

names(leaflet_data) # out 
##  [1] "ZIPCODE"           "BLDGZIP"           "PO_NAME"          
##  [4] "POPULATION"        "AREA"              "STATE"            
##  [7] "COUNTY"            "ST_FIPS"           "CTY_FIPS"         
## [10] "URL"               "SHAPE_AREA"        "SHAPE_LEN"        
## [13] "NEIGHBORHOOD_NAME" "BOROUGH_GROUP"     "label"            
## [16] "lat"               "lon"               "COVID_CASE_COUNT" 
## [19] "COVID_CASE_RATE"   "POP_DENOMINATOR"   "COVID_DEATH_COUNT"
## [22] "COVID_DEATH_RATE"  "PERCENT_POSITIVE"  "TOTAL_COVID_TESTS"
## [25] "FoodStoreNum"      "NursingHomeNum"    "totPop"           
## [28] "elderlyPop"        "malePop"           "femalePop"        
## [31] "whitePop"          "blackPop"          "asianPop"         
## [34] "hispanicPop"       "adultPop"          "citizenAdult"     
## [37] "geom"              "elderlyPct"
# Color pallettes (one func per variable)
#for some reason Ive been having some issues with the leaflet namespace so I just have to manually instantiate...
pal_case <- colorQuantile(
  palette = "YlOrRd",
  domain = leaflet_data$COVID_CASE_RATE,
  n = 7,
  na.color = "transparent"
)

pal_pos <- colorQuantile(
  palette = "PuRd",
  domain = leaflet_data$PERCENT_POSITIVE,
  n = 7,
  na.color = "transparent"
)

pal_elderly <- colorQuantile(
  palette = "Blues",
  domain = leaflet_data$elderlyPct,
  n = 7,
  na.color = "transparent"
)

# Hover labels
label_case <- paste0(
  "ZIP: ", leaflet_data$ZIPCODE,
  "<br>Case rate: ", round(leaflet_data$COVID_CASE_RATE, 1)
)

label_pos <- paste0(
  "ZIP: ", leaflet_data$ZIPCODE,
  "<br>Percent positive: ", round(leaflet_data$PERCENT_POSITIVE, 2), "%"
)

label_elderly <- paste0(
  "ZIP: ", leaflet_data$ZIPCODE,
  "<br>Elderly share: ", round(leaflet_data$elderlyPct, 2), "%"
)

#---------------
# Pop-ups (when user clicks on ZIP polygon)
popup_case <- paste0(
  "<strong>ZIP Code:</strong> ", leaflet_data$ZIPCODE, "<br>",
  "<strong>Neighborhood:</strong> ", leaflet_data$PO_NAME, "<br>",
  "<strong>COVID Case Count:</strong> ", leaflet_data$COVID_CASE_COUNT, "<br>",
  "<strong>COVID Case Rate:</strong> ", round(leaflet_data$COVID_CASE_RATE, 1)
)

popup_pos <- paste0(
  "<strong>ZIP Code:</strong> ", leaflet_data$ZIPCODE, "<br>",
  "<strong>Total COVID Tests:</strong> ", leaflet_data$TOTAL_COVID_TESTS, "<br>",
  "<strong>Percent Positive:</strong> ", round(leaflet_data$PERCENT_POSITIVE, 2), "%"
)

popup_elderly <- paste0(
  "<strong>ZIP Code:</strong> ", leaflet_data$ZIPCODE, "<br>",
  "<strong>Total Population:</strong> ", leaflet_data$totPop, "<br>",
  "<strong>Elderly Population:</strong> ", leaflet_data$elderlyPop, "<br>",
  "<strong>Elderly Share:</strong> ", round(leaflet_data$elderlyPct, 2), "%"
)

#-------------
# Hover style options lol:
poly_highlight <- highlightOptions(
  weight = 2,
  color = "black",
  fillOpacity = 0.9,
  bringToFront = TRUE
)

poly_label <- labelOptions(
  style = list(
    "font-weight" = "normal",
    "padding" = "3px 8px"
  ),
  textsize = "12px",
  direction = "auto"
)

#-------
#Start map:
covid_leaflet <- leaflet(leaflet_data) %>%
  addProviderTiles("CartoDB.Positron", group = "Light") %>% # basemap
  addProviderTiles("CartoDB.DarkMatter", group = "Dark") %>%
  addProviderTiles("OpenStreetMap", group = "OSM")

# First overlay (COVID rate):
covid_leaflet <- covid_leaflet %>%
  addPolygons(
    fillColor = ~pal_case(COVID_CASE_RATE),
    color = "white",
    weight = 0.5,
    fillOpacity = 0.8,
    group = "COVID Case Rate",
    label = lapply(label_case, htmltools::HTML),
    popup = lapply(popup_case, htmltools::HTML),
    highlightOptions = poly_highlight,
    labelOptions = poly_label
  ) %>%
  addLegend(
    position = "bottomright",
    pal = pal_case,
    values = ~COVID_CASE_RATE,
    title = "COVID Case Rate",
    opacity = 0.8
  )

#----------------
# % Positive overlay
covid_leaflet <- covid_leaflet %>%
  addPolygons(
    fillColor = ~pal_pos(PERCENT_POSITIVE),
    color = "white",
    weight = 0.5,
    fillOpacity = 0.8,
    group = "Percent Positive",
    label = lapply(label_pos, htmltools::HTML),
    popup = lapply(popup_pos, htmltools::HTML),
    highlightOptions = poly_highlight,
    labelOptions = poly_label
  ) %>%
  addLegend(
    position = "topleft",
    pal = pal_pos,
    values = ~PERCENT_POSITIVE,
    title = "Percent Positive",
    opacity = 0.8
  )

# Elderly Share Overlay
covid_leaflet <- covid_leaflet %>%
  addPolygons(
    fillColor = ~pal_elderly(elderlyPct),
    color = "white",
    weight = 0.5,
    fillOpacity = 0.8,
    group = "Elderly Share",
    label = lapply(label_elderly, htmltools::HTML),
    popup = lapply(popup_elderly, htmltools::HTML),
    highlightOptions = poly_highlight,
    labelOptions = poly_label
  ) %>%
  addLegend(
    position = "bottomleft",
    pal = pal_elderly,
    values = ~elderlyPct,
    title = "Elderly Share",
    opacity = 0.8
  )

# Layer Controls [*]
covid_leaflet <- covid_leaflet %>%
  addLayersControl(
    baseGroups = c("Light", "Dark", "OSM"),
    overlayGroups = c("COVID Case Rate", "% Positive", "Elderly Share"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  hideGroup(c("Percent Positive", "Elderly Share")) %>%
  addControl(
    html = "<strong>Map Layers</strong><br>Select a basemap and toggle overlays here.", #man
    position = "topright"
  )

covid_leaflet
# To make saveable as individual html:

final_map <- leaflet(leaflet_data) %>%
  addTiles() %>%
  addPolygons(
    fillColor = ~pal_pos(PERCENT_POSITIVE),
    fillOpacity = 0.7,
    color = "black",
    weight = 1
  )

saveWidget(final_map, "final_leaflet_map.html", selfcontained = TRUE)