R Spatial Lab Assignment # 3

Task 1:

Plot at least two high-quality static maps at the zip code level, with one using the COVID-19 testing data and one using a related factor such as the number of nursing homes or the elderly population. You can use either plot method for sf or ggplot method. - These maps must be choropleth maps, i.e., polygon filling colors representing quantities. - Maps must have title, legend, graticule (grid lines), neatline (border), and labels on axes (commonly geography coordinates in degrees). - Bonus point: add scale bar and north arrow.

  1. Setup - acsCovidData reading acsPopByZip.gpkg provided for this week’s data, sand setting layer NYC_COVID_Pop_by_Zip.
  2. Setting the two breaks for breaks_covid and beaks_elderly using classIntervals() via lecture.
    • Used - 0.00001 so that the lowest value is captured in the first bin.
    • Added breaks_elderly$brks[1] <- 0 to manually set the lower bound to 0 so that the legend isn’t a decimal.
  3. I mutated the columns for acsCovidData using cut() fixed scientific notation in legend.
  4. st_bbox for bounding box in WGS84 so coord_sf() axis labels display as degrees for both maps.
  5. Plot of map_covid and map_elderly
    • na.value = "grey80" - remove data with missing zip codes
    • datum = sf::st_crs(4326) - graticule grid via ggplot2 coord_sf
acsCovidData <- sf::st_read('./acsPopByZip.gpkg', 
                            layer = "NYC_COVID_Pop_by_Zip")
## Reading layer `NYC_COVID_Pop_by_Zip' from data source 
##   `C:\Users\moniq\Desktop\SPRING_26\DATA_ANALYSIS_R\Week_9\R-Spatial\acsPopByZip.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 180 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 913129 ymin: 120020.9 xmax: 1067113 ymax: 272710.9
## Projected CRS: NAD83 / New York Long Island (ftUS)
# Choropleth Map breaks - Positive covid casesa and elderly 
breaks_covid <- classIntervals(
  c(min(acsCovidData$Positiv, na.rm = TRUE) - 0.00001, acsCovidData$Positiv),
  n = 7, style = "quantile"
)

breaks_elderly <- classIntervals(
  c(min(acsCovidData$eldrlyP, na.rm = TRUE) - 0.00001, acsCovidData$eldrlyP),
  n = 7, style = "quantile"
)
## Warning in classIntervals(c(min(acsCovidData$eldrlyP, na.rm = TRUE) - 1e-05, :
## var has missing values, omitted in finding classes
# Forcing the first break to 0 so legend isnt in negatives for lower bounds
breaks_elderly$brks[1] <- 0

acsCovidData <- acsCovidData %>%
  mutate(
    covid_cat = cut(Positiv,breaks_covid$brks,dig.lab = 4),
    elderly_cat = cut(eldrlyP,breaks_elderly$brks,dig.lab = 6))

nyc_bbox <- sf::st_bbox(sf::st_transform(acsCovidData, 4326))

# covid map
map_covid <- ggplot(acsCovidData) +
  geom_sf(aes(fill = covid_cat), color = "white", linewidth = 0.1) +
  scale_fill_brewer(palette = "YlOrRd", name = "Positive Cases", na.value = "grey80") +
  coord_sf(
    xlim = c(nyc_bbox["xmin"], nyc_bbox["xmax"]),
    ylim = c(nyc_bbox["ymin"], nyc_bbox["ymax"]),
    default_crs = sf::st_crs(4326),
    datum = sf::st_crs(4326)  
  ) +
  annotation_scale(location = "bl", width_hint = 0.3) +
  annotation_north_arrow(location = "tr", which_north = "true",
  style = north_arrow_fancy_orienteering()) +
  labs(
    x = "Longitude",
    y = "Latitude",
    title = "NYC COVID-19 Positive Cases by ZIP Code",
    caption = "Data Source: NYC DOH, 2020"
  ) +
  theme_bw() +
  theme(
    legend.position = "right",
    panel.border = element_rect(color = "black", linewidth = 1.2),
    panel.grid = element_line(color = "grey80", linewidth = 0.3)
  )

map_covid

# Elderly pop map
map_elderly <- ggplot(acsCovidData) +
  geom_sf(aes(fill = elderly_cat), color = "white", linewidth = 0.1) +
  scale_fill_brewer(palette = "Blues", name = "Elderly Pop (65+)", na.value = "grey80") +
  coord_sf(
    xlim = c(nyc_bbox["xmin"], nyc_bbox["xmax"]),
    ylim = c(nyc_bbox["ymin"], nyc_bbox["ymax"]),
    default_crs = sf::st_crs(4326),
    datum = sf::st_crs(4326)
  ) +
  annotation_scale(location = "bl", width_hint = 0.3) +
  annotation_north_arrow(location = "tr", which_north = "true",
  style = north_arrow_fancy_orienteering()) +
  labs(
    x = "Longitude",
    y = "Latitude",
    title = "NYC Elderly Population (65+) by ZIP Code",
    caption = "Data Source: ACS 2018, US Census"
  ) +
  theme_bw() +
  theme(
    legend.position = "right",
    panel.border = element_rect(color = "black", linewidth = 1.2),
    panel.grid = element_line(color = "grey80", linewidth = 0.3)
  )

