Explanation of the template

R Spatial Lab Assignment # 3

task 0: Import necessary libraries and data path

task 1: Plot the COVID-19 data Using Plot()

## Read COVID-19 file
covid_19_cases_df <- read.csv(file.path(data_dir, 'tests-by-zcta_2021_04_23.csv'))

## Because zip_code column "ZIPCODE" is character type, we convert zip code in covid-19 to char as well
## Reference: Convert method from https://stackoverflow.com/questions/57711323/convert-dataframe-numeric-column-to-character
covid_19_cases_df$MODIFIED_ZCTA<-as.character(covid_19_cases_df$MODIFIED_ZCTA)

## Filter NAs in COVID-19 file
covid_19_cases_df <- covid_19_cases_df %>%
                      ## remove the na values from the list, otherwise sf() cannot work
                      filter(!is.na(lon) , !is.na(lat))


## Read Zip Code Shapefile
zip_code_sf <- st_read(file.path(data_dir,'ZIP_CODE_040114.shp'))
## Reading layer `ZIP_CODE_040114' from data source 
##   `/Users/ziruizheng/Desktop/College Class/Now/Data Analysis & Viz/2025_Spring_week 7-10/Section_09/ZIP_CODE_040114.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 263 features and 12 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 913129 ymin: 120020.9 xmax: 1067494 ymax: 272710.9
## Projected CRS: NAD83 / New York Long Island (ftUS)
zip_code_sf <- st_transform(zip_code_sf, crs = 4326)


## Join COVID-19 data with zipcode spatial data
covid_19_cases_sf <- zip_code_sf %>%
  left_join(covid_19_cases_df, by = c("ZIPCODE" = "MODIFIED_ZCTA"))  


long_text <- '\n\nNotice:\n\nHere is a graph of COVID-19 data with Case Count.'
wrapped_text <- strwrap(long_text, width = 70)
cat(wrapped_text, sep = "\n")
## 
## 
## Notice:
## 
## Here is a graph of COVID-19 data with Case Count.
## Plot spatial map of COVID-19 cases
plot(covid_19_cases_sf["COVID_CASE_COUNT"], 
     main = "COVID-19 Cases in each zipcode", 
     breaks = "quantile", 
     border = NA)

task 1.2: Plot the Retail Food data Using Plot()

## Read file
## The fileEncoding method is from https://stackoverflow.com/questions/14363085/invalid-multibyte-string-in-read-csv
retail_food_df <- read.csv(file.path(data_dir,'nys_retail_food_store_xy.csv'), fileEncoding="latin1")

## Convert zip code from integer to char
retail_food_df$Zip.Code<-as.character(retail_food_df$Zip.Code)

## Count the number of store in each zip code
retail_food_df <- retail_food_df %>%
                  group_by(Zip.Code) %>%
                  summarise(number_of_store = n())

retail_food_sf <- zip_code_sf %>%
  left_join(retail_food_df, by = c("ZIPCODE" = "Zip.Code")) 


long_text <- '\n\nNotice:\n\nHere is a graph of Food store data with store number every zipcode'
wrapped_text <- strwrap(long_text, width = 70)
cat(wrapped_text, sep = "\n")
## 
## 
## Notice:
## 
## Here is a graph of Food store data with store number every zipcode
library(RColorBrewer)
plot(retail_food_sf["number_of_store"], 
     main = "Retail Food Stores number in each zip code",
     breaks = "quantile", nbreaks = 3,
     pal = brewer.pal(3, "OrRd"),
     border = NA)

task 2: ggplot2 to create multi-map figure between COVID-19 confirmed cases and number of food stores

## Merge the covid 19 and food store data
covid_19_cases_retail_food_sf_merge <- covid_19_cases_sf %>%
  left_join(retail_food_df, by = c("ZIPCODE" = "Zip.Code" )) 


# Clear out NAs in covid_19_cases_retail_food_sf_merge
covid_19_cases_retail_food_sf_merge <- covid_19_cases_retail_food_sf_merge %>%
                                      filter(!is.na(number_of_store), !is.na(COVID_CASE_COUNT))

