Packages

require(tidyverse);
## Loading required package: 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
require(sf); 
## Loading required package: sf
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
require(mapview); 
## Loading required package: mapview
require(magrittr); 
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## 
## The following object is masked from 'package:purrr':
## 
##     set_names
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

R Assignment

Reading in the Data

covid<-read_csv("R-Spatial_II_Lab/R-Spatial_II_Lab/tests-by-zcta_2021_04_23.csv", lazy=FALSE)
## Rows: 177 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (3): NEIGHBORHOOD_NAME, BOROUGH_GROUP, label
## dbl (10): MODIFIED_ZCTA, lat, lon, COVID_CASE_COUNT, COVID_CASE_RATE, POP_DE...
## 
## ℹ 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.
zc<-st_read('ZIP_CODE_040114/ZIP_CODE_040114.shp')
## Reading layer `ZIP_CODE_040114' from data source 
##   `C:\Users\OWNER\Downloads\intro to r class\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)
fs<-st_read('R-Spatial_II_Lab/R-Spatial_II_Lab/nycFoodStore.shp')
## Reading layer `nycFoodStore' from data source 
##   `C:\Users\OWNER\Downloads\intro to r class\R-Spatial_II_Lab\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
hf<-read.csv('NYS_Health_Facility.csv')
hf_clean<-hf %>%
  filter(complete.cases(Facility.Longitude, Facility.Latitude))
acs<-read.csv("R-Spatial_II_Lab/R-Spatial_II_Lab/ACSDP5Y2018.DP05_data_with_overlays_2020-04-22T132935.csv")
censust<-st_read('R-Spatial_II_Lab/R-Spatial_II_Lab/2010 Census Tracts/geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a.shp')
## Reading layer `geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a' from data source `C:\Users\OWNER\Downloads\intro to r class\R-Spatial_II_Lab\R-Spatial_II_Lab\2010 Census Tracts\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)

COVID and Zip Codes

Joining covid data and zip code data. Covid data without zip codes will be used later for joining.

covid_zc<-merge(zc,covid, by.x="ZIPCODE", by.y="MODIFIED_ZCTA")
mapview(covid_zc, zcol=c('COVID_CASE_COUNT'))
mapview(covid_zc, zcol=c('TOTAL_COVID_TESTS'))

Aggregate Food Store and Zip Codes

Aggregation of nyc food store data and zip code areas

st_crs(fs)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
st_crs(zc)
## Coordinate Reference System:
##   User input: NAD83 / New York Long Island (ftUS) 
##   wkt:
## PROJCRS["NAD83 / New York Long Island (ftUS)",
##     BASEGEOGCRS["NAD83",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["SPCS83 New York Long Island zone (US survey foot)",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",40.1666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-74,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",41.0333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",40.6666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",984250,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["United States (USA) - New York - counties of Bronx; Kings; Nassau; New York; Queens; Richmond; Suffolk."],
##         BBOX[40.47,-74.26,41.3,-71.8]],
##     ID["EPSG",2263]]
zc84<-st_transform(zc, 4326)
somefs<- fs %>%
  filter(Estbl_T %in% c('JAC','A'))

somefs1<-st_join(zc84, fs, join= st_contains)
  
somefs2<- somefs1 %>%
  group_by(ZIPCODE) %>%
  summarise(PO_NAME=first(PO_NAME),
            COUNTY=first(COUNTY),
            FSNum = n())

mapview(somefs2, zcol=c('FSNum'))

Aggregate Health Facility and Joining to Food Store Data

Aggregation of some nyc health facilities to previous nyc food store zip codes.

hf_sf<-st_as_sf(hf_clean, coords=c('Facility.Longitude','Facility.Latitude'))
st_crs(hf_sf)<-4326

hfsf_dtc<- hf_sf %>%
  filter(Short.Description=='DTC') %>%
  st_join(somefs2, ., join= st_contains)


dtc<- hfsf_dtc %>%
  group_by(ZIPCODE,FSNum) %>%
  summarise(PO_NAME=first(PO_NAME),
            COUNTY=first(COUNTY),
            HFNum = n())
## `summarise()` has grouped output by 'ZIPCODE'. You can override using the
## `.groups` argument.

Joining Covid Data to Health Facilities and Food Stores

Merging the original covid data to the recently produced dtc sf object and removing some unnecessary columns.

dtc1<-merge(dtc,covid,by.x='ZIPCODE',by.y='MODIFIED_ZCTA')
dtc1_clean<- dtc1 %>%
  select(-label,-COVID_CASE_RATE,-COVID_DEATH_COUNT,-COVID_DEATH_RATE)

ACS and Census Tracts

Joining ACS census data and the census tract shapefile data. Renaming the columns in ACS and adding the tractFIPS to census tract data.

acs1<- acs %>%
  select(GEO_ID, 
                totPop = DP05_0001E, elderlyPop = DP05_0024E, 
                malePop = DP05_0002E, femalePop = DP05_0003E,  
                whitePop = DP05_0037E, blackPop = DP05_0038E,
                asianPop = DP05_0067E, hispanicPop = DP05_0071E,
                adultPop = DP05_0021E, citizenAdult = DP05_0087E) %>%
  mutate(censusCode = stringr::str_sub(GEO_ID, -9,-1))

census_fixed<- censust %>%
  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_censust <- merge(census_fixed, acs1, by.x ='tractFIPS', by.y = 'censusCode')

Aggregate ACS and Joining to Rest of the Data

Aggregation of ACS data to the other data that includes covid cases, health facilities, and food stores.

acs_censust84<-st_transform(acs_censust,st_crs(zc84))
st_crs(acs_censust84) 
## 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]]
acs_censust84_centroids <- st_centroid(acs_censust84)


all<-st_join(dtc1_clean,acs_censust84_centroids, join=st_contains)

all1<- all %>%
  group_by(ZIPCODE, PO_NAME, COUNTY, COVID_CASE_COUNT, 
           TOTAL_COVID_TESTS, HFNum, FSNum) %>% 
  summarise(totPop = sum(totPop),
            elderlyPctg = sum(elderlyPop)/totPop*100,
            malePctg = sum(malePop)/totPop*100,
            femalePctg = sum(femalePop)/totPop*100,
            asianPctg = sum(asianPop)/totPop*100,
            blackPctg = sum(blackPop)/totPop*100,
            hispanicPctg = sum(hispanicPop)/totPop*100,
            whitePctg = sum(whitePop)/totPop*100)
## `summarise()` has grouped output by 'ZIPCODE', 'PO_NAME', 'COUNTY',
## 'COVID_CASE_COUNT', 'TOTAL_COVID_TESTS', 'HFNum'. You can override using the
## `.groups` argument.
mapview(all1, zcol=c('HFNum'))
mapview(all1, zcol=c('FSNum'))
mapview(all1, zcol=c('hispanicPctg'))
mapview(all1, zcol=c('COVID_CASE_COUNT'))

Saving Data

save(acs,acs_censust,acs_censust84,acs_censust84_centroids,acs1,all,all1,census_fixed,censust,covid,dtc,
     dtc1,dtc1_clean,fs,hf,hf_clean,hf_sf,hfsf_dtc,somefs,somefs1,somefs2,zc,zc84,
     file="Session_8_data.RData")