R Week 08 Assignment
Veronica Welsh
3/22/2025
Explanation of the template
Update the title with your information. Make sure to include identification information so that we know it is your submission.
Also update the author name and date accordingly.
Check out the Source Code from the top-right corner
</>Code menu.
In the following R code chunk, load_packages is the code
chunk name. include=FALSE suggests that the code chunk will
run, but the code itself and its outputs will not be included in the
rendered HTML. echo=TRUE in the following code chunk
suggests that the code and results from running the code will be
included in the rendered HTML.
R Spatial Lab Assignment # 2
Don’t use a single chunk for the entire assignment. Break it into multiple chunks.
task 1: Join the COVID-19 data to the NYC zip code area data (sf or sp polygons).
zip_sf <- st_read("/Users/veronicawelsh/Desktop/gtech r/Week_08/ZIP_CODE_040114/ZIP_CODE_040114.shp")## Reading layer `ZIP_CODE_040114' from data source
## `/Users/veronicawelsh/Desktop/gtech r/Week_08/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)covid_data<-read.csv("/Users/veronicawelsh/Desktop/gtech r/Week_08/tests-by-zcta_2021_04_23.csv")
colnames(zip_sf)## [1] "ZIPCODE" "BLDGZIP" "PO_NAME" "POPULATION" "AREA"
## [6] "STATE" "COUNTY" "ST_FIPS" "CTY_FIPS" "URL"
## [11] "SHAPE_AREA" "SHAPE_LEN" "geometry"covid_data$MODIFIED_ZCTA<-as.character(covid_data$MODIFIED_ZCTA)
covid_zip_joined <- zip_sf %>%
left_join(covid_data, by = c("ZIPCODE" = "MODIFIED_ZCTA"))
plot(covid_zip_joined["COVID_CASE_COUNT"], main = "COVID Cases by NYC Zip Code")
task 2: Aggregate the NYC food retails store data (points) to the zip code data, so that we know how many retail stores in each zip code area. Note that not all locations are for food retail. And we need to choose the specific types according to the data.
food_retail<-st_read("/Users/veronicawelsh/Desktop/gtech r/Week_08/nycFoodStore.shp")## Reading layer `nycFoodStore' from data source
## `/Users/veronicawelsh/Desktop/gtech r/Week_08/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 84colnames(food_retail)## [1] "ï__Cnty" "Lcns_Nm" "Oprtn_T" "Estbl_T" "Entty_N" "DBA_Nam"
## [7] "Strt_Nmb" "Stret_Nm" "Add_L_2" "Add_L_3" "City" "State"
## [13] "Zip_Cod" "Sqr_Ftg" "Locatin" "Coords" "geometry"unique(food_retail$Oprtn_T) #So they are all stores?## [1] "Store"unique(food_retail$Estbl_T) #Or just "A" are retail stores? I will assume that's the case.## [1] "JAC" "A" "JACD" "JACDK" "JAD" "JABCHK" "JACHK" "JABC"
## [9] "JAZ" "JABCK" "JACK" "JACDHK" "JABH" "JACH" "JACDE" "JABCH"
## [17] "JABCDH" "JABK" "JABHK" "JABCD" "JACG" "JACDH" "JADHK" "JKA"
## [25] "JADK" "JAB" "JAHK" "JABCDK" "JACZ" "JAK" "JADO" "JDA"st_crs(zip_sf)## 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]]st_crs(food_retail)## 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]]food_retail_transformed <- st_transform(food_retail, crs = st_crs(zip_sf))
food_stores<-food_retail_transformed %>%
filter(str_detect(Estbl_T, '[A]')) %>%
st_join(zip_sf, ., join = st_contains) %>%
group_by(ZIPCODE) %>%
summarise(FoodStoreNum = n())
plot(food_stores["FoodStoreNum"], breaks = "jenks", main = "Number of Food Stores by NYC Zip Code")
task 3: Aggregate the NYC health facilities (points) to the zip code data. Similarly, choose appropriate subtypes such as nursing homes from the facilities.
library(sf)
library(dplyr)
zip_sf <- st_read("/Users/veronicawelsh/Desktop/gtech r/Week_08/ZIP_CODE_040114/ZIP_CODE_040114.shp")## Reading layer `ZIP_CODE_040114' from data source
## `/Users/veronicawelsh/Desktop/gtech r/Week_08/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)healthfac_data <- read.csv("/Users/veronicawelsh/Desktop/gtech r/Week_08/NYS_Health_Facility.csv")
nycNursingHome <- healthfac_data %>%
dplyr::filter(Short.Description == 'NH')
nycNursingHome_noNA <- nycNursingHome %>%
filter(!is.na(Facility.Longitude) & !is.na(Facility.Latitude))
nycNursingHome_sf <- st_as_sf(nycNursingHome_noNA, coords = c("Facility.Longitude", "Facility.Latitude"), crs = 2263)
nursing_home_counts <- st_join(zip_sf[, "ZIPCODE"], nycNursingHome_sf, join = st_contains) %>%
group_by(ZIPCODE) %>%
summarise(NursingHomeNum = n(), .groups = "drop")
nursing_home_counts_no_geom <- st_drop_geometry(nursing_home_counts)
zip_sf <- zip_sf %>%
left_join(nursing_home_counts_no_geom, by = "ZIPCODE")
plot(zip_sf["NursingHomeNum"], main = "Number of Nursing Homes by NYC Zip Code")
task 4: Join the Census ACS population, race, and age data to the NYC Planning Census Tract Data.
task 5: Aggregate the ACS census data to zip code area data.
nycCensus<-sf::st_read("/Users/veronicawelsh/Desktop/gtech r/Week_08/2010 Census Tracts/geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a.shp", stringsAsFactors = FALSE)## Reading layer `geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a' from data source `/Users/veronicawelsh/Desktop/gtech r/Week_08/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)nycCensus %<>% dplyr::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='')
)
acsData <- readLines("/Users/veronicawelsh/Desktop/gtech r/Week_08/ACSDP5Y2018.DP05_data_with_overlays_2020-04-22T132935.csv") %>%
magrittr::extract(-2) %>%
textConnection() %>%
read.csv(header=TRUE, quote= "\"") %>%
dplyr::select(GEO_ID,
totPop = DP05_0001E, elderlyPop = DP05_0024E, # >= 65
malePop = DP05_0002E, femalePop = DP05_0003E,
whitePop = DP05_0037E, blackPop = DP05_0038E,
asianPop = DP05_0067E, hispanicPop = DP05_0071E,
adultPop = DP05_0021E, citizenAdult = DP05_0087E) %>%
dplyr::mutate(censusCode = stringr::str_sub(GEO_ID, -9,-1));
acsData %>%
magrittr::extract(1:10,)## GEO_ID totPop elderlyPop malePop femalePop whitePop blackPop
## 1 1400000US36005000100 7080 51 6503 577 1773 4239
## 2 1400000US36005000200 4542 950 2264 2278 2165 1279
## 3 1400000US36005000400 5634 710 2807 2827 2623 1699
## 4 1400000US36005001600 5917 989 2365 3552 2406 2434
## 5 1400000US36005001900 2765 76 1363 1402 585 1041
## 6 1400000US36005002000 9409 977 4119 5290 3185 4487
## 7 1400000US36005002300 4600 648 2175 2425 479 2122
## 8 1400000US36005002400 172 0 121 51 69 89
## 9 1400000US36005002500 5887 548 2958 2929 903 1344
## 10 1400000US36005002701 2868 243 1259 1609 243 987
## asianPop hispanicPop adultPop citizenAdult censusCode
## 1 130 2329 6909 6100 005000100
## 2 119 3367 3582 2952 005000200
## 3 226 3873 4507 4214 005000400
## 4 68 3603 4416 3851 005001600
## 5 130 1413 2008 1787 005001900
## 6 29 5905 6851 6170 005002000
## 7 27 2674 3498 3056 005002300
## 8 14 0 131 42 005002400
## 9 68 4562 4237 2722 005002500
## 10 0 1985 1848 1412 005002701popData <- merge(nycCensus, acsData, by.x ='tractFIPS', by.y = 'censusCode')
sum(popData$totPop)## [1] 8443713popNYC <- sf::st_transform(popData, st_crs(covid_zip_joined))
popData %>% magrittr::extract('elderlyPop') %>% plot(breaks = 'jenks')
covidPopZipNYC <- sf::st_join(covid_zip_joined,
popNYC %>% sf::st_centroid(),
join = st_contains) %>%
group_by(ZIPCODE, PO_NAME, POPULATION, COUNTY, COVID_CASE_COUNT, TOTAL_COVID_TESTS) %>%
summarise(totPop = sum(totPop),
malePctg = sum(malePop)/totPop*100,
asianPop = sum(asianPop),
blackPop = sum(blackPop),
hispanicPop = sum(hispanicPop),
whitePop = sum(whitePop)) ## Warning: st_centroid assumes attributes are constant over geometries## `summarise()` has grouped output by 'ZIPCODE', 'PO_NAME', 'POPULATION',
## 'COUNTY', 'COVID_CASE_COUNT'. You can override using the `.groups` argument.