streetuse.R

# Reproduction ov https://twitter.com/TransAlt/status/1547237573726191620 for
# Boston, using Census data + Boston GIS vectors

library(dplyr)
library(ggplot2)
library(units)
library(sf)
library(tidycensus)

census_api_key(API KEY HERE)

# Retrieve / load geographies ---------------------------------------------

download_and_unzip <- function(url, temp_file = "temp.zip", ...) {
  download.file(url, "temp.zip")
  unzip("temp.zip", ...)
  file.remove("temp.zip")
}

# Census-designated places
download_and_unzip("https://www2.census.gov/geo/tiger/TIGER2020/PLACE/tl_2020_25_place.zip")
places <- st_read("tl_2020_25_place.shp")

# Tracts
download_and_unzip("https://www2.census.gov/geo/tiger/TIGER2020/TRACT/tl_2020_25_tract.zip")
tracts <- st_read("tl_2020_25_tract.shp")

# Sidewalks and street edges downloaded from the city of Boston's ESRI feature
# server from the following layers:
#
# * sidewalk_inventory.geojson: https://gis.cityofboston.gov/arcgis/rest/services/Infrastructure/OpenData/MapServer/0
# * street_edges.geojson: http://gis.cityofboston.gov/arcgis/rest/services/Basemaps/base_map_webmercatorV2/MapServer/17
sidewalks <- st_read("sidewalk_inventory.geojson")
street_edges <- st_read("street_edges.geojson")

# Block groups within boston ----------------------------------------------
# Using simple definition: centroids within the Boston Census-designated place

boston_tracts <- subset(
  tracts,
  st_within(
    st_centroid(geometry),
    subset(places, NAME == "Boston"),
    sparse = FALSE
  )
)

# Proportion of workers in major transit modes ----------------------------

# Census data table definitions:
# https://www.socialexplorer.com/data/ACS2020_5yr/metadata/?ds=ACS20_5yr&table=B08006
acs5 <- get_acs(
  geography = "tract",
  variables = c(
    
    # Workers 16 years and over
    workers = "B08006_001",
    
    # Workers 16 years and over -> Car, Truck, or Van
    automobile = "B08006_002",
    
    # Workers 16 years and over -> Public Transporation (Excluding Taxicab)
    transit = "B08006_008",
    
    # Workers 16 years and over -> Bicycle
    bicycle = "B08006_014",
    
    # Workers 16 years and over -> Walked
    walked = "B08006_015"
  ),
  state = "MA",
  couty = "Suffolk"
)

acs5_boston <- subset(acs5, GEOID %in% boston_tracts$GEOID)

total_workers = sum(subset(acs5_boston, variable == "workers")$estimate)

commuting_table <- acs5_boston %>%
  group_by(variable) %>%
  summarize(total = sum(estimate)) %>%
  mutate(percentage = total / total_workers)

commuting_table

## Plot #1 ----

commuting_table %>%
  filter(variable != "workers") %>%
  ggplot() +
    aes(x = variable, y = total) +
    geom_bar(stat = "identity") +
    labs(
      title = "Local commuters by commuting type",
      x = "Commuting type",
      y = "Commuters"
    ) +
    theme_minimal()

# Proportion of land for streets vs sidewalks -----------------------------

land_areas <- data.frame(
  category = c("sidewalks", "roads"),
  area_m2 = as.numeric(c(
    set_units(sum(st_area(sidewalks)), "m^2"),
    set_units(sum(st_area(street_edges)), "m^2")
  ))
)

land_areas

## Plot #2 ----

ggplot(land_areas) +
  aes(x = category, y = area_m2) +
  geom_bar(stat = "identity") +
  labs(
    title = "Land use by category",
    x = "Category",
    y = expression("Area (m²)")
  ) +
  theme_minimal()