## Loading required package: sf
## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: ggplot2
## Loading required package: ggspatial
## Warning: package 'ggspatial' was built under R version 4.5.3
## Loading required package: mapview
## Loading required package: leaflet
## Loading required package: htmlwidgets
## Loading required package: tmap
## Loading required package: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ lubridate 1.9.5     ✔ tibble    3.3.1
## ✔ purrr     1.2.1     ✔ tidyr     1.3.2
## ✔ readr     2.2.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: cowplot
## 
## 
## Attaching package: 'cowplot'
## 
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
##### Merge NYC zip codes with COVID tests
covid_map <- nyc_zipcodes %>%
  left_join(ZCTA_tests, by = c("ZIPCODE" = "label"))

##### Plotting total COVID tests per zip code
ggplot(covid_map) +
  geom_sf(aes(fill = TOTAL_COVID_TESTS), color = "black", size = 0.2) +
  scale_fill_viridis_c(option = "plasma", name = "Total Tests") +
  labs(title = "COVID-19 Testing by NYC Zip Code",
       subtitle = "Tests by ZCTA",
       x = "Longitude", y = "Latitude") +
  annotation_scale(location = "bl", width_hint = 0.3, pad_y = unit(5, "cm") ) +
  annotation_north_arrow(location = "tl", style = north_arrow_fancy_orienteering()) +
  theme_minimal() +
  theme(legend.position = "right")

##### Removing geometry from ACS data
acs_data <- st_set_geometry(ACS_zipcode, NULL) %>%
  select(county_fips, DP05_0024E) %>%
  distinct()

##### Merging NYC zip codes with ACS Population Data
acs_map <- nyc_zipcodes %>%
  left_join(
    st_set_geometry(ACS_zipcode, NULL) %>% select(county_fips, DP05_0024E),
    by = c("CTY_FIPS" = "county_fips")
  )
## Warning in sf_column %in% names(g): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 14 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
##### Aggregating zip codes to borough polygons
acs_map <- st_join(nyc_zipcodes, ACS_zipcode, join = st_intersects)

acs_joined <- st_join(nyc_zipcodes, ACS_zipcode)

acs_map <- acs_joined %>%
  group_by(ZIPCODE) %>%
  summarise(
    DP05_0024E = sum(DP05_0024E, na.rm = TRUE),
    geometry = st_union(geometry)
  )

##### Plotting map of population under 18 by zip code
ggplot(acs_map) +
  geom_sf(aes(fill = DP05_0024E), color = "black", size = 0.2) +
  scale_fill_viridis_c(option = "magma", name = "Population <18yo") +
  labs(
    title = "Population Under 18 by NYC Zip Code",
    x = "Longitude", y = "Latitude"
  ) +
  annotation_scale(location = "bl", width_hint = 0.3, pad_y = unit(5, "cm")) +
  annotation_north_arrow(location = "tl", style = north_arrow_fancy_orienteering()) +
  theme_minimal()

##### Map 1: Total COVID tests
map_tests <- ggplot(covid_map) +
  geom_sf(aes(fill = TOTAL_COVID_TESTS), color = "black", size = 0.2) +
  scale_fill_viridis_c(option = "plasma", name = "Total Tests", limits = c(0, max(covid_map$TOTAL_COVID_TESTS, na.rm=TRUE))) +
  labs(title = "COVID-19 Tests by Zip Code") +
annotation_scale(
  location = "bl",
  width_hint = 0.3,
  pad_y = unit(6, "cm") ) +
  annotation_north_arrow(location = "tl", style = north_arrow_fancy_orienteering()) +
  theme_minimal() +
  theme(legend.position = "right") +
  scale_x_continuous(n.breaks = 4) +
theme(
  axis.text.x = element_text(angle = 45, hjust = 1)
)

##### Map 2: Population under 18yo
map_population <- ggplot(acs_map) +
  geom_sf(aes(fill = DP05_0024E), color = "black", size = 0.2) +
  scale_fill_viridis_c(option = "magma", name = "Population <18yo", limits = c(0, max(acs_map$DP05_0024E, na.rm=TRUE))) +
  labs(title = "Population Under 18yo by Zip Code") +
 annotation_scale(
  location = "bl",
  width_hint = 0.3,
  pad_y = unit(5, "cm") ) +
  annotation_north_arrow(location = "tl", style = north_arrow_fancy_orienteering()) +
  theme_minimal() +
  theme(legend.position = "right") +
  scale_x_continuous(n.breaks = 4) +
theme(
  axis.text.x = element_text(angle = 45, hjust = 1)
)

##### Multi-map figure creation
plot_grid(map_tests, map_population, ncol = 2, align = "h")

##### Setting up data
covid_map <- nyc_zipcodes %>%
  left_join(ZCTA_tests, by = c("ZIPCODE" = "label"))

nyc_zipcodes <- st_transform(nyc_zipcodes, 4326)
covid_map    <- st_transform(covid_map, 4326)
acs_map      <- st_transform(acs_map, 4326)

##### Aesthetics
pal_covid <- colorNumeric("plasma", covid_map$TOTAL_COVID_TESTS, na.color = "transparent")
pal_acs   <- colorNumeric("magma", acs_map$DP05_0024E, na.color = "transparent")

##### Base map with three layers
m <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  
##### Layer 1: Total COVID tests
addPolygons(
    data = covid_map,
    fillColor = ~pal_covid(TOTAL_COVID_TESTS),
    fillOpacity = 0.7,
    color = "black",
    weight = 1,
    group = "COVID Tests",
    label = ~paste0("ZIP: ", ZIPCODE),
    popup = ~paste0(
      "<b>ZIP:</b> ", ZIPCODE, "<br>",
      "<b>Total Tests:</b> ", TOTAL_COVID_TESTS
    )
  ) %>%
  
##### Layer 2: Population under 18yo
addPolygons(
    data = acs_map,
    fillColor = ~pal_acs(DP05_0024E),
    fillOpacity = 0.7,
    color = "black",
    weight = 1,
    group = "Population <18",
    label = ~paste0("ZIP: ", ZIPCODE),
    popup = ~paste0(
      "<b>ZIP:</b> ", ZIPCODE, "<br>",
      "<b>Population <18:</b> ", DP05_0024E
    )
  ) %>%
  
##### Layer 3: Zip code outlines
  addPolygons(
    data = nyc_zipcodes,
    fill = FALSE,
    color = "blue",
    weight = 1,
    group = "ZIP Boundaries",
    label = ~ZIPCODE
  ) %>%
  
##### Layer control
  addLayersControl(
    overlayGroups = c("COVID Tests", "Population <18", "ZIP Boundaries"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  
##### Legend and Scale Bar
    addLegend(
    pal = pal_covid,
    values = covid_map$TOTAL_COVID_TESTS,
    title = "COVID Tests",
    position = "bottomright"
  ) %>%
    addScaleBar(position = "bottomleft") %>%
  
##### Minimap Overview
    addMiniMap(toggleDisplay = TRUE)

##### Saving the map
htmlwidgets::saveWidget(m, "interactive_map.html", selfcontained = TRUE)
m