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