task 1: Plot at least two high-quality static maps at the zip code level.

# Map 1:
zip_clean <- zipcodes_with_data %>%
  filter(!is.na(overall_positivity_rate)) #Filter out NA values

map1 <- ggplot(zip_clean) +
  geom_sf(aes(fill = overall_positivity_rate), color = "white", size = 0.2) +
  scale_fill_viridis_c(option = "plasma", 
                       name = "Positivity Rate (%)") +
  theme_minimal() +
  labs(title = "Positivity Rate in NYC (4/12-4/23/20)",
  )

map1

# Map 2: 
map2 <- ggplot(zip_clean) +
  geom_sf(aes(fill = residential_health), color = "white", size = 0.2) +
  scale_fill_viridis_c(option = "plasma", 
                       name = "Number of Residential Health Care Facilities") +
  theme_minimal() +
  labs(title = "Residential Healthcare Facilities in NYC",
  )

map2

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.

zipcodes_with_data_and_acs <- base::merge(zipcodes_with_data, aggregated, by = "MODZCTA")

# Map 1:
zipcodes_clean <- zipcodes_with_data %>%
  filter(!is.na(overall_positivity_rate)) #Filter out NA values

map1 <- ggplot(zipcodes_clean) +
  geom_sf(aes(fill = overall_positivity_rate), color = "white", size = 0.2) +
  scale_fill_viridis_c(option = "plasma", 
                       name = "Positivity Rate (%)") +
  theme_minimal() +
  labs(title = "Positivity Rate in NYC (4/12-4/23/20)")

# Map 2: 
mapold <- ggplot(zipcodes_with_data_and_acs) +
  geom_sf(aes(fill = medianage), color = "white", size = 0.2) +
  scale_fill_viridis_c(option = "plasma", 
                       name = "Median Age") +
  theme_minimal() +
  labs(title = "Median Age")

map1 + mapold

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

map_leaflet <- leaflet(zipcodes_clean) %>%
  addTiles() 

map_leaflet <- map_leaflet %>%
  addPolygons(
    fillColor = ~colorNumeric("Reds", overall_positivity_rate)(overall_positivity_rate),
    fillOpacity = 0.7,
    color = "white",
    weight = 1,
    group = "COVID-19 Positive Rate (%)",
    label = ~paste("ZIP:", MODZCTA, "Positivity Rate:", overall_positivity_rate),
    popup = ~paste("COVID-19 Positive Rate (%):", overall_positivity_rate)
  )

map_leaflet <- map_leaflet %>%
  addPolygons(
    fillColor = ~colorNumeric("Blues", residential_health)(residential_health),
    fillOpacity = 0.7,
    color = "white",
    weight = 1,
    group = "# of Residential Health Care Facilities",
    label = ~paste("Residential Health Care Facilities:", residential_health),
    popup = ~paste("Residential Health Care Facilities:", residential_health)
  )

coords <- st_coordinates(facilities_sf)

facilities_clean <- facilities_sf %>%
  mutate(lon = coords[,1],
         lat = coords[,2]) %>%
  filter(lon >= -74.3 & lon <= -73.6) %>%
  filter(lat >= 40.4 & lat <= 41.0) %>%
  select(-lon, -lat) 

map_leaflet <- map_leaflet %>%
  addCircleMarkers(
    data = facilities_clean,
    radius = 3,
    color = "lightgreen",
    fillColor = "lightgreen",
    group = "Locations of Health Care Facilities",
    label = ~Facility.Name,
    popup = ~paste("Name:", Facility.Name, "Description", Description)
  ) 

map_leaflet <- map_leaflet %>%
  addLayersControl(
    overlayGroups = c("COVID-19 Positive Rate (%)", "# of Residential Health Care Facilities", "Locations of Health Care Facilities"),
    options = layersControlOptions(collapsed = FALSE)
  )

map_leaflet <- map_leaflet %>%
  addLegend(
    "bottomleft",
    pal = colorNumeric("Reds", zip_clean$overall_positivity_rate),
    values = zip_clean$overall_positivity_rate,
    title = "COVID-19 Positive Rate",
    group = "COVID-19 Positive Rate"
  ) 

map_leaflet