# Read data created in HW08
zip_data<- st_read("C:/Users/linwa/Downloads/GTECH785_2025/R_Spatial/Data/R-Spatial_II_Lab/HW08_zip_data.gpkg")
## Reading layer `HW08_zip_data' from data source 
##   `C:\Users\linwa\Downloads\GTECH785_2025\R_Spatial\Data\R-Spatial_II_Lab\HW08_zip_data.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 248 features and 21 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25576 ymin: 40.49584 xmax: -73.6996 ymax: 40.91517
## Geodetic CRS:  WGS 84
glimpse(zip_data)
## Rows: 248
## Columns: 22
## $ ZIPCODE           <chr> "00083", "10001", "10002", "10003", "10004", "10005"…
## $ covid_case_count  <int> 0, 1542, 5902, 2803, 247, 413, 164, 379, 3605, 1686,…
## $ total_covid_test  <int> 0, 20158, 48197, 41076, 3599, 6102, 2441, 6342, 3877…
## $ COVID_DEATH_COUNT <int> 0, 35, 264, 48, 2, 0, 1, 4, 118, 37, 62, 14, 36, 32,…
## $ COUNTY            <chr> "New York", "New York", "New York", "New York", "New…
## $ PO_NAME           <chr> "Central Park", "New York", "New York", "New York", …
## $ CHHA              <int> NA, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, …
## $ DTC               <int> NA, 2, 3, 4, 1, 0, 0, 0, 0, 4, 3, 4, 4, 0, 3, 5, 0, …
## $ DTC.EC            <int> NA, 2, 5, 3, 0, 1, 0, 0, 3, 1, 1, 1, 2, 0, 0, 1, 3, …
## $ HOSP.EC           <int> NA, 0, 2, 4, 0, 0, 0, 1, 1, 5, 3, 0, 0, 1, 10, 1, 0,…
## $ HOSP.SB           <int> NA, 0, 2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, …
## $ HOSP              <int> NA, 0, 0, 3, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 2, 0, 0, …
## $ ADHCP             <int> NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ FoodStoreNum      <int> 1, 54, 199, 94, 13, 9, 5, 20, 84, 48, 111, 46, 116, …
## $ totPop            <int> 3, 19146, 74310, 53487, 1731, 8809, 4639, 15510, 586…
## $ malePctg          <dbl> 0.00000, 51.18040, 48.38783, 50.34868, 54.76603, 42.…
## $ asianPop          <int> 0, 4837, 32149, 8027, 358, 1974, 1195, 6040, 8503, 4…
## $ blackPop          <int> 2, 1092, 5969, 3130, 71, 185, 52, 1008, 5745, 1159, …
## $ hispanicPop       <int> 1, 2435, 20065, 4375, 150, 504, 507, 1096, 14672, 24…
## $ whitePop          <int> 1, 12485, 24033, 40503, 1207, 6419, 3130, 7912, 3595…
## $ elderlyPop        <int> 0, 2500, 15815, 6296, 52, 209, 66, 2134, 8907, 3469,…
## $ geom              <POLYGON [°]> POLYGON ((-73.94969 40.7970..., POLYGON ((-7…

Task 1

Plot at least two high-quality static maps with one using the COVID-19 data and one using a related factor.

plot(zip_data["covid_case_count"], main = "Total Covid Case Count by ZIP Code Area", 
     border = NA, 
     axes=TRUE,
     reset=FALSE,
     graticule = TRUE)

plot(zip_data["totPop"], main = "Total Population by ZIP Code Area", 
     border = NA, 
     axes=TRUE,
     reset=FALSE,
     graticule = TRUE)

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.

# Create labels for g1
labelCoords <- zip_data %>% sf::st_centroid() %>%
  dplyr::filter(COVID_DEATH_COUNT > 400) %>%
  sf::st_coordinates()
