Week 9 Assignment

Load packages and previous data:

library(ggplot2)
library(sf)

Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(biscale)
library(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

library(patchwork) 
library(ggspatial)
library(leaflet)
library(htmlwidgets)
library(htmltools)
library(rmapshaper)

# Simplify GeoJSON geometry since file size is too big
original <- st_read("/Users/ayeshanaveed/Downloads/R-Spatial_II_Lab/final_zipcode_data.geojson")

Reading layer `final_zipcode_data' from data source 
  `/Users/ayeshanaveed/Downloads/R-Spatial_II_Lab/final_zipcode_data.geojson' 
  using driver `GeoJSON'
Simple feature collection with 996 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91553
Geodetic CRS:  WGS 84

simplified <- ms_simplify(original, keep=0.02) 
write_sf(simplified, "simplified_zipcodes.geojson")

Warning in CPL_write_ogr(obj, dsn, layer, driver,
as.character(dataset_options), : GDAL Error 6: DeleteLayer() not supported by
this dataset.

final_data <- st_read("simplified_zipcodes.geojson")

Reading layer `simplified_zipcodes' from data source 
  `/Users/ayeshanaveed/Downloads/simplified_zipcodes.geojson' 
  using driver `GeoJSON'
Simple feature collection with 5964 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.25545 ymin: 40.49613 xmax: -73.70001 ymax: 40.9154
Geodetic CRS:  WGS 84

Task 1:

# Filter out rows with zero population
final_data <- final_data %>%
  filter(totPop > 0)  

# Calculate the variables
final_data <- final_data %>%
  mutate(
    elderly_pct = (elderlyPop / totPop) * 100,       # Elderly population percentage
    case_rate = (positive / totPop) * 100000        # COVID-19 case rate per 100k
  )

# Map 1: COVID-19 Case Rates
covid_map <- ggplot(final_data) +
  geom_sf(aes(fill = case_rate), color = NA) +
  scale_fill_viridis_c(
    name = "Cases per 100k",
    option = "magma",
    limits = c(0, max(final_data$case_rate, na.rm = TRUE))
  ) +
  labs(
    title = "COVID-19 Case Rates by ZIP Code",
    subtitle = "New York City",
    caption = "Data Source: NYC Health Department"
  ) +
  coord_sf(datum = sf::st_crs(4326)) +  # Add graticule
  theme_void() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 12),
    plot.caption = element_text(size = 10),
    legend.position = "bottom"
  ) +
  ggspatial::annotation_scale(location = "br") +   # Add scale bar
  ggspatial::annotation_north_arrow(location = "tl", style = ggspatial::north_arrow_fancy_orienteering())

# Map 2: Elderly Population Percentage
elderly_map <- ggplot(final_data) +
  geom_sf(aes(fill = elderly_pct), color = NA) +
  scale_fill_distiller(
    name = "% Elderly Population",
    palette = "PuBuGn",
    direction = 1,
    limits = c(0, max(final_data$elderly_pct, na.rm = TRUE))
  ) +
  labs(
    title = "Elderly Population Percentage by ZIP Code",
    subtitle = "New York City",
    caption = "Data Source: US Census Bureau"
  ) +
  coord_sf(datum = sf::st_crs(4326)) +   # Add graticule
  theme_void() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 12),
    plot.caption = element_text(size = 10),
    legend.position = "bottom"
  ) +
  ggspatial::annotation_scale(location = "br") +   # Scale bar
  ggspatial::annotation_north_arrow(location = "tl", style = ggspatial::north_arrow_fancy_orienteering())

# Save maps as PNG files
ggsave("covid_case_rates_map.png", covid_map, width = 10, height = 8, dpi = 96) # Lower image quality since file size is >10 mb
ggsave("elderly_population_map.png", elderly_map, width = 10, height = 8, dpi = 96)

print(covid_map)

print(elderly_map)

Task 2:

