R Spatial Lab Assignment #3

task 1: Choropleth Maps

load("~/R-Spatial/COVID_POP_ZIP_NYC.RData")

names(COVID_POP_ZIP_NYC)
##  [1] "ZIPCODE"            "PO_NAME"            "POPULATION"        
##  [4] "COUNTY"             "POSITIVE_COVID"     "TOTAL_TESTS"       
##  [7] "FoodStoreNum"       "HealthFacilityNum"  "totPop"            
## [10] "malePctg"           "asianPop"           "blackPop"          
## [13] "hispanicPop"        "whitePop"           "geometry"          
## [16] "POSITIVE_COVID_cat" "positive_per_1000"  "rate_cat"          
## [19] "covid_cat"          "latino_cat"
breaks_covid <- classIntervals(
  c(min(COVID_POP_ZIP_NYC$POSITIVE_COVID, na.rm = TRUE) - .00001,
    COVID_POP_ZIP_NYC$POSITIVE_COVID),
  n = 7,
  style = "quantile"
)
## Warning in classIntervals(c(min(COVID_POP_ZIP_NYC$POSITIVE_COVID, na.rm = TRUE)
## - : var has missing values, omitted in finding classes
map_covid <- ggplot(COVID_POP_ZIP_NYC) +
  geom_sf(aes(fill = POSITIVE_COVID_cat), color = "white", linewidth = 0.2) +
  scale_fill_brewer(palette = "OrRd", direction=1, name = "Positive COVID") +
  coord_sf(default_crs = st_crs(4326)) +
  labs(
    x = "Longitude",
    y = "Latitude",
    title = "Positive COVID-19 Tests by ZIP Code",
    caption = "Data Source: NYC COVID and ZIP data"
  ) +
  annotation_scale(location = "bl", width_hint = 0.3) +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    style = north_arrow_fancy_orienteering,
  ) +
  theme_bw() +
  theme(
    legend.position = "right",
    panel.border = element_rect(color = "black", fill = NA),
    panel.background = element_rect(fill = "gray97", color = NA),
    panel.grid.major = element_line(color = "gray60", linewidth = 0.5),
  panel.grid.minor = element_line(color = "gray80", linewidth = 0.3)
  )
map_covid

Health Facities only include hospitals, hospital extentions and nursing homes

map_health_facilities <- ggplot(COVID_POP_ZIP_NYC) +
  geom_sf(aes(fill = HealthFacilityNum), color = "white", linewidth = 0.2) +
  scale_fill_viridis_c(name = "Health Facilities") +
  coord_sf(default_crs = st_crs(4326)) +
  labs(
    title = "Health Facilities by ZIP Code in NYC",
    x = "Longitude",
    y = "Latitude"
  ) +
  annotation_north_arrow(location = "tr", which_north = "true") +
  annotation_scale(location = "bl") +
  theme_minimal() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
    plot.title = element_text(face = "bold", size = 14),
    panel.grid.major = element_line(color = "gray60", linewidth = 0.5),
  panel.grid.minor = element_line(color = "gray80", linewidth = 0.3)
  )

map_health_facilities

task 2: Multi-Map

covid_vals <- COVID_POP_ZIP_NYC$POSITIVE_COVID
covid_vals <- covid_vals[is.finite(covid_vals)]



breaks_covid <- classIntervals(
  c(min(covid_vals) - 0.00001, covid_vals),
  n = 7,
  style = "quantile"
)

COVID_POP_ZIP_NYC <- COVID_POP_ZIP_NYC %>%
  mutate(
    covid_cat = cut(
      POSITIVE_COVID,
      breaks = breaks_covid$brks,
      include.lowest = TRUE,
      dig.lab = 4
    )
  )

map3 <- ggplot(COVID_POP_ZIP_NYC) +
  geom_sf(aes(fill = covid_cat), color = "white", linewidth = 0.2) +
  scale_fill_brewer(palette = "OrRd", name = "Positive COVID\n(count)") +
  coord_sf(default_crs = NULL, lims_method = "geometry_bbox") +
  labs(
    title = "Positive COVID-19 Tests by ZIP Code",
    x = "Longitude",
    y = "Latitude"
  ) +
  annotation_scale(location = "bl",
  width_hint = 0.25,
  text_cex = 2.5)+
  annotation_north_arrow(location = "tr", which_north = "true",height = unit(4, "cm"),
  width = unit(4, "cm")) +
  theme_bw() +
  theme(
    panel.border = element_rect(color = "black", fill = NA),
    panel.grid.major = element_line(color = "gray80", linewidth = 0.5),
  panel.grid.minor = element_line(color = "gray80", linewidth = 0.3),
   plot.title = element_text(size = 80, face = "bold"),
  axis.title = element_text(size = 60),
  axis.text = element_text(size = 40),
  legend.title = element_text(size = 60),
  legend.text = element_text(size = 45),
  legend.key.size = unit(2.5, "cm")
  )

latino_vals <- COVID_POP_ZIP_NYC$hispanicPop
latino_vals <- latino_vals[is.finite(latino_vals)]

breaks_latino <- classIntervals(
  c(min(latino_vals) - 0.00001, latino_vals),
  n = 7,
  style = "quantile"
)

COVID_POP_ZIP_NYC <- COVID_POP_ZIP_NYC %>%
  mutate(
    latino_cat = cut(
      hispanicPop,
      breaks = breaks_latino$brks,
      include.lowest = TRUE,
      dig.lab = 4
    )
  )
