Task 1

Plot at least two high-quality static maps with one using the COVID-19 data and one using a related factor. You can use either plot method for sf or ggplot method.

sf_zip <- st_read(dsn = here("Module 9/Data/ZIP_CODE_040114.shp")) %>% 
  janitor::clean_names()

sf_zip_2831 <- st_transform(sf_zip, 2831)

d_covid <- read_csv(file = here("Module 10/Data/tests-by-zcta_2021_04_23.csv")) %>% 
  janitor::clean_names()


sf_zip_2831$zipcode <- as.double(sf_zip_2831$zipcode)

j_zip_covid <- inner_join(x = sf_zip_2831, y =  d_covid,
                          by = c("zipcode" = "modified_zcta")) %>% 
  select(-url, -shape_area, -shape_len, -bldgzip) %>%
  mutate(zip_area = as.numeric(st_area(geometry)/1e6)) # squared kilometers


sf_cov_4326 <- st_transform(x = j_zip_covid, crs = 4326) 

## Getting the basemap

mapBound <- sf_cov_4326 %>% 
  st_bbox() %>%
  st_as_sfc() %>%
  st_buffer(0.02) %>%
  st_bbox() %>%
  as.numeric()

nyc_basemap <- ggmap::get_stamenmap(bbox = mapBound, 
                                    zoom = 13, messaging = FALSE,
                                    maptype = 'terrain-background')


ggmap::ggmap(nyc_basemap) +
  geom_sf(data = sf_cov_4326, aes(fill = population),
          inherit.aes = F) +
  labs(title = "NYC Population", subtitle = "Aggregated by zipcode",
       fill = "Population") +
  xlab("Longitude") +
  ylab("Latitude") +
  ggspatial::annotation_north_arrow(pad_y = unit(8.5, "cm"), 
                                    style = ggspatial::north_arrow_fancy_orienteering()) +
  ggspatial::annotation_scale(pad_x = unit(7, "cm")) +
  theme_light() +
  scale_fill_viridis_c(alpha = 0.6, option = "mako", direction = -1) 


ggmap::ggmap(nyc_basemap) +
  geom_sf(data = sf_cov_4326, aes(fill = covid_death_rate),
          inherit.aes = F) +
  labs(title = "NYC COVID-19 Death Rate", subtitle = "Aggregated by zipcode",
       xlab = "Longitude", ylab = "Latitude", fill = "Population") +
  ggspatial::annotation_north_arrow(pad_y = unit(8.5, "cm"), 
                                    style = ggspatial::north_arrow_fancy_orienteering()) +
  ggspatial::annotation_scale(pad_x = unit(7, "cm")) +
  theme_light() +
  scale_fill_viridis_c(alpha = 0.6, option = "inferno", direction = -1)


Task 2

Use ggplot2 and other ggplot-compatible packages to create a multi-map figure illustrating the possible relationship between COVID-19 confirmed cases or rate and another factor (e.g., the number of nursing homes, number of food stores, neighborhood racial composition, elderly population, etc.). The maps should be put side by side on one single page. Add graticule to at least one of those maps and label some of the feature on the map where applicable and appropriate.

sf_pop_zip <- st_read(here("Module 11/1. data/acsPopByZip.shp")) %>% 
  st_transform(crs = 4326)

sf_cov_pop <- st_join(x = sf_cov_4326, y = sf_pop_zip) %>% 
  mutate(eld_pc = eldrlyP/totPop) 

toner_basemap <- ggmap::get_stamenmap(bbox = mapBound, 
                                    zoom = 13, messaging = FALSE,
                                    maptype = 'toner-lite')

p_cov <- ggmap::ggmap(toner_basemap) +
  geom_sf(data = sf_cov_4326,aes(fill = covid_death_rate), colour = "gray50") +
  labs(title = "COVID-19 Death Rate", subtitle = "Aggregated by Zipcode",
       fill = "Death Rate") +
  xlab("Longitude") +
  ylab("Latitude") +
  # ggspatial::annotation_north_arrow(pad_x = unit(16, "cm"), 
  #                                   pad_y = unit(1, "cm"), 
  #                                   style = ggspatial::north_arrow_fancy_orienteering()) +
  # ggspatial::annotation_scale(pad_x = unit(14, "cm")) +
  scale_fill_viridis_c(alpha = 0.4, direction = -1, option = "inferno") +
  theme_light()


p_pop <- ggmap::ggmap(toner_basemap) +
  geom_sf(data = sf_cov_pop, aes(fill = eld_pc), colour = "gray50") +
  labs(title = "Proportion of Elderly Population",
       subtitle = "Aggregated by Zipcode",
       xlab = "Longitude", ylab = "Latitude",
       fill = "Elder Population (%)") +
  ggspatial::annotation_north_arrow(pad_x = unit(6.2, "cm"),
                                    pad_y = unit(0.5, "cm"), 
                                    width = unit(1, "cm"),
                                    height = unit(1, "cm"),
                                    style = ggspatial::north_arrow_fancy_orienteering()) +
  ggspatial::annotation_scale(pad_x = unit(5, "cm")) +
  scale_fill_viridis_c(alpha = 0.2, direction = -1, option = "viridis") +
  theme_light()


Task 3

Create a web-based interactive map for COIVD-19 data using tmap or leaflet package and save it as a HTML file.

pal_cov <- colorQuantile(palette = "YlGnBu", n = 9, alpha = 0.5,
                         domain = NULL)

p_popup <- paste0("<strong>Borough:</strong>",
                  sf_cov_4326$borough_group,
                  " <br/>","<strong>Neighboorhood:</strong>",
                  sf_cov_4326$neighborhood_name,
                  " <br/>",
                  "<strong> Population: </strong>",
                  sf_cov_4326$population,
                  " <br/>",
                  "<strong>COVID Death Rate: </strong>", 
                  sf_cov_4326$covid_death_rate %>% 
                    round(3)%>%
                    format(nsmall = 3), 
                  " /sqkm",
                  sep="")

imap_cov <- leaflet(data = sf_cov_4326) %>% 
  addPolygons(stroke = F, 
              opacity = 0.3, 
              fillColor =  ~pal_cov(x = sf_cov_4326$covid_death_rate),
              fillOpacity = 0.6,
              smoothFactor = 0.6, 
              popup = p_popup) %>% 
  addTiles()

# htmlwidgets::saveWidget(widget = imap_cov, 
#                         file = here("Module 11/html_covid_map.html"))
imap_cov