load packages
Getting census API key
## Your original .Renviron will be backed up and stored in your R HOME directory if needed.## Your API key has been stored in your .Renviron and can be accessed by Sys.getenv("CENSUS_API_KEY").
## To use now, restart R or run `readRenviron("~/.Renviron")`## [1] "aeb63576bd50b4e3e7ce9beabbac2d8ccc2bf6ce"downloading population data, getting 2020 ACS data for population density
nyc_pop <- get_acs(geography = "tract", variables = "B01003_001", state = "NY", county = c("Bronx", "Kings", "New York", "Queens", "Richmond"), geometry = TRUE, year = 2020)## Getting data from the 2016-2020 5-year ACS## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.## | | | 0% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%loading my download NYC building footprint data
build_footprint <- st_read("/Users/elvagao/Desktop/Rstudio_GTECH/Project_GTECH38520/Building_Footprints/geo_export_1aaeed34-8dab-47fe-a8f7-60738ad1c437.shp") ## Reading layer `geo_export_1aaeed34-8dab-47fe-a8f7-60738ad1c437' from data source `/Users/elvagao/Desktop/Rstudio_GTECH/Project_GTECH38520/Building_Footprints/geo_export_1aaeed34-8dab-47fe-a8f7-60738ad1c437.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1082678 features and 16 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.2555 ymin: 40.49843 xmax: -73.70006 ymax: 40.91505
## Geodetic CRS: WGS 84tract geometry
nyc_tracts <- tracts(state = "NY", county = c("New York", "Kings", "Queens", "Bronx", "Richmond"), cb = TRUE, year = 2020) %>% st_transform(2263)simplify geometry and crop to NYC tracts
build_footprint <- st_transform(build_footprint, st_crs(nyc_pop))fixing the geometries before intersecting (fixed after knowing the error)
build_footprint <- st_make_valid(build_footprint)try again intersecting again
build_footprint <- st_intersection(build_footprint, st_union(nyc_pop))## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries*warning meaning that I am intersecting many building polygons with a single, merged geometry. If a single building is split into many pieces during intersection, the attribute values would stay across those pieces.
count buildings per tract
building_count <- st_join(build_footprint, nyc_pop) %>% group_by(GEOID) %>% summarise(num_buildings = n())joining back to population data
nyc_merged <- nyc_pop %>% left_join(st_drop_geometry(building_count), by = "GEOID") %>% mutate(num_buildings = replace_na(num_buildings, 0))suitable score
nyc_merged <- nyc_merged %>% mutate(pop_density_scaled = scale(estimate), building_density_scaled = scale(num_buildings), suitability = -1 * (pop_density_scaled + building_density_scaled))Define a suitability score
nyc_merged <- nyc_merged %>% mutate(suitability = scale(num_buildings) - (scale(estimate)))
nyc_merged$suitability[is.infinite(nyc_merged$suitability) | is.na(nyc_merged$suitability)] <- 0NA
nyc_merged <- nyc_merged %>% mutate(suitability = replace_na(suitability, 0))
nyc_merged <- nyc_merged[!st_is_empty(nyc_merged), ]summary(nyc_merged$suitability)## V1
## Min. :-7.407291
## 1st Qu.:-0.799213
## Median : 0.179096
## Mean :-0.000835
## 3rd Qu.: 0.831006
## Max. : 5.778914sum(is.na(nyc_merged$suitability))## [1] 0table(is.na(nyc_merged$suitability))##
## FALSE
## 2324nyc_merged$suitability[is.na(nyc_merged$suitability) | is.nan(nyc_merged$suitability)] <- 0create a color palette based on quantiles of suitability
pal <- colorQuantile("RdYlGn", nyc_merged$suitability, n = 5, na.color = "#ccc")leaflet can only be in 4326
nyc_merged_map <- st_transform(nyc_merged, 4326)suitability map (change and upgrade from professor’s advice)
leaflet(data = nyc_merged) %>% addProviderTiles("OpenStreetMap") %>% addPolygons(fillColor = ~pal(suitability), color = "gray", weight = 0.3, opacity = 1, fillOpacity = 0.8, popup = ~paste0("<strong>Census Tract:</strong> ", GEOID, "<br>", "<strong>Suitability Index:</strong> ", round(suitability, 2), "<br>", "<strong>Population:</strong> ", estimate, "<br>", "<strong>Building Count:</strong> ", num_buildings), highlightOptions = highlightOptions(weight = 2, color = "black", bringToFront = TRUE)) %>% addLegend(pal = pal, values = ~suitability, title = "Suitability Index", position = "bottomleft") ## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'populations density in people per square kilometer
nyc_merged <- nyc_merged %>% mutate(area_km2 = as.numeric(st_area(geometry)) / 1e6, pop_density = estimate / area_km2) leaflet only display in 4326
nyc_merged_map <- st_transform(nyc_merged, 4326)create color palette based on quantiles of suitability
pal_density <- colorQuantile("YlOrRd", nyc_merged_map$pop_density, n = 5, na.color = "#ccc")leaflet(data = nyc_merged_map) %>% addProviderTiles("OpenStreetMap") %>% addPolygons(fillColor = ~pal_density(pop_density), color = "white", weight = 0.3, opacity = 1, fillOpacity = 0.8, popup = ~paste0("<strong>Census Tract:</strong> ", GEOID, "<br>", "<strong>Population:</strong> ", estimate, "<br>", "<strong>Pop Density (per km²):</strong> ", round(pop_density, 1), "<br>", "<strong>Building Count:</strong> ", num_buildings), highlightOptions = highlightOptions(weight = 2, color = "black", bringToFront = TRUE)) %>% addLegend(pal = pal_density, values = ~pop_density, title = "Population Density<br>(people/km²)", position = "bottomleft")