R Spatial Lab Assignment #3
Task 0: Set up library continuing from Week 8
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.0 ✔ readr 2.1.6
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.1
## ── 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
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
library(mapview)
library(janitor)
##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(ggplot2)
library(dplyr)
library(ggspatial)
## Warning: package 'ggspatial' was built under R version 4.5.3
library(patchwork)
## Warning: package 'patchwork' was built under R version 4.5.3
library(htmlwidgets)
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.5.3
wd <- getwd()
data_dir <- file.path(wd,"data1")
zip_final <- st_read("data1/nyc_data_W8.gpkg", layer="zip_final") #Load Proper Data
## Reading layer `zip_final' from data source
## `C:\Users\wildk\OneDrive\Documents\2026 Spring\Visual in R\R-Spatial-Week9\data1\nyc_data_W8.gpkg'
## using driver `GPKG'
## Simple feature collection with 248 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.25576 ymin: 40.49584 xmax: -73.6996 ymax: 40.91517
## Geodetic CRS: WGS 84
health_fac <- st_read("data1/nyc_data.gpkg", layer="health")
## Reading layer `health' from data source
## `C:\Users\wildk\OneDrive\Documents\2026 Spring\Visual in R\R-Spatial-Week9\data1\nyc_data.gpkg'
## using driver `GPKG'
## Simple feature collection with 3848 features and 34 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -79.6299 ymin: -75.45935 xmax: 43.21162 ymax: 44.97849
## Geodetic CRS: WGS 84
colnames(zip_final) #Identifying the names of the columns
## [1] "ZIPCODE" "Positive" "Total" "food_store_count"
## [5] "health_count" "total_population" "covid_rate" "covid_pct"
## [9] "geom"
Task 1: Create Maps in the style of Static Choropleth Maps
# Map 1:
zip_final <- zip_final %>%
mutate(
covid_pct = ifelse(Total > 0,
(Positive / Total) * 100,
NA)
) #Covert positive count to covid-19 positive percentage for those who took a covid test
zip_final <- st_transform(zip_final, crs = 4326)
zip_clean <- zip_final %>%
filter(!is.na(covid_pct)) #Filter out NA values
map_covid <- ggplot(zip_clean) +
geom_sf(aes(fill = covid_pct), color = "white", size = 0.2) +
scale_fill_viridis_c(name = "COVID Rate Percentage") +
labs(
title = "COVID-19 Positivity Rate by Zip Code in NYC",
x = "Longitude",
y= "Latitude"
) +
coord_sf() +
theme_minimal() +
theme(
panel.grid.major = element_line(color = "gray60"),
panel.border = element_rect(color = "black", fill = NA, size = 1)
)
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Map 2:
map_DTC <- ggplot(zip_clean) +
geom_sf(aes(fill = health_count), color = "white", size = 0.2) +
scale_fill_viridis_c(name = "Diagnostic & Treatment Center Count") +
labs(
title = "Number of Diagnostic & Treatment Centers by Zip Code in NYC",
x = "Longitude",
y = "Latitude"
) +
coord_sf() +
theme_minimal() +
theme(
panel.grid.major = element_line(color = "grey60"),
panel.border = element_rect(color = "black", fill = NA, size = 1)
)
Task 3: Create an interactive map using the leaftlet package
map_leaflet <- leaflet(zip_clean) %>%
addTiles() #Base Map
map_leaflet <- map_leaflet %>%
addPolygons(
fillColor = ~colorNumeric("Reds", covid_pct)(covid_pct),
fillOpacity = 0.7,
color = "white",
weight = 1,
group = "COVID-19 Positive Rate (%)",
label = ~paste("ZIP:", ZIPCODE, "<br>Rate:", covid_pct),
popup = ~paste("COVID-19 Positive Rate (%):", covid_pct)
) #COVID-19 Layer
map_leaflet <- map_leaflet %>%
addPolygons(
fillColor = ~colorNumeric("Blues", health_count)(health_count),
fillOpacity = 0.7,
color = "white",
weight = 1,
group = "Diagnostic & Treatment Centers",
label = ~paste("Centers:", health_count),
popup = ~paste("Diagnostic & Treatment Centers:", health_count)
) #Diagnostic and Treatment Centers Layer
coords <- st_coordinates(health_fac) #Original health facility data points contain outlier points so they must be filterd out
health_fac_clean <- health_fac %>%
mutate(lon = coords[,1],
lat = coords[,2]) %>%
filter(lon >= -74.3 & lon <= -73.6) %>%
filter(lat >= 40.4 & lat <= 41.0) %>%
select(-lon, -lat) #Filter out data to just NYC so outlier points are not included
map_leaflet <- map_leaflet %>%
addCircleMarkers(
data = health_fac_clean,
radius = 3,
color = "lightgreen",
fillColor = "lightgreen",
group = "Facility Locations",
label = ~facility_name,
popup = ~facility_name
) #Facility Health Points Layer from cleaned data file
map_leaflet <- map_leaflet %>%
addLayersControl(
overlayGroups = c("COVID-19 Positive Rate (%)", "Diagnostic & Treatment Centers", "Facility Locations"),
options = layersControlOptions(collapsed = FALSE)
) #Control Layer
map_leaflet <- map_leaflet %>%
addLegend(
"bottomright",
pal = colorNumeric("Reds", zip_clean$covid_pct),
values = zip_clean$covid_pct,
title = "COVID-19 Positive Rate",
group = "COVID-19 Positive Rate"
) #Add Legend
map_leaflet ##Examine Interactive Map for errors
htmlwidgets::saveWidget(map_leaflet, "covid_interactive_map.html") #Save File
