NYC AMI

Acquire, clean, and format ACS median household and family income data with geometries

library(tidycensus)
library(tidyverse)
library(tigris)
library(sf)
library(zctaCrosswalk)
library(htmlwidgets)
options(tigris_use_cache = TRUE)

# variables_ACS2022 <- load_variables(2022, "acs5/profile")
year <- 2022 # year of interest

# GEOIDs of NYC counties/boroughs
NYC_county_GEOIDs <- c(36005, 36047, 36061, 36081, 36085)

# FIPS code of NYC counties/boroughs
NYC_county_FIPS <- c(005, 047, 061, 081, 085)

# get ACS median family income (MFI) and median household income (MHI) by county
nyc_county_income <- get_acs(
  survey = "acs5",
  geography = "county",
  state = 36,
  variables = c(mhi = "DP03_0062", mfi = "DP03_0086"),
  year = year,
  geometry = TRUE,
  cb = FALSE # use core TIGER/Line shapefiles for shoreline mapping
) %>% 
  select(-moe) %>%
  pivot_wider(names_from = variable, values_from = estimate) %>%
  filter(GEOID %in% NYC_county_GEOIDs) %>%
  st_transform(4326) %>%
  erase_water(year = year) %>%
  # convert geomcol to multipolygon
  st_collection_extract(type = c("POLYGON", "POINT", "LINESTRING")) %>% 
  st_make_valid()

# get ACS median family income (MFI) and median household income (MHI) by ZCTA
nyc_zcta_income <- get_acs(
  survey = "acs5",
  geography = "zcta",
  variables = c(mhi = "DP03_0062", mfi = "DP03_0086"),
  year = year,
  geometry = TRUE,
  cb = FALSE # use core TIGER/Line shapefiles for shoreline mapping
) %>% 
  filter(GEOID %in% get_zctas_by_county(NYC_county_GEOIDs)) %>%
  select(-moe) %>%
  pivot_wider(names_from = variable, values_from = estimate) %>%
  st_transform(4326) %>%
  erase_water(year = year) %>%
  # convert geomcol to multipolygon
  st_collection_extract(type = c("POLYGON", "POINT", "LINESTRING")) %>% 
  st_make_valid()

# get ACS median family income (MFI) and median household income (MHI) by Tract
nyc_tract_income <- get_acs(
  survey = "acs5",
  geography = "tract",
  state = 36,
  county = NYC_county_FIPS,
  variables = c(mhi = "DP03_0062", mfi = "DP03_0086"),
  year = year,
  geometry = TRUE,
  cb = FALSE # use core TIGER/Line shapefiles for shoreline mapping
) %>% 
  select(-moe) %>%
  pivot_wider(names_from = variable, values_from = estimate) %>%
  st_transform(4326) %>%
  erase_water(year = year) %>%
  # convert geomcol to multipolygon
  st_collection_extract(type = c("POLYGON", "POINT", "LINESTRING")) %>% 
  st_make_valid()

# read in shapefile of NYC 2020 community district tabulations areas (CDTAs)
nyc_cdtas <- read_sf("../data/vector/raw/nycdta2020_24c/nycdta2020.shp") %>% 
  st_transform(4326)

# read in ACS 5-year 2018-2022 economic data by CDTA, provided by NYC DCP
# filter for identifying info and household/family income
nyc_cdta_income <- 
  read_csv("../data/tabular/raw/NYC-DCP_ACS5_1822_CDTA_Econ.csv") %>%
  select(GeoType, CDTAType, GeogName, GeoID, Borough, MdHHIncE, MdFamIncE) %>%
  rename(mhi = MdHHIncE, 
         mfi = MdFamIncE,
         NAME = GeogName)

# join the CDTA geometries to the income data
nyc_cdta_income <- left_join(nyc_cdta_income,
                             nyc_cdtas[, c("CDTA2020", "geometry")], 
                             by = join_by(GeoID == CDTA2020)) %>%
  st_as_sf()

Use income data to create leaflet maps for each scale

library(sf)
library(leaflet)
library(classInt)

# Function for creating a leaflet map of ACS median household income at
# various scales
create_lfmap <- function(data, title, borderWeight) {
  
  # Define a color palette using Fisher-Jenks
  pal <- colorBin(palette = "Greens",
                   bins = classIntervals(data$mhi,
                                         n = 5,
                                         style = "fisher")$brks,
                   pretty = FALSE)
  
  # Hover label text
  labels <- sprintf(
    "<span style='font-size: 20px'><b>%s</b></span><br/>
    2022 Median household income: <strong>$%s</strong><br/>
    MFI - MHI = $%s",
    data$NAME,
    prettyNum(data$mhi, big.mark = ","),
    prettyNum(data$mfi - data$mhi, big.mark = ",")
  ) %>% lapply(htmltools::HTML)
  
  # Create the leaflet map
  leaflet(data = data,
          options = leafletOptions(
            zoomControl = TRUE,
            zoomSnap = 0.25,
            zoomDelta = 1)) %>%
    addProviderTiles("CartoDB.DarkMatter") %>%
    addPolygons(fillColor = ~pal(mhi),
                weight = borderWeight,
                opacity = 1,
                color = "white",
                dashArray = "7",
                fillOpacity = 0.7,
                highlightOptions = highlightOptions(
                  weight = 5,
                  color = "white",
                  dashArray = "",
                  fillOpacity = 0.7,
                  bringToFront = TRUE),
                label = labels,
                labelOptions = labelOptions(
                  style = list("font-weight" = "normal", padding = "3px 8px"),
                  textsize = "15px",
                  direction = "auto")) %>%
    addLegend(pal = pal, 
              values = ~mhi, 
              opacity = 0.7, 
              title = title,
              position = "bottomright")
}

# Create a leaflet map for each scale
nyc_county_map <-
  create_lfmap(nyc_county_income, "NYC Median Household Income by County", 2)
nyc_cdta_map <-
  create_lfmap(nyc_cdta_income, "NYC Median Household Income by CDTA", 1)
nyc_zcta_map <-
  create_lfmap(nyc_zcta_income, "NYC Median Household Income by ZCTA", 1)
nyc_tract_map <-
  create_lfmap(nyc_tract_income, "NYC Median Household Income by Tract", 0.75)

Leaflet map outputs

nyc_county_map

nyc_cdta_map

nyc_zcta_map

nyc_tract_map