R Spatial Lab Assignment # 9

task 1:

nys_health_sf <- st_read("nys_health_sf.gpkg", quiet = TRUE)

covid_data <-  st_read("covidPopZipNYC.gpkg", quiet = TRUE)

str(covid_data)
## Classes 'sf' and 'data.frame':   180 obs. of  20 variables:
##  $ ZIPCODE             : chr  "10001" "10002" "10003" "10004" ...
##  $ PO_NAME             : chr  "New York" "New York" "New York" "New York" ...
##  $ POPULATION          : num  22413 81305 55878 2187 8107 ...
##  $ COUNTY              : chr  "New York" "New York" "New York" "New York" ...
##  $ NEIGHBORHOOD_NAME   : chr  "Chelsea/NoMad/West Chelsea" "Chinatown/Lower East Side" "East Village/Gramercy/Greenwich Village" "Financial District" ...
##  $ COVID_CASE_COUNT    : num  1542 5902 2803 247 413 ...
##  $ TOTAL_COVID_TESTS   : num  20158 48197 41076 3599 6102 ...
##  $ COVID_DEATH_COUNT   : num  35 264 48 2 0 1 4 118 37 62 ...
##  $ FoodStoreNum        : int  57 192 87 10 10 4 24 84 52 100 ...
##  $ healthfacilities_NUM: int  7 20 16 1 1 1 2 5 12 13 ...
##  $ totPop              : int  19146 74310 53487 NA 8809 4639 15510 58683 28129 56433 ...
##  $ adultPop            : int  17658 64499 49620 NA 7853 4096 12954 52081 25531 51006 ...
##  $ malePctg            : num  51.2 48.4 50.3 NA 42.8 ...
##  $ femalePop           : num  48.8 51.6 49.7 NA 57.2 ...
##  $ elderlyPop          : int  2500 15815 6296 NA 209 66 2134 8907 3469 8848 ...
##  $ asianPop            : int  4837 32149 8027 NA 1974 1195 6040 8503 4702 7301 ...
##  $ blackPop            : int  1092 5969 3130 NA 185 52 1008 5745 1159 3220 ...
##  $ hispanicPop         : int  2435 20065 4375 NA 504 507 1096 14672 2434 6828 ...
##  $ whitePop            : int  12485 24033 40503 NA 6419 3130 7912 35952 21438 43728 ...
##  $ geom                :sfc_MULTIPOLYGON of length 180; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:110, 1:2] 981959 981980 981988 981993 982052 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geom"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:19] "ZIPCODE" "PO_NAME" "POPULATION" "COUNTY" ...
hist(covid_data$COVID_CASE_COUNT,
     main = "Histogram of New COVID-19 Cases",
     col = "skyblue",
     border = "white")

breaks_j <- classIntervals(c(min(covid_data$COVID_CASE_COUNT) - .00001,
                              covid_data$COVID_CASE_COUNT), n = 7, style = "jenks")

covid_data <- covid_data %>% 
  mutate(covid_data_cut = cut(COVID_CASE_COUNT, breaks_j$brks, 
                             include.lowest = TRUE, dig.lab = 4))

clean_labels <- sapply(1:(length(breaks_j$brks) - 1), function(i) {
  paste0(
    format(breaks_j$brks[i], big.mark = ",", scientific = FALSE, trim = TRUE),
    " – ",
    format(breaks_j$brks[i + 1], big.mark = ",", scientific = FALSE, trim = TRUE)
  )
})

ggplot(covid_data) + 
  geom_sf(aes(fill = covid_data_cut)) +
  scale_fill_brewer(palette = "Purples", 
                    name = ' Covid Cases per Zipcode', 
                    labels = clean_labels) +
  labs(
    x = 'Longitude', y = 'Latitude',
    title = 'New York City Covid Cases in 2021',
    caption = 'Data Source: The City of New York'
  ) +
  theme(
    plot.caption = element_text(
      hjust = 0.5,    
      color = "grey60",   
      size = 8,           
      margin = margin(t = 5) 
    )
  ) 

breaks_j_pop <- classIntervals(c(min(covid_data$elderlyPop) - .00001,
                              covid_data$elderlyPop), n = 5, style = "jenks")
## Warning in classIntervals(c(min(covid_data$elderlyPop) - 1e-05,
## covid_data$elderlyPop), : var has missing values, omitted in finding classes
covid_data <- covid_data %>% 
  mutate(elderlyPop_cut = cut(elderlyPop, breaks_j_pop$brks, 
                             include.lowest = TRUE, dig.lab = 4))
    