map_elderly

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, neighborhood racial composition, elderly population, etc.). - The maps should be put side by side on one single page. Both maps must have the exactly same size and scale.

  1. Modified both maps to reomve the tite and caption, setting both legends below the maps. Also removed the y-axis since they share the same latitude range.
  2. combined map is both plots where plot_layout() to control the grid, and plot_annotation() is the shared title/caption. 
map_covid_left <- map_covid +
  labs(title = NULL, caption = NULL) +
  theme(legend.position = "bottom", 
        legend.title.position = "top")

map_elderly_right <- map_elderly +
  labs(title = NULL, caption = NULL, y = "") +
  theme(legend.position = "bottom",    
    legend.title.position = "top", 
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank())

# combining maps covid positive left and elderly pop to the right
combined <- map_covid_left + map_elderly_right +
  plot_layout(ncol = 2, nrow = 1) +
  plot_annotation(
    title = "NYC COVID-19 Positive Cases and Elderly Population by ZIP Code",
    caption = "Data Sources: NYC DOH 2020; ACS 2018, US Census",
    theme = theme(plot.title = element_text(size = 14, hjust = 0.5))
  )

combined

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. - The interactive map must have at least three layers (three layers in a single interactive map, not three layers in three separate interactive maps.) - At least one of the three layers has legend (opens and closes together with the corresponding layer). - Layers should have tips (show up when a user moves mouse over a feature like a point or polygon) and popups (when a user click on a feature). - Bonus point: the interactive map is implemented with leaflet and/or extra mapping features (explicitly state what extra features you added). Features implemented in packages by default would not count. For example, mapview would show the current map zoom level and the coordinates of the mouse by default.

  1. Reproject data to EPSG:4326 as acsCovidData_4326, added health_sf_nyc from lab 2 as a point layer because I wanted to see points with the polygons, added a filter for hospital data.
  2. Set palette colors for all four maps.
  3. Using paste0 to display specific pop-up data for all four datasets per map.
  4. Leaflet map creation - event loop
      1. Added addEasyButton using a font awesome crosshair icon to re-center map using an onClick() handler function and map.setView(). R docs easyButton
      1. addPolygons() - changed palette, popup and column names per dataset.
      1. Basemap options via leaflet-providers
      1. Legends - Added a legend for each map and positioned them in separate areas to prevent overlaps.
      1. addLayersControl() via lecture but also added hideGroup() to set covid layer by default and a toggle for the other two maps visibility.
# Adding health facilities as point layer
health_sf_nyc <- st_read("nyc_data.gpkg", layer = "health_facilities")
## Reading layer `health_facilities' from data source 
##   `C:\Users\moniq\Desktop\SPRING_26\DATA_ANALYSIS_R\Week_9\R-Spatial\nyc_data.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1293 features and 34 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 929529 ymin: 127612.4 xmax: 1066447 ymax: 271059.5
## Projected CRS: NAD83 / New York Long Island (ftUS)
hospitals <- health_sf_nyc %>%
  filter(Short.Description == 'HOSP')

# Reproject data for leaflet 
acsCovidData_4326 <- sf::st_transform(acsCovidData, 4326)
health_sf_nyc_4326 <- sf::st_transform(hospitals, 4326)

# Palettes for covid positive, elderly pop and asian pop
pal_covid_pos <- colorQuantile("YlOrRd", NULL, n = 5)

pal_elderly <- colorQuantile("YlGnBu", 
      domain = c(min(acsCovidData_4326$eldrlyP, na.rm = TRUE),
      max(acsCovidData_4326$eldrlyP, na.rm = TRUE)),
      alpha = TRUE)

pal_asian <- colorQuantile("Purples",
      domain = c(min(acsCovidData_4326$asianPp, na.rm = TRUE),
      max(acsCovidData_4326$asianPp, na.rm = TRUE)),
      alpha = TRUE)

# Pop-up labels 
covid_popup <- paste0(
  "<strong>ZIP Code: </strong>",acsCovidData_4326$ZIPCODE,"<br/>",
  "<strong>Place Name: </strong>",acsCovidData_4326$PO_NAME,"<br/>",
  "<strong>Confirmed Cases: </strong>",acsCovidData_4326$Positiv,"<br/>",
  "<strong>Total Tests: </strong>",acsCovidData_4326$Total)

