Code setup

R Spatial Lab Assignment # 3

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

#Reading in files
nyczip <- st_read("data/R-Spatial_I_Lab/ZIP_CODE_040114/ZIP_CODE_040114.shp")
## Reading layer `ZIP_CODE_040114' from data source 
##   `/Users/joelbeckwith/Documents/Hunter_College/RVizGeospatial/R-spatial/data/R-Spatial_I_Lab/ZIP_CODE_040114/ZIP_CODE_040114.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 263 features and 12 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 913129 ymin: 120020.9 xmax: 1067494 ymax: 272710.9
## Projected CRS: NAD83 / New York Long Island (ftUS)
covidcsv <- readr::read_csv("data/R-Spatial_II_Lab/tests-by-zcta_2020_04_12.csv", lazy = FALSE)
## Rows: 178 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (4): MODZCTA, Positive, Total, zcta_cum.perc_pos
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
covidzip_merged <- base::merge(nyczip, covidcsv, by.x = "ZIPCODE", by.y = "MODZCTA")
nyczip <- st_transform(nyczip, st_crs(4326))

nyshealth <- readr::read_csv("data/R-Spatial_I_Lab/NYS_Health_Facility.csv",
                             show_col_types = FALSE, 
                             lazy = FALSE)
nychealth_clean <- nyshealth %>% 
  filter(!is.na(`Facility Longitude`), !is.na(`Facility Latitude`),`Regional Office` == "Metropolitan Area Regional Office - New York City",!`Facility ID` %in% c(10383, 10324,10306),`Short Description` == "NH")
nychealth_sf <- st_as_sf(nychealth_clean,
                         coords = c("Facility Longitude","Facility Latitude"),
                         crs = 4326)

#plotting
library(RColorBrewer)
pal <- brewer.pal(5, "Reds")
plot(covidzip_merged["Total"],
     main="Covid Cases by Census Area, Week of April 22, 2020",
     breaks="quantile", nbreaks = 5,
     pal= pal,
     border = NA,
     graticule = st_crs(4326),
     axes = TRUE,
     reset = TRUE)

ptColor <- rgb(0.1, 0.6, 1, 0.5)
plot(st_geometry(nyczip),
     main="Nursing Home Facilities in NYC",
     graticule = st_crs(4326),
     axes=TRUE,
     lwd=0.3,
     reset=FALSE)
plot(nychealth_sf["Facility ID"], pch=19, col=ptColor,cex=0.3, add=TRUE)

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.

library(ggplot2)
require(classInt)
## Loading required package: classInt
library(dplyr)
library(grid)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
#covid map
breaks_qt <- classIntervals(c(min(covidzip_merged$Total) - .00001,
                              covidzip_merged$Total), n = 7, style = "quantile")

covidzip_merged <- mutate(covidzip_merged, 
                            covid_rate_cat = cut(Total, breaks_qt$brks, 
                                                 dig.lab = 4, digits=1))

covidmap <- ggplot(covidzip_merged) + 
  geom_sf(aes(fill=covid_rate_cat)) +
  scale_fill_brewer(palette = "OrRd", name='Cases by Zip Code') +
  coord_sf(xlim = c(-74.3, -73.65), default_crs = sf::st_crs(4326)) +
  labs(x='Longitude', y='Latitude', 
       title='NYC Covid Totals by Zip Code',
       caption = 'Data Source: [NYC Department of Health]') +
  theme(
    panel.grid.major = element_line(color = "gray80", size = 0.5),
    panel.grid.minor = element_line(color = "gray90", size = 0.5) 
  )
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#nursing home map
nychealth_sf <- nychealth_sf %>%
   rename(ZIPCODE = `Facility Zip Code`)

nh_count_by_zip <- nychealth_sf %>%
  st_drop_geometry() %>%
  group_by(ZIPCODE) %>%
  summarise(nh_count = n())

nyczip_nh <- nyczip %>%
  left_join(nh_count_by_zip, by = "ZIPCODE") %>%
  mutate(nh_count = replace_na(nh_count, 0))


nhmap <- ggplot(nyczip_nh) + 
  geom_sf(aes(fill = nh_count)) +
  scale_fill_viridis_c(option = "plasma", name = "Nursing Homes") +
  coord_sf(xlim = c(-74.3, -73.65)) +
  labs(title = "Nursing Homes by ZIP Code in NYC",
       x = "Longitude", y = "Latitude",
       caption = "Data Source: NYS DOH") 

grid.arrange(covidmap, nhmap, ncol = 2)
grid.text("Multi-map: Covid Cases & Nursing Homes in NYC", x = 0.5, y = 1, gp = gpar(fontsize = 16, fontface = "bold"))

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.

library(webshot)

mvcovid_map <- mapview(covidzip_merged["Total"], layer.name = "Covid Cases")

mapshot(mvcovid_map, file = "covid_map.html")