Pratt Center Test Analysis
Alaina Turnquist
Project Description
The goal of this project is to create test maps demonstrating spatial visualization skills.
There are 3 sections of maps for Brooklyn Community Board 8:
- Basic and Population Density
- Rent Burden
- Evictions
Set Up
- Add Libraries
# Set the proxy
Sys.setenv(http_proxy="bcpxy.nycnet:8080")
Sys.setenv(https_proxy="bcpxy.nycnet:8080")
## Add libraries
library(tidyverse)
library(leaflet)
library(tidycensus)
library(ggplot2)
library(terra)
library(tidygeocoder)
library(scales)
library(ggspatial)
library(prettymapr)
library(ggmap)
library(lubridate)
library(RColorBrewer)
setwd("C:/Users/aturnquist/OneDrive - NYC O365 HOSTED/Y-Drive/Pratt Center Work")- Import Community District Tabulation Areas and their census tracts. Import Evictions data
# import DCP Community District Tabulation Areas (CDTAs and JIAs)
cdta_polygons <- st_read("C:/Users/aturnquist/OneDrive - NYC O365 HOSTED/Y-Drive/Pratt Center Work/data/nycdta2020_25c/nycdta2020_25c/nycdta2020.shp", quiet = TRUE)
# import list of census 2020 tracts within each CDTA
bk_cb8_cdta_tracts <- st_read("C:/Users/aturnquist/OneDrive - NYC O365 HOSTED/Y-Drive/Pratt Center Work/data/2020_Census_Tracts_to_2020_NTAs_and_CDTAs_Equivalency_20250917.csv", quiet = TRUE)
# import evictions data
evictions <- st_read("C:/Users/aturnquist/OneDrive - NYC O365 HOSTED/Y-Drive/Pratt Center Work/data/Evictions_20250917.csv", quiet = TRUE)- Import Census Data
## Set the variables wanted from ACS
variables = c(vars <- c(
total_population = "B01003_001",
median_income = "B19013_001",
total_units = "B25002_001",
occupied_units = "B25002_002",
vacant_units = "B25002_003",
median_rent = "B25064_001",
median_value = "B25077_001",
rent_30_34 = "B25070_007",
rent_35_39 = "B25070_008",
rent_40_49 = "B25070_009",
rent_50_plus = "B25070_010"
))
## Get the tracts from ACS 2023 5-year and remove progress bars from the knitted html
nyc_tracts <- suppressMessages(
get_acs(
geography = "tract",
variables = vars,
state = "NY",
county = c("Bronx", "Kings", "New York", "Queens", "Richmond"),
year = 2023,
geometry = TRUE,
cache_table= TRUE
))## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |====================== | 31% | |======================= | 33% | |======================== | 34% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%print("NOTE: Progress bars are inevitable here (although very ugly), since the output of TidyCensus prints progress directly in the console. For a cleaner look on the page, I would typically do the download in a separate file and then read in the file for the Census data. However, I wanted to show my work on how I did the pull" ) ## [1] "NOTE: Progress bars are inevitable here (although very ugly), since the output of TidyCensus prints progress directly in the console. For a cleaner look on the page, I would typically do the download in a separate file and then read in the file for the Census data. However, I wanted to show my work on how I did the pull"Manipulate Census Data
Manipulate the Census Data to make it spatial and showing variables by tract
## Pivot so each observation contains the ACS variables for that tract
nyc_tracts <- nyc_tracts %>%
select(GEOID, NAME, variable, estimate,geometry) %>%
pivot_wider(
names_from = variable,
values_from = estimate
)
# complete a left_join to add the cdta names/codes and nta names to nyc_tracts
nyc_tracts <- nyc_tracts %>%
left_join(
bk_cb8_cdta_tracts %>% select(GEOID, CDTAName, CDTACode, NTAName),
by = "GEOID"
)
# drop everything that is not in BK08
BK_CB8_tracts <- nyc_tracts %>%
filter(CDTACode == "BK08")
# create a BK08 border polygon
BK08_border <- cdta_polygons %>%
filter(CDTA2020 == "BK08")Set Coordinate System
## Make sure all layers are in the same coordinate system (WGS 84)
BK_CB8_tracts <- st_transform(BK_CB8_tracts, crs = 4326)
BK08_border <- st_transform(BK08_border, crs = 4326)Create Initial Maps
- Basic District Map Static Version
## create alternative projections
## bk_cb8_tracts_proj <- st_transform(BK_CB8_tracts, 2263)
## cb_bk8_outline_proj <- st_transform(BK08_border, 2263)
## Make the map
CB_Boundary_Map <- ggplot() +
annotation_map_tile(type = "cartolight", zoomin = -1) + # OSM basemap
geom_sf(data = BK08_border, fill = "blue", color = "blue",
alpha = 0.1, size = 0.5) +
theme_minimal() +
labs(title = "Brooklyn Community Board 8 Boundary") +
theme_minimal() +
theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank(),
panel.grid = element_blank()
)
# let's see the map
CB_Boundary_Map
# Now name and save the map to the proper file
ggsave(
filename = "C:/Users/aturnquist/OneDrive - NYC O365 HOSTED/Y-Drive/Pratt Center Work/Maps/cb8_map.png",
plot = CB_Boundary_Map,
width = 8,
height = 6,
dpi = 300
)- Basic District Map Interactive Version
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
addPolygons(data = BK08_border ,
fillColor = "blue",
color = "blue",
weight = 1,
opacity = 1,
fillOpacity = 0.1) %>%
addControl(
html = "<h3 style='text-align:center;'>Brooklyn Community Board 8 Boundary</h3>",
position = "topleft") - District Map by Tract Population Density (static version)
## Calculate km2 area of each tract
### first change it into Ny Plane before doing this transformation
BK_CB8_tracts <- st_transform(BK_CB8_tracts, 2263)
# Calculate area in m2 and then convert to km2, then create a population density variable
BK_CB8_tracts <- BK_CB8_tracts %>%
mutate(
area_m2 = st_area(geometry),
area_km2 = as.numeric(area_m2) / 1e6,
pop_dens_km2 = total_population / area_km2
)# transform it back into WGS 84 coordinate system since we are done with the calculation
BK_CB8_tracts <- st_transform(BK_CB8_tracts, 4326) Now we need to remove Lincoln Terrace/Arthur Somers Park from the data frame, so it will not affect the coloring of the map.
bk8_nopark <- BK_CB8_tracts %>%
filter(total_population > 0)
ggplot() +
annotation_map_tile(type = "cartolight") +
geom_sf(
data = bk8_nopark,
aes(fill = pop_dens_km2),
color = "white",
size = 0.2
) +
geom_sf(
data = BK08_border,
fill = NA,
color = "blue",
size = 1.2
) +
scale_fill_viridis_c(
option = "plasma",
trans = "sqrt",
name = "Pop density (people/km²)",
labels = label_number(big.mark = ",", accuracy = 1) # <- inside scale
) +
labs(title = "Brooklyn CB8 Population Density by Tract") +
theme_minimal() +
theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank(),
panel.grid = element_blank()
)
- District Map by Tract Population Density Interactive Version
## set the color pallette
pal <- colorNumeric(
palette = "plasma",
domain = BK_CB8_tracts$pop_dens_km2[BK_CB8_tracts$pop_dens_km2 > 0],
na.color = "grey"
)
# Create the map
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
addPolygons(
data = bk8_nopark,
fillColor = ~pal(pop_dens_km2),
color = "white",
weight = 0.5,
opacity = 1,
fillOpacity = 0.7,
highlightOptions = highlightOptions(
weight = 2,
color = "black",
fillOpacity = 0.9,
bringToFront = TRUE
),
label = ~paste0("Population density: ", round(pop_dens_km2, 1), " people/km²")
) %>%
addPolygons(
data = BK08_border,
fill = FALSE,
color = "blue",
weight = 3
) %>%
addLegend(
pal = pal,
values = bk8_nopark$pop_dens_km2,
title = "Pop density (people/km²)",
position = "bottomright"
) %>%
addControl(
html = "<h3>Brooklyn Community Board 8 Population Density</h3>",
position = "topright"
)Rent Burden
Let’s look at where renters are paying more than 30% of their income in gross rent by Census Tract.
## we need to start by summing all of the units within each census tract that are above 30%
BK_CB8_tracts <- BK_CB8_tracts %>%
mutate(
rent_burdened_units = rent_30_34 + rent_35_39 + rent_40_49 + rent_50_plus
)
## Now let's create a % of rent burdened units variable
BK_CB8_tracts <- BK_CB8_tracts %>%
mutate(
rent_burden_percent = (rent_burdened_units/occupied_units)*100
)pal <- colorNumeric(
palette = "YlOrRd",
domain = BK_CB8_tracts$rent_burden_percent
)
# Clean popup content
BK_CB8_tracts <- BK_CB8_tracts %>%
mutate(popup_text = paste0(
"<strong>Tract: </strong>", GEOID, "<br>",
"<strong>Rent-burdened: </strong>", round(rent_burden_percent, 1), "%"
))
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>% # Light background
addPolygons(
data = BK_CB8_tracts,
fillColor = ~pal(rent_burden_percent),
color = "white",
weight = 0.5,
opacity = 1,
fillOpacity = 0.7,
label = ~paste0(round(rent_burden_percent, 1), "%"
),
highlightOptions = highlightOptions(
weight = 2,
color = "black",
fillOpacity = 0.9,
bringToFront = TRUE
)
) %>%
addPolygons(
data = BK08_border,
fill = FALSE,
color = "blue",
weight = 2
) %>%
addLegend(
pal = pal,
values = BK_CB8_tracts$rent_burden_percent,
title = "% Rent-Burdened Units",
position = "bottomright"
)Evictions Data Manipulation
## first filter for the CB
evictions_cb8 <- evictions %>%
filter(Community.Board == "8", BOROUGH == "BROOKLYN", Residential.Commercial == "Residential")
## first make lat and long numeric to check for invalid coordinates that will cause problems when mapping
evictions_cb8 <- evictions_cb8 %>%
mutate(
Longitude = as.numeric(Longitude),
Latitude = as.numeric(Latitude)
)
## check for missing values in lat and long that might cause issues when mapping
sum(is.na(evictions_cb8$Longitude))## [1] 0sum(is.na(evictions_cb8$Latitude))## [1] 0## tell r which column are lat and long and then set the coordinate system
evictions_cb8 <- evictions_cb8 %>%
st_as_sf(
coords = c("Longitude", "Latitude"),
crs = 4326,
remove = FALSE
)## Parse the date and create a year column
evictions_cb8 <- evictions_cb8 %>%
mutate(
Executed.Date = mdy(Executed.Date),
Year = year(Executed.Date))
## remove everything that is na in the year column
evictions_by_year <- evictions_cb8 %>%
group_by(Year) %>%
summarize(num_evictions = n()) %>%
arrange(Year)
evictions_by_year## Simple feature collection with 9 features and 2 fields
## Geometry type: MULTIPOINT
## Dimension: XY
## Bounding box: xmin: -73.97429 ymin: 40.66739 xmax: -73.92192 ymax: 40.68212
## Geodetic CRS: WGS 84
## # A tibble: 9 × 3
## Year num_evictions geometry
## <dbl> <int> <MULTIPOINT [°]>
## 1 2017 336 ((-73.97273 40.67847), (-73.97072 40.67805), (-73.9707 40…
## 2 2018 335 ((-73.97283 40.68122), (-73.97255 40.67939), (-73.97117 4…
## 3 2019 264 ((-73.97086 40.67713), (-73.96805 40.67412), (-73.96775 4…
## 4 2020 44 ((-73.97088 40.68145), (-73.96479 40.67869), (-73.96192 4…
## 5 2021 2 ((-73.93573 40.66999), (-73.92634 40.67124))
## 6 2022 83 ((-73.97232 40.68029), (-73.96417 40.67989), (-73.96367 4…
## 7 2023 210 ((-73.97429 40.68131), (-73.97326 40.68212), (-73.96672 4…
## 8 2024 196 ((-73.9665 40.67566), (-73.96567 40.67415), (-73.96461 40…
## 9 2025 188 ((-73.97056 40.67801), (-73.97031 40.67703), (-73.96957 4…Map Residential Evictions in Community Board 8
## Create eviction counts data frame counting evictions by coordinates and address
eviction_counts <- evictions_cb8 %>%
group_by(Longitude, Latitude, Eviction.Address) %>%
summarize(num_evictions = n()) %>%
ungroup()## Reaffirm the coordinate system we want for the border polygon
BK08_border <- st_transform(BK08_border, 4326)
# Convert eviction points to sf for mapping,
eviction_counts_sf <- eviction_counts %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, remove = FALSE)- Map of Evictions from 2017-2025 in Brooklyn Community Board 8 (static)
# Create the map
ggplot() +
annotation_map_tile(type = "cartolight") +
geom_sf(data = BK08_border, fill = NA, color = "purple", size = 1) +
geom_sf(
data = eviction_counts_sf,
aes(size = num_evictions),
color = "orange",
alpha = 0.4
) +
scale_size_continuous(range = c(1, 6)) +
theme_minimal() +
labs(
title = "Residential Evictions in Brooklyn CB8 2017-2025",
size = "Number of evictions"
) +
theme_minimal() +
theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank(),
panel.grid = element_blank()
)
## Here we can see that the most evictions since 2017 are happening in Eastern Crown Heights - Map of Evictions for 2025 to date in Community Board 8 (Static)
## Now let's look at the same thing but just for this year so far in 2025
evictions_2025 <- evictions_cb8 %>%
filter(Year == 2025)
## count evictions
eviction_counts_2025 <- evictions_2025 %>%
group_by(Longitude, Latitude, Eviction.Address) %>%
summarize(num_evictions = n()) %>%
ungroup()## Create the map
ggplot() +
annotation_map_tile(type = "cartolight") +
geom_sf(data = BK08_border, fill = NA, color = "purple", size = 1.2) +
geom_sf(
data = eviction_counts_2025,
aes(size = num_evictions),
color = "orange",
alpha = 0.4
) +
scale_size_continuous(range = c(1, 6)) +
labs(
title = "Residential Evictions in Brooklyn CB8 (2025)",
size = "Number of evictions"
) +
theme_minimal() +
theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank(),
panel.grid = element_blank()
)
- Map of 2025 Evictions to date in Brooklyn Community Board 8 (Interactive)
pal <- colorNumeric("Reds", domain = eviction_counts_2025$num_evictions)
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>% # Basemap
addPolygons(
data = BK08_border,
fill = FALSE,
color = "blue",
weight = 2
) %>%
addCircleMarkers(
data = eviction_counts_2025,
radius = ~sqrt(num_evictions) * 2, # scale size by eviction count
color = "red",
fillOpacity = 0.4,
stroke = FALSE,
label = ~paste0(Eviction.Address, ": ", num_evictions, " evictions")
) %>%
addLegend(
position = "bottomright",
title = "Evictions per Address",
pal = pal,
values = eviction_counts_2025$num_evictions
) %>%
addControl(
"Brooklyn Community Board 8 Evictions (2025)",
position = "topright"
)