## Warning: st_centroid assumes attributes are constant over geometries
labelData <- zip_data %>% sf::st_centroid() %>%
  dplyr::filter(COVID_DEATH_COUNT > 400) %>%
  dplyr::mutate(x = labelCoords[,1], y=labelCoords[,2])
## Warning: st_centroid assumes attributes are constant over geometries
# Covid19 death count; you can also use geom_sf_label()
g1 <- ggplot(zip_data) + 
  geom_sf(aes(fill = COVID_DEATH_COUNT), show.legend = TRUE) +
  scale_fill_gradient(low = "yellow", high = "red", name = 'Covid Death Count') + 
  coord_sf(default_crs = sf::st_crs(4326)) +  
  geom_label_repel(data = labelData,
                   aes(x=x,y=y,label = COVID_DEATH_COUNT),
                   label.size = .09,
                   size = 3,
                   segment.color = rgb(0.5,0.5,0.9),
                   segment.size = 0.8)+
  labs(x = 'Longitude', y = 'Latitude', 
       title = 'COVID Death Count by ZIP Code') +
  theme(legend.position = "bottom",panel.grid.major = element_line(color = "gray", size = 0.5),  
        panel.grid.minor = element_line(color = "gray", size = 0.25),
        panel.grid = element_line(size = 0.5))
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Elderly population
g2 <- ggplot(zip_data) + 
  geom_sf(aes(fill = elderlyPop), show.legend = TRUE) +
  scale_fill_gradient(low = "yellow", high = "red", name = 'Elderly Population') + 
  coord_sf(default_crs = sf::st_crs(4326)) +  
  labs(x = 'Longitude', y = 'Latitude', 
       title = 'Elderly Population t by ZIP Code') +
  theme(legend.position = "bottom", 
        panel.grid.major = element_line(color = "gray", size = 0.5),  
        panel.grid.minor = element_line(color = "gray", size = 0.25),
        panel.grid = element_line(size = 0.5))

g1 + g2 + plot_layout(ncol = 2, nrow = 1)

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.

# Pop up
covid_popup <- paste0("<strong>Zip Code: </strong>", 
                  zip_data$ZIPCODE,
                   " <br/>",
                  "<strong> Place Name: </strong>",
                  zip_data$PO_NAME,
                  " <br/>",
                  "<strong> Total Population: </strong>",
                  zip_data$totPop,
                  " <br/>",
                   "<strong> COVID-19 Case Count: </strong>",
                  zip_data$covid_case_count,
                  " <br/>",
                  sep="")
# Breaks for the legend, using quantiles
breaks_qt <- classIntervals(c(min(zip_data$COVID_DEATH_COUNT) - .00001, zip_data$COVID_DEATH_COUNT), 
                            n = 5, style = "quantile")

# Define color palette based on the quantile breaks
color_pal <- colorBin(palette = "YlOrRd", 
                      domain = zip_data$COVID_DEATH_COUNT, 
                      bins = breaks_qt$brks, 
                      na.color = "transparent")

# Create map
htmlMap <- leaflet(zip_data %>% sf::st_transform(4326)) %>%
        addTiles(group = "OSM") %>%
        addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
        addPolygons(
                stroke = TRUE, # show the polygon boundary
                fillOpacity = 0.8,
                fillColor = ~color_pal(COVID_DEATH_COUNT), # Use color palette based on COVID_DEATH_COUNT
                color = 'grey', # color of the boundary
                weight = 1, # width of the boundary
                label = paste(zip_data$COVID_DEATH_COUNT),
                group = "Death Count",
                popup = covid_popup) %>%
  
        addLegend("bottomright", 
            pal = color_pal,  # Pass the color palette to the legend
            values = zip_data$COVID_DEATH_COUNT,  # Use actual values for the legend
            title = 'COVID-19 Death Count', 
            opacity = 1) %>%
        addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
        addLayersControl(baseGroups = c("OSM", "Carto"), 
                         overlayGroups = c("Death Count"))  

htmlMap