libraries

library(sf)
## Warning: package 'sf' was built under R version 4.2.3
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(tidyverse)
## Warning: package 'tibble' was built under R version 4.2.3
## Warning: package 'dplyr' was built under R version 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1
## ── 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
library(mapview)
## Warning: package 'mapview' was built under R version 4.2.3
library(ggplot2)
library(dplyr)
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.2.3
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.2.3
library(leaflet.providers)
## Warning: package 'leaflet.providers' was built under R version 4.2.3
library(RColorBrewer)
# Set working directory
wd <- dirname(rstudioapi::getActiveDocumentContext()$path)
setwd(wd)
acsPopByZip <- st_read("Data/R-Spatial_III_Lab/acsPopByZip.gpkg")
## Reading layer `NYC_COVID_Pop_by_Zip' from data source 
##   `C:\Users\amyca\OneDrive\Documents\GTECH7852_R\R-spatial\GTECH7852_HW\Data\R-Spatial_III_Lab\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)
HW11 <- st_read("Data/R-Spatial_III_Lab/HW11.gpkg")
## Reading layer `NYC_Demographics' from data source 
##   `C:\Users\amyca\OneDrive\Documents\GTECH7852_R\R-spatial\GTECH7852_HW\Data\R-Spatial_III_Lab\HW11.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 175 features and 35 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25576 ymin: 40.49584 xmax: -73.70098 ymax: 40.91517
## Geodetic CRS:  WGS 84
healthfacilities <- st_read("Data/R-Spatial_III_Lab/healthfacilities.gpkg")
## Reading layer `NYC_HealthFacilties_SF' from data source 
##   `C:\Users\amyca\OneDrive\Documents\GTECH7852_R\R-spatial\GTECH7852_HW\Data\R-Spatial_III_Lab\healthfacilities.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 830 features and 37 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.19681 ymin: 0 xmax: 0 ymax: 40.90462
## Geodetic CRS:  WGS 84

Task 1:

Task 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.

GGPLOT METHOD: Multimap 1

Positive COVID 19 Cases (4/12/2020 vs. 4/19/2020)

# Creating breaks
br <- c(0, 25, 50, 75, 100) 

HW11$X4_12_20_PP<- cut(HW11$X4_12_20_cum_percpositive,
                              breaks = br,
                              dig.lab = 3)
HW11$X4_19_20_PP<- cut(HW11$X4_19_20_cum_percpositive,
                              breaks = br,
                              dig.lab = 3)
HW11$X4_23_21_PP<- cut(HW11$X4_23_21_cum_percpositive,
                              breaks = br,
                              dig.lab = 3)
# Creating a palette
pal <- hcl.colors(6, "Inferno", rev = TRUE, alpha = 0.7)


plot1 <- ggplot() +
  geom_sf(data = HW11, 
          aes(fill= X4_12_20_PP)) +
          labs(
               title='Percent Positive 4/12/2020') +
          theme(legend.position = "none") +
          # Custom palette
          scale_fill_manual(values = pal,
                            drop = FALSE,
                            na.value = "grey80") 

plot2 <-  ggplot() + geom_sf(data = HW11, 
          aes(fill= X4_19_20_PP)) +
          labs( 
               title='Percent Positive 4/19/2020')+
          theme(legend.position = "none") +
          # Custom palette
          scale_fill_manual(values = pal,
                            drop = FALSE,
                            na.value = "grey80")+
          geom_sf_label(data=HW11 %>%
                filter(X4_12_20_cum_percpositive > 75),
                aes(label = X4_12_20_cum_percpositive),
                label.size = .01,
                size = 1.5)

grid.arrange(
  plot1,
  plot2,
  ncol = 2) 
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

### GGPLOT METHOD: Multimap 2 #### COVID-19 Cases and Health Facilties

#clean
HW11$Total_HF <- HW11$Total_HealthFacilities
#palette
pal <- hcl.colors(6, "Inferno", rev = TRUE, alpha = 0.7)


p1 <-  ggplot() +
  geom_sf(data = acsPopByZip,
          aes(fill= Positiv)) +
          labs(
               title='Positive COVID-19 Cases')+
          #theme(legend.position = "none") +
          scale_fill_viridis_c(direction = -1)
          # Custom palette
          #scale_fill_manual(values = pal,
                            #drop = FALSE,
                            #na.value = "grey80") 
p2 <-  ggplot() + 
  geom_sf(data = HW11, 
          aes(fill= Total_HF)) +
          labs(
               title='Health Facilties') +
          #theme(legend.title = HF) + 
          # Custom palette
          scale_fill_viridis_c(direction = -1) 


grid.arrange(
  p1,
  p2,
  ncol = 2)

Task 3:

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

Interactive Web Map (Leaflet)

COVID 19 Positive Cases & Health Facilties

# Creating a palete
pal_fun <- colorQuantile("OrRd", NULL, n = 5)

# Quick data clean-up
healthfacilities <- st_transform(healthfacilities, crs=4326) %>%
  filter(facility_latitude > 0)
acsPopByZip <- st_transform(acsPopByZip, crs=4326)

# Leaflet
polyHighlightOption <- leaflet::highlightOptions(opacity = 1.0, fillColor = 'black')
polyLabelOption <- leaflet::labelOptions(opacity = 0.6)
p_tip <- paste0("Positive Cases: ", acsPopByZip$Positiv);

htmlMap <- leaflet(acsPopByZip) %>%
  addPolygons(
    stroke = FALSE, 
    fillColor = ~pal_fun(Positiv),
    fillOpacity = 0.8, smoothFactor = 0.5,
    highlightOptions = polyHighlightOption,
    label = p_tip,
    labelOptions = polyLabelOption) %>%
  addCircles(radius = 0.02, weight = 1,
             data = healthfacilities %>% 
               sf::st_set_crs(4326) %>% sf::st_cast("POINT"),
             group = "Health Fac.") %>%
  addTiles(group = "OSM") %>%
  addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
  addLayersControl(baseGroups = c("Carto", "OSM"), 
                   overlayGroups = c("COVID19", "HF")) %>%
  addLegend("bottomright",  # location
          pal=pal_fun,    # palette function
          values=~Positiv,  # value to be passed to palette function
          title = 'COVID-19 Positive Cases <br> & Health Facilities') # legend title


htmlMap