R Week 08 Assignment

R Spatial Lab Assignment # 2

task 1:

nyc_zip_codes <- st_read(
  dsn = "R-Spatial_I_Lab/NYC_COVID_Assignment_Week7.gpkg",
  layer = "nyc_zip_codes",
  quiet = TRUE
)

covid_data <- read_csv("R-Spatial_II_Lab/tests-by-zcta_2020_04_12.csv")

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

names(nyc_zip_codes)

##  [1] "ZIPCODE"    "BLDGZIP"    "PO_NAME"    "POPULATION" "AREA"      
##  [6] "STATE"      "COUNTY"     "ST_FIPS"    "CTY_FIPS"   "URL"       
## [11] "SHAPE_AREA" "SHAPE_LEN"  "geom"

names(covid_data)

## [1] "MODZCTA"           "Positive"          "Total"            
## [4] "zcta_cum.perc_pos"

covid_data <- covid_data %>%
  mutate(MODZCTA = as.character(MODZCTA))

zpNYC <- nyc_zip_codes %>%
  left_join(covid_data, by = c("ZIPCODE" = "MODZCTA"))

task 2:

nycFoodStoreSF <- st_read("R-Spatial_II_Lab/nycFoodStore.shp")

## Reading layer `nycFoodStore' from data source 
##   `/Users/giannaselleck/Documents/RStudio/Week_8/R-Spatial_II_Lab/nycFoodStore.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 11300 features and 16 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.2484 ymin: 40.50782 xmax: -73.67061 ymax: 40.91008
## Geodetic CRS:  WGS 84

food_sel <- nycFoodStoreSF %>%
  filter(str_detect(Estbl_T, "[AJD]"))

zpNYC_food <- st_join(zpNYC, food_sel, join = st_contains) %>%
  group_by(ZIPCODE, PO_NAME, COUNTY, Positive, Total) %>%
  summarise(FoodStoreNum = n(), .groups = "drop")

plot(zpNYC_food["FoodStoreNum"], 
     breaks = "jenks", 
     main = "Number of Selected Food Stores per ZIP")

task 3:

nycHealthFacilities <- read_csv("R-Spatial_II_Lab/NYS_Health_Facility.csv")

## Rows: 3990 Columns: 36
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (28): Facility Name, Short Description, Description, Facility Open Date,...
## dbl  (8): Facility ID, Facility Phone Number, Facility Fax Number, Facility ...
## 
## ℹ 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.

nycHealthFacilities_sf <- nycHealthFacilities %>%
  filter(!is.na(`Facility Longitude`), !is.na(`Facility Latitude`)) %>%
  st_as_sf(coords = c("Facility Longitude", "Facility Latitude"), crs = 4326)

nh_sel <- nycHealthFacilities_sf %>%
  filter(`Short Description` == "NH")

zpNYC_health <- st_join(zpNYC_food, nh_sel, join = st_contains) %>%
  group_by(ZIPCODE, PO_NAME, COUNTY, Positive, Total, FoodStoreNum) %>%
  summarise(NH_Count = n(), .groups = "drop")

task 4:

nycCensus <- st_read("R-Spatial_II_Lab/geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a.shp")

## Reading layer `geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a' from data source `/Users/giannaselleck/Documents/RStudio/Week_8/R-Spatial_II_Lab/geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2165 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91553
## Geodetic CRS:  WGS84(DD)

nycCensus <- nycCensus %>%
  mutate(
    cntyFIPS = case_when(
      boro_name == "Bronx" ~ "005",
      boro_name == "Brooklyn" ~ "047",
      boro_name == "Manhattan" ~ "061",
      boro_name == "Queens" ~ "081",
      boro_name == "Staten Island" ~ "085"
    ),
    tractFIPS = paste(cntyFIPS, ct2010, sep = "")
  )

acs_raw <- readLines("./R-Spatial_II_Lab/ACSDP5Y2018.DP05_data_with_overlays_2020-04-22T132935.csv")

acsData <- acs_raw %>%
  magrittr::extract(-2) %>%
  textConnection() %>%
  read.csv(header = TRUE) %>%
  select(
    GEO_ID,
    totPop = DP05_0001E,
    elderlyPop = DP05_0024E,
    whitePop = DP05_0037E,
    blackPop = DP05_0038E,
    asianPop = DP05_0067E,
    hispanicPop = DP05_0071E
  ) %>%
  mutate(censusCode = stringr::str_sub(GEO_ID, -9, -1))

popData <- merge(nycCensus, acsData,
                 by.x = "tractFIPS",
                 by.y = "censusCode")

task 5:

popNYC <- st_transform(popData, st_crs(zpNYC_health))

pop_points <- st_centroid(popNYC)

## Warning: st_centroid assumes attributes are constant over geometries

pop_zip <- st_join(zpNYC_health, pop_points, join = st_contains) %>%
  group_by(ZIPCODE, PO_NAME, COUNTY, Positive, Total, FoodStoreNum, NH_Count) %>%
  summarise(
    totPop = sum(totPop, na.rm = TRUE),
    elderlyPop = sum(elderlyPop, na.rm = TRUE),
    whitePop = sum(whitePop, na.rm = TRUE),
    blackPop = sum(blackPop, na.rm = TRUE),
    asianPop = sum(asianPop, na.rm = TRUE),
    hispanicPop = sum(hispanicPop, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  rename(
    COVID_CASE_COUNT = Positive,
    TOTAL_COVID_TESTS = Total
  )

final_sf <- pop_zip
names(final_sf)

##  [1] "ZIPCODE"           "PO_NAME"           "COUNTY"           
##  [4] "COVID_CASE_COUNT"  "TOTAL_COVID_TESTS" "FoodStoreNum"     
##  [7] "NH_Count"          "totPop"            "elderlyPop"       
## [10] "whitePop"          "blackPop"          "asianPop"         
## [13] "hispanicPop"       "geom"

library(classInt)
library(sf)

breaks <- classIntervals(final_sf$totPop, n = 5, style = "jenks")

library(RColorBrewer)
colors <- brewer.pal(5, "YlOrRd")

plot(final_sf["totPop"], breaks = breaks$brks, pal = colors,
     main = "Total Population per ZIP")