#read in required packages
require(tidyverse); require(magrittr);
## Loading required package: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: magrittr
## 
## 
## Attaching package: 'magrittr'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     set_names
## 
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
require(sp);require(sf);require(rgdal);
## Loading required package: sp
## Loading required package: sf
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
## Loading required package: rgdal
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'rgdal'
require(classInt);require(RColorBrewer);
## Loading required package: classInt
## Loading required package: RColorBrewer
require(ggplot2);require(ggmap);require(leaflet);
## Loading required package: ggmap
## ℹ Google's Terms of Service: <https://mapsplatform.google.com>
##   Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service/>
##   OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles/>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
## Attaching package: 'ggmap'
## 
## The following object is masked from 'package:magrittr':
## 
##     inset
## 
## Loading required package: leaflet
require(ggrepel); require(ggsn);
## Loading required package: ggrepel
## Loading required package: ggsn
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'ggsn'
library(ggpubr); library(devtools);
## Loading required package: usethis
require(htmlwidgets); require(mapview)
## Loading required package: htmlwidgets
## Loading required package: mapview
library(classInt)
acsCovidDat <- sf::st_read('./Data/acsPopByZip.gpkg', layer="NYC_COVID_Pop_by_Zip")
## Reading layer `NYC_COVID_Pop_by_Zip' from data source 
##   `/Users/jamesreo/Desktop/R Independent Study/Weeks 7-9/Session 9/R-spatial/Data/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)
  1. Plot at least two high-quality static maps with one using the COVID-19 data and one using a related factor. You can use either plot method for sf or ggplot method.
#simple plot for visual aid
plot(acsCovidDat)
## Warning: plotting the first 9 out of 13 attributes; use max.plot = 13 to plot
## all

ggplot(data = acsCovidDat) +
  geom_sf(aes(fill = Positiv)) +
  scale_fill_viridis_c(option = "viridis", na.value = "white", guide = guide_legend(title = "COVID-19 Cases")) +
  labs(title = "COVID-19 BY NYC Zipcode")

ggplot(data = acsCovidDat) +
  geom_sf(aes(fill = eldrlyP)) +
  scale_fill_viridis_c(option = "magma", na.value = "white", guide = guide_legend(title = "Population")) +
  labs(title = "Elderly Population BY NYC Zipcode")

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 (e.g., the number of nursing homes, number of food stores, neighborhood racial composition, elderly population, etc.). The maps should be put side by side on one single page. Add graticule to at least one of those maps and label some of the feature on the map where applicable and appropriate.

# Convert acsCovidData to long format for 'Positiv' and 'eldrlyP'
acsCovidDat_long <- acsCovidDat %>% 
  pivot_longer(cols = c(Positiv, eldrlyP), names_to = "variable", values_to = "value") %>% 
  st_as_sf()

# Calculate breaks for classification
breaks_lmh <- classIntervals(var = acsCovidDat_long$value, n = 3, style = "quantile", na.rm = TRUE)$brks
## Warning in classIntervals(var = acsCovidDat_long$value, n = 3, style =
## "quantile", : var has missing values, omitted in finding classes
# Adjust the first break slightly 
breaks_lmh[1] <- breaks_lmh[1] - 0.0001

# Classify 'value' into Low, Medium, High categories
acsCovidDat_long$valueCategory <- cut(acsCovidDat_long$value, breaks = breaks_lmh, labels = c("Low", "Medium", "High"))

# Now, create the plot with facets for 'Positiv' and 'eldrlyP', colored by 'valueCategory'

ggplot(acsCovidDat_long) +
  geom_sf(aes(fill = valueCategory)) +
  scale_fill_brewer(palette = "OrRd", name = "Category") +
  labs(title = 'COVID-19 Positive Cases and Elderly Population by Zip Code', 
       caption = 'Data Source: [...]',
       x = "Longitude",  
       y = "Latitude") +  
  facet_wrap(~variable) +
  coord_sf(datum = st_crs(4326)) +  
  theme_minimal() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1),  
        axis.text.y = element_text())

  1. Create a web-based interactive map for COIVD-19 data using tmap, mapview, or leaflet package and save it as a HTML file.
pal_fun <- colorQuantile("YlOrRd", NULL, n = 5)

pal_totalpop <- colorQuantile("YlGnBu", 
                        domain = c(min(acsCovidDat$totPop, na.rm = T), max(acsCovidDat$totPop, na.rm = T)),
                        alpha = TRUE)

covid_popup <- paste0("<strong>Zip Code: </strong>", 
                  acsCovidDat$ZIPCODE,
                   " <br/>",
                  "<strong> Place Name: </strong>",
                  acsCovidDat$PO_NAME,
                  " <br/>",
                  "<strong> Confirmed COVID-19 Cases: </strong>",
                  acsCovidDat$Positiv,
                  sep="")

polyHighlightOption <- leaflet::highlightOptions(opacity = 1.0, fillColor = 'black')
polyLabelOption <- leaflet::labelOptions(opacity = 1.0, textsize = '11px')

htmlMap <- leaflet(acsCovidDat %>% sf::st_transform(4326)) %>%
        addPolygons(
                stroke = TRUE, # show the polygon boundary
                fillOpacity = 0.8,
                fillColor = ~pal_totalpop(totPop),
                color = 'grey', # color of the boundary
                weight = 1, # width of the boundary
                highlightOptions = polyHighlightOption,
                label = paste(acsCovidDat$totPop),
                group = "Total Population") %>%
        addPolygons(
                stroke = FALSE,
                fillColor = ~pal_fun(Positiv),
                fillOpacity = 0.8, smoothFactor = 0.5,
                popup = covid_popup,
                group = "NYC-COVID",
                label = paste0(acsCovidDat$Positiv, " at ", acsCovidDat$PO_NAME, ", ", acsCovidDat$ZIPCODE),
                highlightOptions = polyHighlightOption,
                labelOptions = polyLabelOption) %>%
        addTiles(group = "OSM") %>%
        addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
        addLayersControl(baseGroups = c("OSM", "Carto"), 
                         overlayGroups = c("NYC-COVID", "Total Population"))  

htmlMap
htmlwidgets::saveWidget(htmlMap, 'nyc_acs_covid_leaflet.html')