# Calculate variables
final_data <- final_data %>%
  filter(totPop > 0) %>%  # Remove rows where total population is zero
  mutate(case_rate = (positive / totPop) * 100000) %>%  # Calculate case rate
  filter(!is.infinite(case_rate)) 

# Map 1: COVID-19 Case Rates
covid_map1 <- ggplot(final_data) +
  geom_sf(aes(fill = case_rate), color = NA) +
  scale_fill_gradient(
    name = "Cases per 100k",
    low = "#FEE0D2",   
    high = "#DE2D26",  
    limits = c(0, max(final_data$case_rate, na.rm = TRUE)) 
  ) +
  labs(
    title = "COVID-19 Case Rates by ZIP Code",
    subtitle = "New York City",
    caption = "Data Source: NYC Health Department"
  ) +
  coord_sf(datum = sf::st_crs(4326)) + 
  theme_void() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 10, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 8),
    plot.caption = element_text(size = 6),
    legend.position = "bottom"
  ) +
  ggspatial::annotation_scale(location = "br") + 
  ggspatial::annotation_north_arrow(location = "tl", style = ggspatial::north_arrow_fancy_orienteering())

# Map 2: Health Facility Counts
health_map <- ggplot(final_data) +
  geom_sf(aes(fill = health_facility_count), color = NA) +
  scale_fill_distiller(
    name = "Health Facilities",
    palette = "Blues",
    direction = 1,
    limits = c(0, max(final_data$health_facility_count, na.rm = TRUE))
  ) +
  labs(
    title = "Health Facility Counts by ZIP Code",
    subtitle = "New York City",
    caption = "Data Source: NYS Health Facility Data"
  ) +
  coord_sf(datum = sf::st_crs(4326)) +   # Graticule
  theme_void() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 10, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 8),
    plot.caption = element_text(size = 6),
    legend.position = "bottom"
  ) +
  ggspatial::annotation_scale(location = "br") +   # Add scale bar
  ggspatial::annotation_north_arrow(location = "tl", style = ggspatial::north_arrow_fancy_orienteering())

# Combine the maps 
multi_map <- covid_map1 + health_map + 
             plot_annotation(
               title = "COVID-19 Rates and Health Facility Distribution",
               subtitle = "Spatial Distribution Across NYC ZIP Codes",
               caption = "Data Sources: NYC Health Department & NYS Health Facility Data",
               theme = theme(
                 plot.title = element_text(hjust = 0.5, size = 18, face = "bold"),
                 plot.subtitle = element_text(hjust = 0.5, size = 14)
               )
             )

# Save as PNG
ggsave("multi_map_health_facilities.png", multi_map, width = 14, height = 8, dpi = 300)

print(multi_map)

Task 3:

final_data <- final_data %>%
  filter(totPop > 0 & positive >= 0) %>%  # Remove rows with zero population or negative cases
  mutate(case_rate = (positive / totPop) * 100000)  

covid_map2 <- leaflet(final_data) %>%
  addTiles() %>%  
  addPolygons(
    fillColor = ~colorNumeric("YlOrRd", case_rate)(case_rate),  
    weight = 1, 
    opacity = 1, 
    color = "white", 
    fillOpacity = 0.7,
    highlightOptions = highlightOptions(
      weight = 2,
      color = "#666",
      fillOpacity = 0.9,
      bringToFront = TRUE
    ),
    label = ~paste0(
      "ZIP Code: ", ZIPCODE, "<br>",
      "Case Rate: ", round(case_rate, 2), " per 100k"
    ),
    labelOptions = labelOptions(
      style = list("font-weight" = "bold", "color" = "black"),
      textsize = "12px",
      direction = "auto"
    )
  ) %>%
  addLegend(
    pal = colorNumeric("YlOrRd", final_data$case_rate),
    values = final_data$case_rate,
    title = "COVID-19 Case Rate<br>(per 100k)",
    position = "bottomright"
  )


saveWidget(covid_map2, file = "covid_interactive_map.html", selfcontained = FALSE)

# Display map 
covid_map2