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.

##Using the same code from the last assignment to pull the same points used for the aggregate layer

Article28_HealthFacilities_sf <- st_read(dsn = r"(H:\Hunter Class\GISr\Section_07\data\NYCData.gpkg)", 
                  layer = 'NYC Health Facilities') %>%
  filter(Short.Description == c('HOSP-EC', 'DTC', 'NH'))%>%
                   st_transform(2831)
## Reading layer `NYC Health Facilities' from data source 
##   `H:\Hunter Class\GISr\Section_07\data\NYCData.gpkg' using driver `GPKG'
## Simple feature collection with 1292 features and 34 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.19681 ymin: 40.51677 xmax: -73.69169 ymax: 40.91062
## Geodetic CRS:  WGS 84
## Warning: There was 1 warning in `stopifnot()`.
## ℹ In argument: `Short.Description == c("HOSP-EC", "DTC", "NH")`.
## Caused by warning in `Short.Description == c("HOSP-EC", "DTC", "NH")`:
## ! longer object length is not a multiple of shorter object length
mapview(st_geometry(Article28_HealthFacilities_sf))
breaks0423_eq <- classIntervals(CovidZipRate_sf$Week_2021_04_23PositiveRate, n = 7, style = "equal")
breaks0419_eq <- classIntervals(CovidZipRate_sf$Week_2020_04_19PositiveRate, n = 7, style = "equal")
breaks0412_eq <- classIntervals(CovidZipRate_sf$Week_2020_04_12PositiveRate, n = 7, style = "equal")

breaks0423_eq$brks[1] <- 0
breaks0419_eq$brks[1] <- 0
breaks0412_eq$brks[1] <- 0

CovidZipRate_sf <- mutate(CovidZipRate_sf, 
                          Week_2021_04_23PositiveRate_Cat = cut(Week_2021_04_23PositiveRate,
                                                                  breaks0423_eq$brks,
                                                                  dig.lab = 4, digits=1),
                          Week_2020_04_19PositiveRate_Cat = cut(Week_2020_04_19PositiveRate,
                                                                  breaks0419_eq$brks,
                                                                  dig.lab = 4, digits=1),
                          Week_2020_04_12PositiveRate_Cat = cut(Week_2020_04_12PositiveRate,
                                                                  breaks0412_eq$brks,
                                                                  dig.lab = 4, digits=1),
                          
)
gCovidMap <- ggplot(CovidZipRate_sf) + 
  geom_sf(aes(fill=Week_2021_04_23PositiveRate_Cat)) +
  geom_sf(data=Article28_HealthFacilities_sf, colour="green3", size=0.5) +
  scale_fill_brewer(palette = "OrRd", name='Postitive Rate') +
  coord_sf( default_crs = sf::st_crs(2831) ) +
  labs(x='Longitude', y='Latitude', 
       title='NYC Postive Covid Test Rate (April 2021)',
       caption = 
      'Points refers to location of article 28 health facilities')

gCovidMap

gHF28Map <- ggplot(HealthFacilities_sf) + 
  geom_sf(aes(fill=n_HF_cat)) +
  scale_fill_brewer(palette = "OrRd", name="Number of Facilities") +
  coord_sf(default_crs = sf::st_crs(2831) ) +
  labs(x="Longitude", y="Latitude", 
       title="Number of Health Facilities within a Zip Code")

gHF28Map

gCovid_HFMap <- ggarrange(gCovidMap +
                           theme(plot.title = element_text(size = 10),
                                 axis.title.x = element_blank(),
                                 axis.title.y = element_blank(),
                                 axis.text.x = element_blank(),
                                 axis.text.y = element_blank(),
                                 plot.caption = element_blank(),
                                 axis.ticks = element_blank(),
                                 legend.text = element_text(size = 8),
                                 legend.title = element_text(size = 9)
                                 ), 
                          gHF28Map +
                            theme(plot.title = element_text(size = 10),
                                 axis.title.x = element_blank(),
                                 axis.title.y = element_blank(),
                                 axis.text.x = element_blank(),
                                 axis.text.y = element_blank(),
                                 axis.ticks = element_blank(),
                                 legend.text = element_text(size = 8),
                                 legend.title = element_text(size = 9)
                                 ),
                          nrow = 1, ncol = 2, align="v") 

gCovid_HFMap

3. Create a web-based interactive map for COIVD-19 data using tmap, mapview, or leaflet package and save it as a HTML file.

CovidRate_3Week_MapView <- 
    mapView(x = CovidZipRate_sf, zcol="Week_2021_04_23PositiveRate_Cat", 
          layer.name = "Positive Test Rate (April 23 2021)", 
          homebutton=FALSE,
          col.regions = brewer.pal(9, "OrRd"),
          label="April 23 2021",
          popup = "April 23 2021") +
    mapView(x = CovidZipRate_sf, zcol="Week_2020_04_19PositiveRate_Cat", 
          layer.name = "Positive Test Rate (April 19 2020)", 
          homebutton=FALSE, 
          col.regions = brewer.pal(9, "OrRd"),
          label="April 19 2020",
          popup = "April 19 2020") +
  mapView(x = CovidZipRate_sf, zcol="Week_2020_04_12PositiveRate_Cat", 
          layer.name = "Positive Test Rate (April 12 2020)", 
          homebutton=FALSE, 
          col.regions = brewer.pal(9, "OrRd"),
          label="April 12 2020",
          popup = "April 12 2020")

CovidRate_3Week_MapView@map <- CovidRate_3Week_MapView@map %>%
  hideGroup(c("Positive Test Rate (April 19 2020)", "Positive Test Rate (April 12 2020)"))

CovidRate_3Week_MapView
mapview::mapshot(CovidRate_3Week_MapView,'3 Week Covid Rate.html')