clean_labels_pop <- sapply(1:(length(breaks_j_pop$brks) - 1), function(i) {
  paste0(
    format(breaks_j_pop$brks[i], big.mark = ",", scientific = FALSE, trim = TRUE),
    " – ",
    format(breaks_j_pop$brks[i + 1], big.mark = ",", scientific = FALSE, trim = TRUE))
    }
  )



ggplot(covid_data) + 
  geom_sf(aes(fill = elderlyPop_cut)) +
  scale_fill_brewer(palette = "OrRd", 
                    name = 'Elderly Population',
                    labels = clean_labels_pop) +
  labs(
    x = 'Longitude', y = 'Latitude',
    title = 'Elderly Population in New York City',
    caption = 'Data Source: American Community Survey'
  ) +
  theme(
    plot.caption = element_text(
      hjust = 0.5,    
      color = "grey60",   
      size = 8,           
      margin = margin(t = 5) 
    )
  ) 

task 2:

g1 <- ggplot(covid_data) + 
  geom_sf(aes(fill = covid_data_cut), show.legend = TRUE) +
  scale_fill_brewer(palette = "OrRd", name = 'Covid Cases per Zipcode', labels = clean_labels) +
  coord_sf(
    xlim = c(-74.26, -73.70),
    ylim = c(40.48, 40.92),
    default_crs = sf::st_crs(4326)
  ) + 
  geom_sf_label(data = covid_data %>% dplyr::filter(COVID_CASE_COUNT > 10000),
                aes(label = COVID_CASE_COUNT),
                label.size = .07,
                size = 2) +
  labs(
    x = 'Longitude', 
    y = 'Latitude', 
    title = 'New York City Covid Cases in 2021',
    caption = 'Data Source: The City of New York'
  ) +
  theme(
    legend.position = "bottom",  
    legend.justification = "center",
    legend.box.just = "center",  
    legend.title = element_text(hjust = 1), 
    plot.caption = element_text(
      hjust = 0,    
      color = "grey60",   
      size = 8,           
      margin = margin(t = 5) 
    )) + 
  guides(fill = guide_legend(
    title.position = "top", 
  )) 


g2 <- ggplot(nys_health_sf) + 
  geom_sf(colour='blue', size=0.1) +
  scale_fill_brewer(palette = "OrRd", name='Noise Complaints') +
  coord_sf(xlim = c(-74.26, -73.70),
           ylim = c(40.48, 40.92), default_crs = sf::st_crs(4326))  +
  labs(x='Longitude', y='Latitude', 
       title='NYC Health Facilities')

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

task 3:

covid_data <- st_transform(covid_data, 4326)

pal_fun <- colorQuantile("Purples", NULL, n = 5)

p_popup <- paste0(
  "<strong>Covid Cases: </strong>", covid_data$COVID_CASE_COUNT, 
  "<br/>",
  "<strong>Covid Deaths: </strong>", covid_data$COVID_DEATH_COUNT
)

polyHighlightOption <- leaflet::highlightOptions(opacity = 0.8, fillColor = 'grey')
polyLabelOption <- leaflet::labelOptions(opacity = 0.6)

htmlMap <- leaflet(covid_data) %>%
  addTiles(group = "OSM") %>%
  addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
  addPolygons(
    stroke = FALSE,
    fillColor = ~pal_fun(COVID_CASE_COUNT),
    fillOpacity = 0.8,
    group = "Covid Cases",
    smoothFactor = 0.5,
    popup = p_popup,
    highlightOptions = polyHighlightOption,
    labelOptions = polyLabelOption
  ) %>%
  addCircleMarkers(radius = 0.03, weight = 1,
             data = nys_health_sf %>% 
               sf::st_set_crs(4326) %>% sf::st_cast("POINT"),
             group = "Health Facilities")%>%
  addLegend(
    "bottomright",
    colors = brewer.pal(7, "Purples"),
    labels = paste0("up to ", format(breaks_j$brks[-1], digits = 1)),
    title = "Covid Cases Per Zipcode"
  ) %>%
  addLayersControl(baseGroups = c("OSM", "Carto"), 
                   overlayGroups = c("Covid Cases", "Health Facilities")) %>%
  hideGroup(c("Health Facilities")) %>%
  setView(lng = -73.93, lat = 40.73, zoom = 10)

htmlMap
htmlwidgets::saveWidget(htmlMap, 
                       file = "nyc_covid_map.html",
                       selfcontained = TRUE)