## Set up for covid 19 multi map graph
require(classInt)
breaks_qt_food <- classIntervals(
  covid_19_cases_retail_food_sf_merge$number_of_store, 
  n = 7, 
  style = "quantile"
)


breaks_qt_covid <- classIntervals(
  covid_19_cases_retail_food_sf_merge$COVID_CASE_COUNT, 
  n = 7, 
  style = "quantile"
)

## Set up "store_number_cat"
covid_19_cases_retail_food_sf_merge <- mutate(covid_19_cases_retail_food_sf_merge, 
                            store_number_cat = cut(number_of_store, breaks_qt_food$brks, 
                                                 dig.lab = 4, digits=1)) 

## Set up "covid_cat"
covid_19_cases_retail_food_sf_merge <- mutate(covid_19_cases_retail_food_sf_merge, 
                            covid_cat = cut(COVID_CASE_COUNT, breaks_qt_covid$brks, 
                                                 dig.lab = 4, digits=1))


# Clear out NAs in covid_19_cases_retail_food_sf_merge, so NA map will not be generated in the fianl result
covid_19_cases_retail_food_sf_merge <- covid_19_cases_retail_food_sf_merge %>%
                                      filter(!is.na(covid_cat), !is.na(store_number_cat))


## Reference: Find color choices from Official website in r: https://ggplot2.tidyverse.org/reference/scale_brewer.html
## Reference: Change position of legend title from: https://r-graph-gallery.com/239-custom-layout-legend-ggplot2.html#:~:text=You%20can%20place%20the%20legend,the%20x%20and%20y%20coordinates.
## Reference: Find adding more rows to the legend from: https://stackoverflow.com/questions/27130610/legend-on-bottom-two-rows-wrapped-in-ggplot2-in-r

## Plot the multi-map graph
covid_19_map <- ggplot(covid_19_cases_retail_food_sf_merge) + 
  geom_sf(aes(fill=covid_cat), color = NA) +
  scale_fill_brewer(palette = "Purples", name='Covid 19 cases level') +
  guides(fill = guide_legend(nrow = 3))+
  labs(title='COVID-19 cases Per zip code',
      caption = 'Data Source: [...]') +
  theme_minimal() +
  theme(axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        legend.title.position = "top",
        legend.position = 'bottom')


food_store_map <- ggplot(covid_19_cases_retail_food_sf_merge) + 
  geom_sf(aes(fill=store_number_cat), color = NA) +
  scale_fill_brewer(palette = "Greens", name='Food Store Level') +
  guides(fill = guide_legend(nrow = 3))+
  labs(title='Food Store Level per zip code',
      caption = 'Data Source: [...]') +
  theme_minimal() +
  theme(axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        legend.title.position = "top",
        legend.position = 'bottom')


## Combine these two maps side by side
library(ggpubr)

combine_map <- ggarrange(covid_19_map, food_store_map, nrow = 1)

ggexport(combine_map, filename = file.path(getwd(), 'covid_food_stores_side_by_side.pdf'),
                       width=10, height=5, pointsize = 9); 
## file saved to /Users/ziruizheng/Desktop/College Class/Now/Data Analysis & Viz/2025_Spring_week 7-10/Section_09/covid_food_stores_side_by_side.pdf
combine_map

Task 3:

## Create a color palette for COVID-19 cases
pal_fun <- colorQuantile("YlOrRd", NULL, n = 5)


leaflet_map <- leaflet(covid_19_cases_sf) %>%
  addPolygons(
    stroke = FALSE, 
    fillColor = ~pal_fun(COVID_CASE_COUNT),
    fillOpacity = 0.8, smoothFactor = 0.5,
    popup = ~paste("Zipcode: ", ZIPCODE, "<br>COVID-19 Cases: ", COVID_CASE_COUNT)
  ) %>%
  addTiles() %>%
  addLegend(
    position = "bottomright",
    pal = pal_fun,
    values = ~COVID_CASE_COUNT,
    title = "COVID-19 Cases"
  ) 

leaflet_map
htmlwidgets::saveWidget(leaflet_map, "covid_leaflet.html")