R Spatial Lab Assignment #3

task 1: Plot at least two high-quality static maps at the zip code level (COVID-19 & Elderly Population)

names(final_data)
##  [1] "ZIPCODE"           "BLDGZIP"           "PO_NAME"          
##  [4] "POPULATION"        "AREA"              "STATE"            
##  [7] "COUNTY"            "ST_FIPS"           "CTY_FIPS"         
## [10] "URL"               "SHAPE_AREA"        "SHAPE_LEN"        
## [13] "NEIGHBORHOOD_NAME" "BOROUGH_GROUP"     "label"            
## [16] "lat"               "lon"               "COVID_CASE_COUNT" 
## [19] "COVID_CASE_RATE"   "POP_DENOMINATOR"   "COVID_DEATH_COUNT"
## [22] "COVID_DEATH_RATE"  "PERCENT_POSITIVE"  "TOTAL_COVID_TESTS"
## [25] "food_store_count"  "nursinghome_count" "total_pop"        
## [28] "elderly_pop"       "white_pop"         "black_pop"        
## [31] "asian_pop"         "hispanic_pop"      "geometry"
final_data <- final_data %>%mutate(covid_quantile = cut(COVID_CASE_COUNT,
    breaks = quantile(COVID_CASE_COUNT, probs = seq(0, 1, 0.2), na.rm = TRUE),
    include.lowest = TRUE))

ggplot(final_data) + 
  geom_sf(aes(fill=covid_quantile))+
  scale_fill_brewer(palette = 'Blues', name='Per Zip Code') +
  coord_sf(xlim = c(-74.25, -73.70), default_crs = sf::st_crs(4326) ) +
  labs(x='Longitude', y='Latitude', 
       title='NYC COVID-19 Cases By Zip Code',
       caption = 'Data Source: NYC OpenData') +
  theme(legend.position = "right",panel.border = element_rect(colour = "black",fill=NA, size=1))

final_data <- final_data %>%
  mutate(elderly_quantile = cut(
    elderly_pop,
    breaks = unique(quantile(elderly_pop, probs = seq(0, 1, 0.2), na.rm = TRUE)),
    include.lowest = TRUE))

ggplot(final_data) + 
  geom_sf(aes(fill=elderly_quantile))+
  scale_fill_brewer(palette = 'Blues', name='Per Zip Code') +
  coord_sf(xlim = c(-74.25, -73.70), default_crs = sf::st_crs(4326) ) +
  labs(x='Longitude', y='Latitude', 
       title='NYC Elderly Population Count By Zipcode',
       caption = 'Data Source: NYC OpenData') +
  theme(legend.position = "right",panel.border = element_rect(colour = "black",fill=NA, size=1))

task 2: Create a multi-map illustrating the relationship between COVID-19 confirmed cases and Nursing Home Counts

final_data <- final_data %>%mutate(covid_quantile = cut(COVID_CASE_COUNT,
    breaks = quantile(COVID_CASE_COUNT, probs = seq(0, 1, 0.2), na.rm = TRUE),
    include.lowest = TRUE))

g3<-ggplot(final_data) + 
  geom_sf(aes(fill=covid_quantile))+
  scale_fill_brewer(palette = 'Blues', name='Per Zip Code') +
  coord_sf(xlim = c(-74.25, -73.70), default_crs = sf::st_crs(4326) ) +
  labs(x='Longitude', y='Latitude', 
       title='NYC COVID-19 Cases By Zip Code',
       caption = 'Data Source: NYC OpenData') +
  theme(legend.position = "right",panel.border = element_rect(colour = "black",fill=NA, size=1))

final_data <- final_data %>%
  mutate(NH_quantile = cut(
    nursinghome_count,
    breaks = unique(quantile(nursinghome_count, probs = seq(0, 1, 0.2), na.rm = TRUE)),
    include.lowest = TRUE))

g4<-ggplot(final_data) + 
  geom_sf(aes(fill = nursinghome_count)) +
  scale_fill_gradient(low = "lightblue",high = "darkblue",name = "Nursing Homes") +
  coord_sf(xlim = c(-74.25, -73.70), default_crs = sf::st_crs(4326)) +labs(
    x = 'Longitude', y = 'Latitude', 
    title = 'Nursing Homes by ZIP Code in NYC',
    caption = 'Data Source: NYC OpenData') +theme(legend.position = "right",
    panel.border = element_rect(colour = "black", fill = NA, size = 1))


g34 <- ggarrange(g3, g4, nrow = 1, ncol = 2, aligh="h") 
ggexport(g34, filename = file.path(getwd(), 'task2.pdf'),
           width=18.00, height=7.00, # the unit is inch
           pointsize = 4)

