Packages and Loading

library(tidyverse);
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ 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
library(sf); 
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(mapview); 
library(magrittr);
## 
## Attaching package: 'magrittr'
## 
## The following object is masked from 'package:purrr':
## 
##     set_names
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
library(RColorBrewer)
library(classInt)
library(leaflet)

load("Session_8_data.RData")

First Map

COVID cases in NYC map

pal<-brewer.pal(8,"YlOrRd")

plot(all1["COVID_CASE_COUNT"],
     main="COVID Cases in NYC",
     breaks="jenks",nbreaks=8,
     pal=pal,
     graticule = st_crs(4326),
     axes=TRUE,
     reset=FALSE)

Next Map

Map of the number of Diagnostic and Treatment Centers in NYC

pal2<-brewer.pal(5,"Blues")

plot(all1["HFNum"],
     main="Number of Diagnostic and Treatment Centers in NYC",
     breaks="jenks", nbreaks=5,
     pal=pal2,
     graticule = st_crs(4326),
     axes=TRUE,
     reset=FALSE)

GGPlot Maps

Map of COVID cases and black population including locations of diagnostic centers

ptColor <- rgb(0.1,0.1,0.9,0.7)
dtcenters<- hf_clean %>%
  filter(Short.Description=='DTC',
         Facility.County %in% c('New York','Kings','Richmond','Bronx','Queens')) %>%
  st_as_sf(.,coords=c('Facility.Longitude','Facility.Latitude'))

st_crs(dtcenters)<-4326
st_crs(dtcenters)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
st_crs(all1)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
breaksjk<-classIntervals(c(min(all1$COVID_CASE_COUNT)-.00001,all1$COVID_CASE_COUNT),
                         n=8, style="jenks")
all1<-mutate(all1, COVID_CASES_JK=cut(COVID_CASE_COUNT, breaksjk$brks,
                                      dig.lab=6, digits=1))

ggplot(all1)+
  geom_sf(aes(fill=COVID_CASES_JK))+
  geom_sf(data=dtcenters, colour=ptColor, size=1.2) +
  geom_sf_label(data = all1 %>% dplyr::filter(COVID_CASE_COUNT > 9700),
                aes(label = COVID_CASE_COUNT),
                label.size = .09,
                size = 2.5) +
   coord_sf(xlim = c(-74.3, -73.65), ylim=c(40.9,40.5)) +
  scale_fill_brewer(palette="YlOrRd", name="Case Count")+
  labs(x='Longitude', y='Latitude', 
       title='COVID Cases per NYC Zip Code',
       caption = 'Data Source: [2021 COVID Data]')
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
breaksjk2<-classIntervals(c(min(all1$blackPctg)-.00001,all1$blackPctg),
                          n=7, style="jenks")
## Warning in classIntervals(c(min(all1$blackPctg) - 1e-05, all1$blackPctg), : var
## has missing values, omitted in finding classes
all1<-mutate(all1, blackPctgJK=cut(blackPctg,breaksjk2$brks,
                                   dig.lab=4,digits=1))

ggplot(all1)+
  geom_sf(aes(fill=blackPctgJK))+
  geom_sf(data=dtcenters, colour=ptColor, size=1.2) +
  coord_sf(xlim = c(-74.3, -73.65), ylim=c(40.9,40.5)) +
  scale_fill_brewer(palette="YlGn", name="Percentage of Black People")+
  labs(x='Longitude', y='Latitude', 
       title='Percentage of Black People per NYC Zip Code',
       caption = 'Data Source: [2020 ACS Data]')

Interactive Map

Interactive map of COVID cases and Black populations

pal_black<-brewer.pal(7,"YlGnBu")

Map<- all1 %>%
  rename(COVID_Cases=COVID_CASES_JK, Black_Percentage=blackPctgJK)

mapview(Map, zcol=c("COVID_Cases","Black_Percentage"), 
        col.regions=list(pal,pal_black))
## Warning: Found more colors (8) than zcol values (7)! 
## Trimming colors to match number of zcol values.