GGPLOT Maps

For this section I wanted to compare Covid cases with food store locations based on the work from last lab.

# For this lab I created a RData file from the previous lab's work
# file is called zip_fstores_hfacilities_demographics_sf
# to preserve this original data I create a duplicate of the same name
# lab2_data_sf protects the rdata file from any issues that might arise over the analysis
load("Lab2Data.rdata")
# I start with an examination of the variables to reaquiant myself with the data
# list_of_columns <- names(zip_covid_fstore_hfacility_demographics_sf)
# print(list_of_columns)
# create the duplicate sf that will be the working object for the rest of the lab
lab2_data_sf <- zip_covid_fstore_hfacility_demographics_sf

# I need to clean the NA values here so I replace them with zero,
# this is for areas like central park with no residents
# I also prepare some new variables to get percentages of elderly by zip code
# my imported sf did not have rates, so I generate some just in case I need them later
lab2_data_sf <- lab2_data_sf %>%
  mutate(across(.cols = where(is.numeric), ~replace(., is.na(.), 0)),
         total_population_65_years_and_over = total_population_65_years_and_over_male + total_population_65_years_and_over_female,
         percent_total_population_65_years_and_over = if_else(total_population == 0, 0, (total_population_65_years_and_over / total_population) * 100),
         covid_death_rate_1k = if_else(POPULATION == 0, 0, (covid_death_count / POPULATION) * 1000),
         covid_case_rate_1k = if_else(POPULATION == 0, 0, (covid_case_count / POPULATION) * 1000),
         estimate_race_total_population_nonwhite = estimate_race_total_population_one_race_american_indian_and_alaska_native + estimate_race_total_population_one_race_asian + estimate_race_total_population_one_race_black_or_african_american + estimate_race_total_population_one_race_native_hawaiian_and_other_pacific_islander + estimate_hispanic_or_latino_and_race_total_population)

# now to make maps

# # Map 1: death rates by zip code choropleth
# ggplot(data = lab2_data_sf) +
#   geom_sf(aes(fill = covid_death_rate_1k), color = NA) +
#   scale_fill_viridis_c(option = "C", direction = -1, guide = guide_legend(title = "Covid Death Rate")) +
#   labs(title = "Covid Death Rates per 1000 by Zip Code") +
#   theme_minimal() +
#   coord_sf(datum = NA)
# 
# # Map 2: mapping covid case rates by zip code choropleth
# ggplot(data = lab2_data_sf) +
#   geom_sf(aes(fill = covid_death_rate_1k), color = NA) +
#   scale_fill_viridis_c(option = "C", direction = -1, guide = guide_legend(title = "Covid Case Rate")) +
#   labs(title = "Covid Case Rates per 1000 by Zip Code") +
#   theme_minimal() +
#   coord_sf(datum = NA)
# 
# # Map 3: mapping elderly population density by zip code choropleth
# ggplot(data = lab2_data_sf) +
#   geom_sf(aes(fill = percent_total_population_65_years_and_over), color = NA) +
#   scale_fill_viridis_c(option = "C", direction = -1, guide = guide_legend(title = "% Population Over 65")) +
#   labs(title = "Percentage of Population Over 65 by Zip Code") +
#   theme_minimal() +
#   coord_sf(datum = NA)

# Map 4: case counts by zip code choropleth
covid_map <- ggplot(data = lab2_data_sf) +
  geom_sf(aes(fill = covid_case_count), color = NA) +
  scale_fill_viridis_c(option = "C", direction = -1, guide = guide_legend(title = "Covid case Count")) +
  labs(title = "Covid Case Counts by Zip Code") +
  theme_minimal() +
  coord_sf(datum = NA)

# Map 5: mapping food stores by zip code choropleth
food_map <- ggplot(data = lab2_data_sf) +
  geom_sf(aes(fill = store_count), color = NA) +
  scale_fill_viridis_c(option = "C", direction = -1, guide = guide_legend(title = "Food Store Count")) +
  labs(title = "Food Store Density by Zip Code") +
  theme_minimal() +
  coord_sf(datum = NA)

ggarrange(covid_map, food_map, ncol = 2, nrow = 1)

Leaflet Map

In this section I wanted to see if zip codes with a higher population of non-white ethnicities faced greater covid case rates than those with higher populations of white people.

# Leaflet Interactive Map
pal_fun <- colorQuantile("YlOrRd", NULL, n = 5)

pal_eth <- colorQuantile("YlGnBu", 
                             domain = c(min(lab2_data_sf$estimate_race_total_population_nonwhite, na.rm = T), max(lab2_data_sf$estimate_race_total_population_nonwhite, na.rm = T)),
                             alpha = TRUE)

covid_popup <- paste0("<strong>Zip Code: </strong>", 
                      lab2_data_sf$ZIPCODE,
                      " <br/>",
                      "<strong> Place Name: </strong>",
                      lab2_data_sf$PO_NAME,
                      " <br/>",
                      "<strong> Confirmed COVID-19 Cases: </strong>",
                      lab2_data_sf$covid_case_count,
                      sep="")

polyHighlightOption <- leaflet::highlightOptions(opacity = 1.0, fillColor = 'black')
polyLabelOption <- leaflet::labelOptions(opacity = 1.0, textsize = '11px')

htmlMap <- leaflet(lab2_data_sf %>% sf::st_transform(4326)) %>%
  addPolygons(
    stroke = TRUE, # show the polygon boundary
    fillOpacity = 0.8,
    fillColor = ~pal_eth(estimate_race_total_population_nonwhite),
    color = 'grey', # color of the boundary
    weight = 1, # width of the boundary
    highlightOptions = polyHighlightOption,
    label = paste(lab2_data_sf$estimate_race_total_population_nonwhite),
    group = "Nonwhite Population") %>%
  addPolygons(
    stroke = FALSE,
    fillColor = ~pal_fun(covid_case_count),
    fillOpacity = 0.8, smoothFactor = 0.5,
    popup = covid_popup,
    group = "NYC-COVID",
    label = paste0(lab2_data_sf$covid_case_count, " at ", lab2_data_sf$PO_NAME, ", ", lab2_data_sf$ZIPCODE),
    highlightOptions = polyHighlightOption,
    labelOptions = polyLabelOption) %>%
  addTiles(group = "OSM") %>%
  addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
  addLayersControl(baseGroups = c("OSM", "Carto"), 
                   overlayGroups = c("NYC-COVID", "Nonwhite Population"))  

# adding legends based on which layer is selected
htmlMap <- htmlMap %>%
  addLegend(pal = pal_eth, values = ~estimate_race_total_population_nonwhite, opacity = 0.7,
            title = "Non-white Population Density",
            position = "bottomright",
            group = "Nonwhite Population") %>%
  addLegend(pal = pal_fun, values = ~covid_case_count, opacity = 0.7,
            title = "COVID-19 Case Counts",
            position = "bottomright",
            group = "NYC-COVID")

htmlMap
# commented out the save step
# mapshot(htmlMap, "nyc_covid_nonwhite_pop_burden.html")