Nursing Homes and COVID Case Counts in NYC ### task 3: Create a web-based interactive map for COVID-19 data

names(final_data)
##  [1] "ZIPCODE"           "BLDGZIP"           "PO_NAME"          
##  [4] "POPULATION"        "AREA"              "STATE"            
##  [7] "COUNTY"            "ST_FIPS"           "CTY_FIPS"         
## [10] "URL"               "SHAPE_AREA"        "SHAPE_LEN"        
## [13] "NEIGHBORHOOD_NAME" "BOROUGH_GROUP"     "label"            
## [16] "lat"               "lon"               "COVID_CASE_COUNT" 
## [19] "COVID_CASE_RATE"   "POP_DENOMINATOR"   "COVID_DEATH_COUNT"
## [22] "COVID_DEATH_RATE"  "PERCENT_POSITIVE"  "TOTAL_COVID_TESTS"
## [25] "food_store_count"  "nursinghome_count" "total_pop"        
## [28] "elderly_pop"       "white_pop"         "black_pop"        
## [31] "asian_pop"         "hispanic_pop"      "geometry"         
## [34] "covid_quantile"    "elderly_quantile"  "NH_quantile"
final_data <- st_transform(final_data, 4326)

#Color palettes
covid_pal <- colorQuantile("YlOrRd", final_data$COVID_CASE_COUNT, n = 5)
death_pal <- colorQuantile("Reds", final_data$COVID_DEATH_RATE, n = 5)
positive_pal <- colorQuantile("Purples", final_data$PERCENT_POSITIVE, n = 5)

#Map tips
p_tip_cases <- paste0("ZIP: ", final_data$ZIPCODE)
p_tip_death <- paste0("ZIP: ", final_data$ZIPCODE)
p_tip_positive <- paste0("ZIP: ", final_data$ZIPCODE)

#Popups
p_popup_cases <- paste0("<strong>ZIP Code: </strong>", final_data$ZIPCODE, "<br/>",
  "<strong>COVID Cases: </strong>", final_data$COVID_CASE_COUNT)
p_popup_death <- paste0("<strong>ZIP Code: </strong>", final_data$ZIPCODE, "<br/>",
  "<strong>Death Rate: </strong>", round(final_data$COVID_DEATH_RATE, 2))
p_popup_positive <- paste0("<strong>ZIP Code: </strong>", final_data$ZIPCODE, "<br/>",
  "<strong>Percent Positive: </strong>", 
  paste0(round(final_data$PERCENT_POSITIVE, 1), "%"))

polyHighlightOption <- highlightOptions(opacity = 1,fillColor = "black")
polyLabelOption <- labelOptions(opacity = 0.6)

htmlMap <- leaflet(final_data) %>%
#COVID Cases Layer
  addPolygons(
    stroke = FALSE,
    fillColor = ~covid_pal(COVID_CASE_COUNT),
    fillOpacity = 0.8,
    smoothFactor = 0.5,
    popup = p_popup_cases,
    label = p_tip_cases,
    highlightOptions = polyHighlightOption,
    labelOptions = polyLabelOption,
    group = "COVID Cases") %>%
#COVID Death Rates Layer
  addPolygons(
    stroke = FALSE,
    fillColor = ~death_pal(COVID_DEATH_RATE),
    fillOpacity = 0.8,
    smoothFactor = 0.5,
    popup = p_popup_death,
    label = p_tip_death,
    highlightOptions = polyHighlightOption,
    labelOptions = polyLabelOption,
    group = "Death Rate") %>%
#Percent Positive Cases Layer
  addPolygons(
    stroke = FALSE,
    fillColor = ~positive_pal(PERCENT_POSITIVE),
    fillOpacity = 0.8,
    smoothFactor = 0.5,
    popup = p_popup_positive,
    label = p_tip_positive,
    highlightOptions = polyHighlightOption,
    labelOptions = polyLabelOption,
    group = "Percent Positive Cases") %>%
#Base Map
  addTiles(group = "OSM") %>%
  addProviderTiles("CartoDB.Positron", group = "Positron") %>%
  addProviderTiles("Esri.WorldImagery", group = "ESRI Imagery") %>%
#Legend for COVID Cases
  addLegend(
    "bottomright",
    pal = covid_pal,
    values = ~COVID_CASE_COUNT,
    title = "COVID Cases") %>%
  addLayersControl(
    baseGroups = c("OSM", "Positron", "ESRI Imagery"),
    overlayGroups = c("COVID Cases", "Death Rate", "Percent Positive"),
    options = layersControlOptions(collapsed = FALSE))

htmlwidgets::saveWidget(htmlMap, 'COVID_Data_leaflet.html', 
                        selfcontained = FALSE, libdir = 'widget-lib')
htmlMap

```