map4 <- ggplot(COVID_POP_ZIP_NYC) +
  geom_sf(aes(fill = latino_cat), color = "white", linewidth = 0.2) +
  scale_fill_brewer(palette = "Blues", name = "Latino Population") +
  coord_sf(default_crs = NULL, lims_method = "geometry_bbox") +
  labs(
    title = "Latino Population by ZIP Code",
    x = "Longitude",
    y = "Latitude"
  ) +
  annotation_scale(location = "bl",
  width_hint = 0.25,
  text_cex = 2.5) +
  annotation_north_arrow(location = "tr", which_north = "true",height = unit(4, "cm"),
  width = unit(4, "cm")) +
  theme_bw() +
  theme(
    panel.border = element_rect(color = "black", fill = NA),
    panel.grid.major = element_line(color = "gray80", linewidth = 0.5),
  panel.grid.minor = element_line(color = "gray80", linewidth = 0.3),
  plot.title = element_text(size = 80, face = "bold"),
  axis.title = element_text(size = 60),
  axis.text = element_text(size = 40),
  legend.title = element_text(size = 60),
  legend.text = element_text(size = 45),
  legend.key.size = unit(2.5, "cm")
  )
library(patchwork)

final_maps3_4 <- map3 + map4
final_maps3_4

TASK 3: Web-based Interactive Map

covid_map_dat <- COVID_POP_ZIP_NYC %>%
  st_transform(4326)

# color palettes
pal_covid <- colorNumeric("YlOrRd", domain = covid_map_dat$POSITIVE_COVID, na.color = "transparent")
pal_latino <- colorNumeric("Blues", domain = covid_map_dat$hispanicPop, na.color = "transparent")
pal_health <- colorNumeric("Greens", domain = covid_map_dat$HealthFacilityNum, na.color = "transparent")

polyLabelOption <- labelOptions(
  style = list("font-weight" = "normal", padding = "3px 8px"),
  textsize = "13px",
  direction = "auto"
)

# highlight options
polyHighlightOption <- highlightOptions(
  weight = 3,
  color = "black",
  bringToFront = TRUE
)

# popups
covid_popup <- paste0(
  "<strong>ZIP Code:</strong> ", covid_map_dat$ZIPCODE,
  "<br><strong>Neighborhood:</strong> ", covid_map_dat$PO_NAME,
  "<br><strong>Positive COVID Tests:</strong> ", covid_map_dat$POSITIVE_COVID
)

latino_popup <- paste0(
  "<strong>ZIP Code:</strong> ", covid_map_dat$ZIPCODE,
  "<br><strong>Neighborhood:</strong> ", covid_map_dat$PO_NAME,
  "<br><strong>Latino Population:</strong> ", covid_map_dat$hispanicPop
)

health_popup <- paste0(
  "<strong>ZIP Code:</strong> ", covid_map_dat$ZIPCODE,
  "<br><strong>Neighborhood:</strong> ", covid_map_dat$PO_NAME,
  "<br><strong>Health Facilities:</strong> ", covid_map_dat$HealthFacilityNum
)

htmlMap <- leaflet(covid_map_dat) %>%
  addTiles(group = "OpenStreetMap") %>%
  addProviderTiles("CartoDB.DarkMatter", group = "Carto Dark") %>%
  
  addPolygons(
    stroke = TRUE,
    color = "gray40",
    weight = 1,
    fillOpacity = 0.8,
    fillColor = ~pal_covid(POSITIVE_COVID),
    group = "COVID-19 Positive Tests",
    label = ~paste0(PO_NAME, " (", ZIPCODE, "): ", POSITIVE_COVID),
    popup = lapply(covid_popup, HTML),
    highlightOptions = polyHighlightOption,
    labelOptions = polyLabelOption
  ) %>%
  
  addPolygons(
    stroke = TRUE,
    color = "gray40",
    weight = 1,
    fillOpacity = 0.8,
    fillColor = ~pal_latino(hispanicPop),
    group = "Latino Population",
    label = ~paste0(PO_NAME, " (", ZIPCODE, "): ", hispanicPop),
    popup = lapply(latino_popup, HTML),
    highlightOptions = polyHighlightOption,
    labelOptions = polyLabelOption
  ) %>%
  
  addPolygons(
    stroke = TRUE,
    color = "gray40",
    weight = 1,
    fillOpacity = 0.8,
    fillColor = ~pal_health(HealthFacilityNum),
    group = "Health Facilities",
    label = ~paste0(PO_NAME, " (", ZIPCODE, "): ", HealthFacilityNum),
    popup = lapply(health_popup, HTML),
    highlightOptions = polyHighlightOption,
    labelOptions = polyLabelOption
  ) %>%
  
  addLayersControl(
    baseGroups = c("OpenStreetMap", "Carto Dark"),
    overlayGroups = c("COVID-19 Positive Tests", "Latino Population", "Health Facilities"),
    options = layersControlOptions(collapsed = FALSE),
    position="topright"
  ) %>%
  
  addLegend(
    pal = pal_covid,
    values = ~POSITIVE_COVID,
    title = "Positive COVID Tests",
    position = "topleft",
    group = "COVID-19 Positive Tests"
  ) %>%
  
  addLegend(
    pal = pal_latino,
    values = ~hispanicPop,
    title = "Latino Population",
    position = "bottomright",
    group = "Latino Population"
  ) %>%
  
  addLegend(
    pal = pal_health,
    values = ~HealthFacilityNum,
    title = "Health Facilities",
    position = "bottomleft",
    group = "Health Facilities"
  ) %>%
  
### EXTRA FEATURE
  addControl(
    html = "<b>NYC COVID-19 & Demographics Map</b>",
    position = "topright"
  )
tmlMap <- htmlMap %>%
  hideGroup(c("Latino Population", "Health Facilities"))

htmlMap