elderly_popup <- paste0(
  "<strong>ZIP Code: </strong>",acsCovidData_4326$ZIPCODE,"<br/>",
  "<strong>Place Name: </strong>",acsCovidData_4326$PO_NAME,"<br/>",
  "<strong>Elderly Pop (65+): </strong>",acsCovidData_4326$eldrlyP)

asian_popup <- paste0(
  "<strong>ZIP Code: </strong>",acsCovidData_4326$ZIPCODE,"<br/>",
  "<strong>Place Name: </strong>",acsCovidData_4326$PO_NAME,"<br/>",
  "<strong>Asian Population: </strong>",acsCovidData_4326$asianPp)

hospital_popup <- paste0("<strong>Name: </strong>",health_sf_nyc_4326$Facility.Name,"<br/>",
  "<strong>Address: </strong>",health_sf_nyc_4326$Facility.Address.1,"<br/>",
  "<strong>Phone Number: </strong>",health_sf_nyc_4326$Facility.Phone.Number)

# Highlights and labels
polyHighlightOption <- leaflet::highlightOptions(opacity = 1.0, fillColor = "black")
polyLabelOption <- leaflet::labelOptions(opacity = 1.0, textsize = "11px")

# The Leaflet map

htmlMap <- leaflet(acsCovidData_4326) %>%
# re-center button at top left
   addEasyButton(easyButton(
    icon = "fa-crosshairs",
    title = "Reset View",
    onClick = JS("function(btn, map){ map.setView([40.70, -73.98], 10); }")
  )) %>%

# Covid positive layer
  addPolygons(
    stroke = FALSE,
    fillColor = ~pal_covid_pos(Positiv),
    fillOpacity = 0.8, smoothFactor = 0.5,
    popup = covid_popup,
    label = paste0(
      acsCovidData_4326$Positiv, " cases at ",
      acsCovidData_4326$PO_NAME, ", ",
      acsCovidData_4326$ZIPCODE),
    highlightOptions = polyHighlightOption,
    labelOptions = polyLabelOption,
    group = "NYC-COVID"
  ) %>%

# Elderly Pop layer
  addPolygons(
    stroke = TRUE, color = "grey", weight = 1,
    fillColor = ~pal_elderly(eldrlyP),
    fillOpacity = 0.8,
    popup = elderly_popup,
    label = paste0("Elderly Pop: ", acsCovidData_4326$eldrlyP,
                         " | ", acsCovidData_4326$PO_NAME, ", ",
                         acsCovidData_4326$ZIPCODE),
    highlightOptions = polyHighlightOption,
    labelOptions = polyLabelOption,
    group = "Elderly Population"
  ) %>%

# Asian Pop layer
  addPolygons(
    stroke = TRUE, color = "grey", weight = 1,
    fillColor = ~pal_asian(asianPp),
    fillOpacity = 0.8,
    popup = asian_popup,
    label = paste0(
      "Asian Pop: ", acsCovidData_4326$asianPp,
      " | ", acsCovidData_4326$PO_NAME, ", ",
      acsCovidData_4326$ZIPCODE),
    highlightOptions = polyHighlightOption,
    labelOptions = polyLabelOption,
    group = "Asian Population"
  ) %>%

    # Hospital point layer
  addCircleMarkers(
      data = health_sf_nyc_4326,
      radius = 4,
      color = "blue",
      fillOpacity = 0.7,
      stroke = FALSE,
      popup = hospital_popup,
      label = paste0("Hospital Name: ", health_sf_nyc_4326$Facility.Name),
      group = "Hospitals") %>%
  
# Basemaps
  addTiles(group = "OSM") %>%
  addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
  addProviderTiles("CartoDB.Positron", group = "Positron") %>%

# Legend Covid Positive - Bottom right
  addLegend("bottomright",
      pal = pal_covid_pos,
      values = ~Positiv,
      title = "COVID-19 Positive <br> Cases by ZIP Code",
      group = "NYC-COVID") %>%

# Legend Elderly Pop - Bottom left
  addLegend("bottomleft",
      pal = pal_elderly,
      values = ~eldrlyP,
      title = "Elderly Population <br> (65+)",
      group = "Elderly Population") %>%

# Legend Asian Pop - Top right
  addLegend("topright",
      pal = pal_asian,
      values = ~asianPp,
      title = "Asian Population",
      group = "Asian Population") %>%
  

  addLayersControl(
    baseGroups = c("Positron", "OSM", "Carto"),
    overlayGroups = c("NYC-COVID", "Hospitals", "Elderly Population", "Asian Population"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  hideGroup(c("Hospitals", "Elderly Population", "Asian Population"))

htmlMap
htmlwidgets::saveWidget(htmlMap, "nyc_covid_interactive.html", selfcontained = TRUE)