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:

library(sf)
## Warning: package 'sf' was built under R version 4.4.3
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
library(RColorBrewer)
library(classInt)
## Warning: package 'classInt' was built under R version 4.4.3
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.4.3
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.4.3
library(htmlwidgets)
## Warning: package 'htmlwidgets' was built under R version 4.4.3
acs <- st_read("acsPopByZip.gpkg")
## Reading layer `NYC_COVID_Pop_by_Zip' from data source 
##   `C:\Users\alyss\Documents\R\week9\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)
names(acs)
##  [1] "ZIPCODE" "PO_NAME" "POPULAT" "COUNTY"  "Positiv" "Total"   "totPop" 
##  [8] "malPctg" "asianPp" "blackPp" "hspncPp" "whitePp" "eldrlyP" "geom"
head(acs)
## Simple feature collection with 6 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 971132.6 ymin: 188447.3 xmax: 992172.8 ymax: 215324.9
## Projected CRS: NAD83 / New York Long Island (ftUS)
##   ZIPCODE  PO_NAME POPULAT   COUNTY Positiv Total totPop  malPctg asianPp
## 1   10001 New York   22413 New York     260   571  25645 51.89706    6788
## 2   10002 New York   81305 New York     712  1358  76815 48.23277   32502
## 3   10003 New York   55878 New York     347   830  54181 50.34053    8027
## 4   10004 New York    2187 New York      24    64     NA       NA      NA
## 5   10005 New York    8107 New York      44   137   8809 42.81984    1974
## 6   10006 New York    3011 New York      14    54   2354 39.59218     640
##   blackPp hspncPp whitePp eldrlyP                           geom
## 1    1897    3679   16725    3239 MULTIPOLYGON (((981958.6 21...
## 2    7882   21621   27156   15770 MULTIPOLYGON (((991339.9 20...
## 3    3443    4375   42259    6391 MULTIPOLYGON (((989830.5 20...
## 4      NA      NA      NA      NA MULTIPOLYGON (((977542.5 18...
## 5     291     504    6603     206 MULTIPOLYGON (((982595.7 19...
## 6      48     245    1693      66 MULTIPOLYGON (((981136.3 19...
covid_map <- ggplot(acs) +
  geom_sf(aes(fill = Positiv), color = "gray60", linewidth = 0.2) +
  scale_fill_distiller(palette = "Reds", direction = 1, name = "COVID Cases") +
  labs(
    title = "COVID-19 Confirmed Cases by ZIP Code in NYC",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
    panel.grid.major = element_line(color = "gray85")
  )

covid_map

elderly_map <- ggplot(acs) +
  geom_sf(aes(fill = eldrlyP), color = "gray60", linewidth = 0.2) +
  scale_fill_distiller(palette = "Blues", direction = 1, name = "Elderly Pop.") +
  labs(
    title = "Elderly Population by ZIP Code in NYC",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
    panel.grid.major = element_line(color = "gray85")
  )

elderly_map

bb <- st_bbox(acs)

xmin <- bb["xmin"]
xmax <- bb["xmax"]
ymin <- bb["ymin"]
ymax <- bb["ymax"]
g1 <- ggplot(acs) +
  geom_sf(aes(fill = Positiv), color = "gray60", linewidth = 0.2) +
  scale_fill_distiller(
    palette = "Reds",
    direction = 1,
    name = "COVID Cases"
  ) +
  coord_sf(
    xlim = c(xmin, xmax),
    ylim = c(ymin, ymax),
    expand = FALSE
  ) +
  labs(
    title = "COVID-19 Confirmed Cases",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 13, hjust = 0.5, face = "bold"),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 9),
    legend.text = element_text(size = 8),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
    panel.grid.major = element_line(color = "gray85")
  )

g2 <- ggplot(acs) +
  geom_sf(aes(fill = eldrlyP), color = "gray60", linewidth = 0.2) +
  scale_fill_distiller(
    palette = "Blues",
    direction = 1,
    name = "Elderly Pop."
  ) +
  coord_sf(
    xlim = c(xmin, xmax),
    ylim = c(ymin, ymax),
    expand = FALSE
  ) +
  labs(
    title = "Elderly Population",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 13, hjust = 0.5, face = "bold"),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    legend.title = element_text(size = 9),
    legend.text = element_text(size = 8),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
    panel.grid.major = element_line(color = "gray85")
  )

# Put maps side by side
side_by_side_map <- ggarrange(
  g1, g2,
  ncol = 2,
  nrow = 1
)

side_by_side_map

ggsave(
  "side_by_side_map.png",
  plot = side_by_side_map,
  width = 12,
  height = 6,
  dpi = 300
)
acs4326 <- st_transform(acs, 4326)

pal_covid <- colorQuantile("Reds", acs4326$Positiv, n = 5)
pal_elderly <- colorQuantile("Blues", acs4326$eldrlyP, n = 5)

covid_popup <- paste0(
  "<strong>ZIP Code:</strong> ", acs4326$ZIPCODE,
  "<br><strong>Place Name:</strong> ", acs4326$PO_NAME,
  "<br><strong>Confirmed COVID Cases:</strong> ", acs4326$Positiv,
  "<br><strong>Total Tests:</strong> ", acs4326$Total
)

elderly_popup <- paste0(
  "<strong>ZIP Code:</strong> ", acs4326$ZIPCODE,
  "<br><strong>Place Name:</strong> ", acs4326$PO_NAME,
  "<br><strong>Elderly Population:</strong> ", acs4326$eldrlyP
)

map1 <- leaflet(acs4326) %>%
  addTiles(group = "OpenStreetMap") %>%
  addProviderTiles("CartoDB.DarkMatter", group = "Dark Map") %>%
  
  addPolygons(
    fillColor = ~pal_covid(Positiv),
    fillOpacity = 0.8,
    color = "gray50",
    weight = 1,
    popup = covid_popup,
    label = ~paste0("ZIP ", ZIPCODE, ": ", Positiv, " COVID cases"),
    group = "COVID Cases"
  ) %>%
  
  addPolygons(
    fillColor = ~pal_elderly(eldrlyP),
    fillOpacity = 0.6,
    color = "navy",
    weight = 1,
    popup = elderly_popup,
    label = ~paste0("ZIP ", ZIPCODE, ": ", eldrlyP, " elderly residents"),
    group = "Elderly Population"
  ) %>%
  
  addPolygons(
    fill = FALSE,
    color = "black",
    weight = 1,
    label = ~paste0("ZIP Code: ", ZIPCODE),
    group = "ZIP Boundaries"
  ) %>%
  
  addLegend(
    "bottomright",
    pal = pal_covid,
    values = ~Positiv,
    title = "COVID Cases",
    group = "COVID Cases"
  ) %>%
  
  addLayersControl(
    baseGroups = c("OpenStreetMap", "Dark Map"),
    overlayGroups = c("COVID Cases", "Elderly Population", "ZIP Boundaries"),
    options = layersControlOptions(collapsed = FALSE)
  )
highlightOptions = highlightOptions(
  weight = 2,
  color = "black",
  bringToFront = TRUE
)

map1
htmlwidgets::saveWidget(map1, "lab9_interactive_map.html", selfcontained = FALSE)