R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

###########################################################################
# 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.

library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
library(patchwork)
## Warning: package 'patchwork' was built under R version 4.5.3
library(classInt)
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(ggplot2)
library(leaflet)

acsCovidDat <- sf::st_read('R-Spatial_III_Lab/acsPopByZip.gpkg', layer="NYC_COVID_Pop_by_Zip")
## Reading layer `NYC_COVID_Pop_by_Zip' from data source 
##   `C:\Users\ranae\Documents\Hunter\Spring 2026\GTECH 78520 - Data Analysis & Viz\R-Spatial\R-Spatial_III_Lab\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)
pal_bg <- RColorBrewer::brewer.pal(3, "Blues")

plot(acsCovidDat["Positiv"], 
     main = "Positive COVID Cases",
     breaks = "quantile", nbreaks = 3,
     graticule = st_crs(4326),
     axes = TRUE,
     pal = pal_bg,
     border = 'black',
     key.pos = 1)

plot(acsCovidDat["blackPp"], 
     main = "Black Population",
     breaks = "quantile", nbreaks = 3,
     graticule = st_crs(4326),
     axes = TRUE,
     pal = pal_bg,
     border = 'black',
     key.pos = 1) 

##########################################################################
# 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.).

# ggplot for covid cases
# creating quantiles
breaks_qt <- classIntervals(
  acsCovidDat$Positiv, 
  n = 3, 
  style = "quantile"
)

acsCovidDat <- acsCovidDat %>%
  mutate(
    Positiv_cat = cut(Positiv, breaks_qt$brks, include.lowest = TRUE)
  )

# plot
gg_covid <- ggplot(acsCovidDat) + 
  geom_sf(aes(fill = Positiv_cat), color = "black", linewidth = 0.2) +
  scale_fill_brewer(palette = "GnBu", name='Positive COVID Cases') +
  coord_sf(xlim = c(-74.05, -73.85), 
  ylim=c(40.68, 40.88), 
  default_crs = sf::st_crs(4326) ) +
  labs(x='Longitude', y='Latitude', 
       title='Positive COVID Cases')+
  theme(legend.position = "right");

gg_covid

# ggplot for black population
# creating quantiles
breaks_qt_bp <- classIntervals(
  acsCovidDat$blackPp, 
  n = 3, 
  style = "quantile"
)
## Warning in classIntervals(acsCovidDat$blackPp, n = 3, style = "quantile"): var
## has missing values, omitted in finding classes
acsCovidDat <- acsCovidDat %>%
  mutate(
    blackPp_cat = cut(blackPp, breaks_qt_bp$brks, include.lowest = TRUE)
  )

# plot
gg_blackpp <- ggplot(acsCovidDat) + 
  geom_sf(aes(fill = blackPp_cat), color = "black", linewidth = 0.2) +
  scale_fill_brewer(palette = "GnBu", name='Black Population') +
  coord_sf(xlim = c(-74.05, -73.85), 
  ylim=c(40.68, 40.88), 
  default_crs = sf::st_crs(4326) ) +
  labs(x='Longitude', y='Latitude', 
       title='Black Population')+
  theme(legend.position = "right");

gg_blackpp

gg_covid + gg_blackpp

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

# create palettes for each layer
covid_p <- colorQuantile('PuRd',NULL,n=3)
blackpp_p <- colorQuantile('BuPu',NULL,n=3)
elderlypp_p <- colorQuantile('BuPu',NULL,n=3)

# higlights
selectHighlight <- leaflet::highlightOptions(opacity = 0, fillColor = 'transparent')
hoverLabel <- leaflet::labelOptions(opacity = 1.0, textsize = '14px')

# popups
covid_popup <- paste0("<strong> Place: </strong>",
                  acsCovidDat$PO_NAME,
                  " <br/>",
                  "<strong> Confirmed COVID-19 Cases: </strong>",
                  acsCovidDat$Positiv,
                  sep="")
      
black_popup <- paste0("<strong> Place: </strong>",
                  acsCovidDat$PO_NAME,
                  " <br/>",
                  "<strong> Black Population: </strong>",
                  acsCovidDat$blackPp,
                  " <br/>",
                  "<strong> Confirmed COVID-19 Cases: </strong>",
                  acsCovidDat$Positiv,
                  sep="")
                  
elderly_popup <- paste0("<strong> Place: </strong>",
                  acsCovidDat$PO_NAME,
                  " <br/>",
                  "<strong> Elderly Population: </strong>",
                  acsCovidDat$eldrlyP,
                  " <br/>",
                  "<strong> Confirmed COVID-19 Cases: </strong>",
                  acsCovidDat$Positiv,
                  sep="")                 

# create map + layers
map <- leaflet(acsCovidDat %>% sf::st_transform(4326)) %>%
                addPolygons(
                      stroke=TRUE,
                      fillOpacity=1,
                      fillColor=~blackpp_p(blackPp),
                      color='black',
                      weight=1,
                      popup=black_popup,
                      highlightOptions=selectHighlight,
                      label=paste0(acsCovidDat$blackPp," black people in ",acsCovidDat$ZIPCODE),
                      group="Black Population",
                      labelOptions=hoverLabel) %>%
                addPolygons(
                      stroke=TRUE,
                      fillOpacity=1,
                      fillColor=~elderlypp_p(eldrlyP),
                      color='black',
                      weight=1,
                      popup=elderly_popup,
                      highlightOptions = selectHighlight,
                      label=paste0(acsCovidDat$eldrlyP, " elderly people in ", 
                      acsCovidDat$ZIPCODE),
                      group="Elderly Population",
                      labelOptions=hoverLabel) %>%
                addPolygons(
                      stroke=TRUE,
                      fillOpacity=1,
                      fillColor=~covid_p(Positiv),
                      color='black',
                      weight=1,
                      popup=covid_popup,
                      highlightOptions = selectHighlight,
                      label=paste0(acsCovidDat$eldrlyP, " positive COVID cases in ", 
                      acsCovidDat$ZIPCODE),
                      group="COVID cases",
                      labelOptions=hoverLabel) %>%      
              addProviderTiles("CartoDB.DarkMatter", group="DarkBG") %>%
              addProviderTiles("Esri.WorldStreetMap", group="WorldStreetMap") %>%
              addLayersControl(baseGroups=c("DarkBG", "WorldStreetMap"), 
                               overlayGroups = 
                               c("Black Population", 
                                 "Elderly Population",
                                 "COVID Polpulation")) 
                                 
map