Reading layer `ZIP_CODE_040114' from data source
`/Users/ayeshanaveed/Downloads/R-Spatial_I_Lab/ZIP_CODE_040114 2/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)
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.
Rows: 29389 Columns: 18
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (11): ï..County, Operation.Type, Establishment.Type, Entity.Name, DBA.Na...
dbl (4): License.Number, Zip.Code, Y, X
num (1): Square.Footage
lgl (2): Address.Line.2, Address.Line.3
ℹ 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.
# Get the missing values in X and Y columnssum(is.na(retail_food_stores$X))
[1] 5417
sum(is.na(retail_food_stores$Y))
[1] 5417
# Filter out the rows with NA in X and Yretail_food_stores <- retail_food_stores[!is.na(retail_food_stores$X) &!is.na(retail_food_stores$Y), ]retail_food_stores_sf <-st_as_sf(retail_food_stores, coords =c("X", "Y"), crs =4326)retail_food_filtered <- retail_food_stores_sf %>%filter(Establishment.Type %in%c("A(3)", "B"))retail_food_filtered <-st_transform(retail_food_filtered, st_crs(nyc_postal_areas))st_crs(nyc_postal_areas) ==st_crs(retail_food_filtered)
[1] "Hospice"
[2] "Residential Health Care Facility - SNF"
[3] "Diagnostic and Treatment Center"
[4] "Certified Home Health Agency"
[5] "Hospital Extension Clinic"
[6] "School Based Diagnostic and Treatment Center Extension Clinic"
[7] "Diagnostic and Treatment Center Extension Clinic"
[8] "School Based Hospital Extension Clinic"
[9] "Hospital"
[10] "Primary Care Hospital - Critical Access Hospital Extension Clinic"
[11] "Adult Day Health Care Program - Offsite"
[12] "Mobile Diagnostic and Treatment Center Extension Clinic"
[13] "Long Term Home Health Care Program"
[14] "Mobile Hospital Extension Clinic"
[15] "Primary Care Hospital - Critical Access Hospital"
[16] "School Based Primary Care Hospital - Critical Access Extension Clinic"
health_facilities_filtered <- health_facilities_sf %>%filter(Description =="Residential Health Care Facility - SNF")health_facilities_filtered <-st_transform(health_facilities_filtered, st_crs(nyc_postal_areas))health_facility_counts <-st_join(nyc_postal_areas, health_facilities_filtered) %>%group_by(ZIPCODE) %>%summarize(health_facility_count =n())# Make sure the merge was successfulhead(health_facility_counts)
Rows: 2168 Columns: 358
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (358): GEO_ID, NAME, DP05_0031PM, DP05_0032E, DP05_0032M, DP05_0032PE, D...
ℹ 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.
nyc_census_tracts <-st_read("/Users/ayeshanaveed/Downloads/R-Spatial_II_Lab_2/2010 Census Tracts/geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a.shp")
Reading layer `geo_export_1dc7b645-647b-4806-b9a0-7b79660f120a' from data source `/Users/ayeshanaveed/Downloads/R-Spatial_II_Lab_2/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)
# Spatial join: Add ZIP Code information to Census Tract datanyc_postal_areas <-st_transform(nyc_postal_areas, st_crs(popData))st_crs(nyc_postal_areas) ==st_crs(popData) # Should return TRUE
# Convert columns to numericcols_to_convert <-c("DP05_0001E", "DP05_0024E", "DP05_0002E", "DP05_0003E", "DP05_0037E", "DP05_0038E", "DP05_0070E", "DP05_0044E")popData_with_zip <- popData_with_zip %>%mutate(across(all_of(cols_to_convert), # Apply to these columns~ readr::parse_number(.x) # Clean and convert to numeric ))popData_with_zip <- popData_with_zip %>%mutate(DP05_0004E = readr::parse_number(DP05_0004E) # Clean and convert to numeric again for this one )
Warning: There was 1 warning in `stopifnot()`.
ℹ In argument: `DP05_0004E = readr::parse_number(DP05_0004E)`.
Caused by warning:
! 147 parsing failures.
row col expected actual
113 -- a number -
165 -- a number -
166 -- a number -
167 -- a number -
177 -- a number -
... ... ........ ......
See problems(...) for more details.
# Ensure both objects have the same CRSretail_store_counts <-st_transform(retail_store_counts, st_crs(zip_aggregated_data))st_crs(zip_aggregated_data) ==st_crs(retail_store_counts)
[1] TRUE
# Spatial join between zip_aggregated_data and retail_store_countsfinal_sf_object <-st_join(zip_aggregated_data, retail_store_counts, join = st_intersects)# Ensure objects have the same CRS and then perform the spatial joinhealth_facility_counts <-st_transform(health_facility_counts, st_crs(final_sf_object))final_sf_object <-st_join(final_sf_object, health_facility_counts, join = st_intersects)# Merge COVID-19 data covid_data <- covid_data %>%mutate(ZIPCODE =as.character(ZIPCODE))final_sf_object <- final_sf_object %>%left_join(covid_data, by ="ZIPCODE")final_sf_object <- final_sf_object %>%group_by(ZIPCODE) %>%# Group by ZIP Codesummarise(totPop =first(totPop), elderlyPop =first(elderlyPop),malePop =first(malePop),femalePop =first(femalePop),whitePop =first(whitePop),blackPop =first(blackPop),hispanicPop =first(hispanicPop),asianPop =first(asianPop),retail_store_count =sum(retail_store_count, na.rm =TRUE), health_facility_count =sum(health_facility_count, na.rm =TRUE), positive =sum(positive, na.rm =TRUE), total =sum(total, na.rm =TRUE), zcta_cum_perc_pos =mean(zcta_cum_perc_pos, na.rm =TRUE), geometry =st_union(geometry) )# Verify aggregation glimpse(final_sf_object)
# Save as GeoJSONst_write(final_sf_object, "final_zipcode_data.geojson", append=TRUE)
Updating layer `final_zipcode_data' to data source `final_zipcode_data.geojson' using driver `GeoJSON'
Updating existing layer final_zipcode_data
Writing 249 features with 14 fields and geometry type Unknown (any).