R Week 08 Assignment
Jiada
3/29/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 # 1
Don’t use a single chunk for the entire assignment. Break it into multiple chunks.
You can name the code chunk and also set options.
task 1: joining the COVID-19 data to the NYC zip code area data
nyc_tests <- readr::read_csv("Data/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.str(nyc_tests)## spc_tbl_ [178 × 4] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ MODZCTA : num [1:178] NA 10001 10002 10003 10004 ...
## $ Positive : num [1:178] 1934 211 539 279 23 ...
## $ Total : num [1:178] 2082 448 1024 662 59 ...
## $ zcta_cum.perc_pos: num [1:178] 92.9 47.1 52.6 42.1 39 ...
## - attr(*, "spec")=
## .. cols(
## .. MODZCTA = col_double(),
## .. Positive = col_double(),
## .. Total = col_double(),
## .. zcta_cum.perc_pos = col_double()
## .. )
## - attr(*, "problems")=<externalptr>str(nyc_sf)## Classes 'sf' and 'data.frame': 263 obs. of 13 variables:
## $ ZIPCODE : chr "11436" "11213" "11212" "11225" ...
## $ BLDGZIP : chr "0" "0" "0" "0" ...
## $ PO_NAME : chr "Jamaica" "Brooklyn" "Brooklyn" "Brooklyn" ...
## $ POPULATION: num 18681 62426 83866 56527 72280 ...
## $ AREA : num 22699295 29631004 41972104 23698630 36868799 ...
## $ STATE : chr "NY" "NY" "NY" "NY" ...
## $ COUNTY : chr "Queens" "Kings" "Kings" "Kings" ...
## $ ST_FIPS : chr "36" "36" "36" "36" ...
## $ CTY_FIPS : chr "081" "047" "047" "047" ...
## $ URL : chr "http://www.usps.com/" "http://www.usps.com/" "http://www.usps.com/" "http://www.usps.com/" ...
## $ SHAPE_AREA: num 0 0 0 0 0 0 0 0 0 0 ...
## $ SHAPE_LEN : num 0 0 0 0 0 0 0 0 0 0 ...
## $ geometry :sfc_POLYGON of length 263; first list element: List of 1
## ..$ : num [1:159, 1:2] -73.8 -73.8 -73.8 -73.8 -73.8 ...
## ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:12] "ZIPCODE" "BLDGZIP" "PO_NAME" "POPULATION" ...test_zip_merge <- base::merge(nyc_sf, nyc_tests, by.x = "ZIPCODE", by.y = "MODZCTA")
names(test_zip_merge) ## [1] "ZIPCODE" "BLDGZIP" "PO_NAME"
## [4] "POPULATION" "AREA" "STATE"
## [7] "COUNTY" "ST_FIPS" "CTY_FIPS"
## [10] "URL" "SHAPE_AREA" "SHAPE_LEN"
## [13] "Positive" "Total" "zcta_cum.perc_pos"
## [16] "geometry"task 2: Aggregate NYC ret data to zip code data
install.packages("lwgeom")## Error in install.packages : Updating loaded packageslibrary(lwgeom)
sf_use_s2(FALSE)
nyc_sf %>% mutate(tract_area = st_area(geometry)) %>%
head()## Simple feature collection with 6 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -73.99193 ymin: 40.63029 xmax: -73.78805 ymax: 40.6863
## Geodetic CRS: WGS 84
## ZIPCODE BLDGZIP PO_NAME POPULATION AREA STATE COUNTY ST_FIPS CTY_FIPS
## 1 11436 0 Jamaica 18681 22699295 NY Queens 36 081
## 2 11213 0 Brooklyn 62426 29631004 NY Kings 36 047
## 3 11212 0 Brooklyn 83866 41972104 NY Kings 36 047
## 4 11225 0 Brooklyn 56527 23698630 NY Kings 36 047
## 5 11218 0 Brooklyn 72280 36868799 NY Kings 36 047
## 6 11226 0 Brooklyn 106132 39408598 NY Kings 36 047
## URL SHAPE_AREA SHAPE_LEN geometry
## 1 http://www.usps.com/ 0 0 POLYGON ((-73.80585 40.6829...
## 2 http://www.usps.com/ 0 0 POLYGON ((-73.9374 40.67973...
## 3 http://www.usps.com/ 0 0 POLYGON ((-73.90294 40.6708...
## 4 http://www.usps.com/ 0 0 POLYGON ((-73.95797 40.6706...
## 5 http://www.usps.com/ 0 0 POLYGON ((-73.97208 40.6506...
## 6 http://www.usps.com/ 0 0 POLYGON ((-73.9619 40.65487...
## tract_area
## 1 2108844 [m^2]
## 2 2752823 [m^2]
## 3 3899350 [m^2]
## 4 2201683 [m^2]
## 5 3425228 [m^2]
## 6 3661185 [m^2]#making both datasets the same CRS to avoid error later
st_crs(nys_ret_sf) <- 4326
st_crs(nys_ret_sf)## 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)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## 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(nyc_sf) <- 4326
nyc_sf <- st_transform(nyc_sf, 4326)
st_crs(nyc_sf)## 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)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## 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]]#joining the two datasets
sf_use_s2(FALSE)
nys_ret_sf %>% dplyr::filter(stringr::str_detect(Establishment.Type, '[A]')) %>%
sf::st_join(nyc_sf) %>%
group_by(ZIPCODE) %>%
summarise(FoodStoreNum = n()) %>%
magrittr::extract('FoodStoreNum') %>%
plot(breaks = "jenks", main="Number of Food Stores")## although coordinates are longitude/latitude, st_intersects assumes that they are
## planar
## although coordinates are longitude/latitude, st_intersects assumes that they are
## planar## although coordinates are longitude/latitude, st_union assumes that they are
## planar
##Task 3: Aggregate the NYC health facilities data w/ zipcodes
nyc_sf %>% mutate(tract_area = st_area(geometry)) %>%
head()## Simple feature collection with 6 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -73.99193 ymin: 40.63029 xmax: -73.78805 ymax: 40.6863
## Geodetic CRS: WGS 84
## ZIPCODE BLDGZIP PO_NAME POPULATION AREA STATE COUNTY ST_FIPS CTY_FIPS
## 1 11436 0 Jamaica 18681 22699295 NY Queens 36 081
## 2 11213 0 Brooklyn 62426 29631004 NY Kings 36 047
## 3 11212 0 Brooklyn 83866 41972104 NY Kings 36 047
## 4 11225 0 Brooklyn 56527 23698630 NY Kings 36 047
## 5 11218 0 Brooklyn 72280 36868799 NY Kings 36 047
## 6 11226 0 Brooklyn 106132 39408598 NY Kings 36 047
## URL SHAPE_AREA SHAPE_LEN geometry
## 1 http://www.usps.com/ 0 0 POLYGON ((-73.80585 40.6829...
## 2 http://www.usps.com/ 0 0 POLYGON ((-73.9374 40.67973...
## 3 http://www.usps.com/ 0 0 POLYGON ((-73.90294 40.6708...
## 4 http://www.usps.com/ 0 0 POLYGON ((-73.95797 40.6706...
## 5 http://www.usps.com/ 0 0 POLYGON ((-73.97208 40.6506...
## 6 http://www.usps.com/ 0 0 POLYGON ((-73.9619 40.65487...
## tract_area
## 1 2108844 [m^2]
## 2 2752823 [m^2]
## 3 3899350 [m^2]
## 4 2201683 [m^2]
## 5 3425228 [m^2]
## 6 3661185 [m^2]# filtering the data to only get nursing homes
nycNursingHome <- nys_health_sf %>%
dplyr::filter(Description == 'Residential Health Care Facility - SNF')
nys_health_sf %>% sf::st_join(nyc_sf) %>%
group_by(ZIPCODE) %>%
summarise(NursingHomeNum= n()) %>%
magrittr::extract('NursingHomeNum') %>%
plot(breaks = "jenks", main="Number of Nursing Homes")## although coordinates are longitude/latitude, st_intersects assumes that they are
## planar
## although coordinates are longitude/latitude, st_intersects assumes that they are
## planar## although coordinates are longitude/latitude, st_union assumes that they are
## planar
Task 4: Joining ACS pop, race, age data to census tract data
#removing the first row of info
acsData <- readLines("Data/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));
#making the census shp data into an sf object
nyc_census_tracts_sf <- st_read('Data/2010 Census Tracts/geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a.shp')## Reading layer `geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a' from data source
## `/Users/jiadavalenza/Desktop/spring 2025/GTECH 385/Spatial R/Week 8/Data/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)str(nyc_census_tracts_sf)## Classes 'sf' and 'data.frame': 2165 obs. of 12 variables:
## $ boro_code : chr "5" "1" "1" "1" ...
## $ boro_ct201: chr "5000900" "1009800" "1010000" "1010200" ...
## $ boro_name : chr "Staten Island" "Manhattan" "Manhattan" "Manhattan" ...
## $ cdeligibil: chr "E" "I" "I" "I" ...
## $ ct2010 : chr "000900" "009800" "010000" "010200" ...
## $ ctlabel : chr "9" "98" "100" "102" ...
## $ ntacode : chr "SI22" "MN19" "MN19" "MN17" ...
## $ ntaname : chr "West New Brighton-New Brighton-St. George" "Turtle Bay-East Midtown" "Turtle Bay-East Midtown" "Midtown-Midtown South" ...
## $ puma : chr "3903" "3808" "3808" "3807" ...
## $ shape_area: num 2497010 1906016 1860938 1860993 1864600 ...
## $ shape_leng: num 7729 5534 5692 5688 5693 ...
## $ geometry :sfc_MULTIPOLYGON of length 2165; first list element: List of 1
## ..$ :List of 1
## .. ..$ : num [1:28, 1:2] -74.1 -74.1 -74.1 -74.1 -74.1 ...
## ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:11] "boro_code" "boro_ct201" "boro_name" "cdeligibil" ...#removing leading zeroes
tract_ACS_merge <- base::merge(nyc_census_tracts_sf, acsData, by.x = "boro_ct201", by.y = "censusCode")
names(tract_ACS_merge) ## [1] "boro_ct201" "boro_code" "boro_name" "cdeligibil" "ct2010"
## [6] "ctlabel" "ntacode" "ntaname" "puma" "shape_area"
## [11] "shape_leng" "GEO_ID" "totPop" "elderlyPop" "malePop"
## [16] "femalePop" "whitePop" "blackPop" "asianPop" "hispanicPop"
## [21] "adultPop" "citizenAdult" "geometry"##Task 5: Aggregating ACS Data to Zipcode data
# making sure the CRS matches
st_crs(tract_ACS_merge)## Coordinate Reference System:
## User input: WGS84(DD)
## wkt:
## GEOGCRS["WGS84(DD)",
## DATUM["WGS84",
## ELLIPSOID["WGS84",6378137,298.257223563,
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]]]ACS_2263 <- st_transform(tract_ACS_merge, 2263)
st_crs(ACS_2263)## Coordinate Reference System:
## User input: EPSG:2263
## 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(test_zip_merge)## 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)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## 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]]zip_2263 <- st_transform(test_zip_merge, 2263)
PopZipNYC <- sf::st_join(zip_2263,
ACS_2263 %>% sf::st_centroid(),
join = st_contains) %>%
group_by(ZIPCODE, PO_NAME, POPULATION, COUNTY) %>%
summarise(totPop = sum(totPop),
malePctg = sum(malePop)/totPop*100,
asianPop = sum(asianPop),
blackPop = sum(blackPop),
hispanicPop = sum(hispanicPop),
whitePop = sum(whitePop)) ## `summarise()` has grouped output by 'ZIPCODE', 'PO_NAME', 'POPULATION'. You can
## override using the `.groups` argument.install.packages("knitr")##
## The downloaded binary packages are in
## /var/folders/4x/4fs4vwpx6hl7r_rqmtc_1t0c0000gn/T//RtmpuKfhNz/downloaded_packageslibrary(knitr)
install.packages("rmarkdown")## Error in install.packages : Updating loaded packageslibrary(rmarkdown)