Percentage Redlined

pct_redlined.R

require('tigris')
require('tidyr')
require('dplyr')
require('sf')

options(tigris_use_cache=FALSE)

tracts <- tracts("MA", year = 2013) %>%
  st_as_sf() %>%
  mutate(cens_land = as.numeric(ALAND) * 10.76391) %>%
  select(-c(NAME, NAMELSAD, MTFCC, FUNCSTAT, ALAND, AWATER, INTPTLAT, INTPTLON)) %>%
  rename_all(tolower) %>%
  st_transform(2249)

# Create bounding box from Census tracts to narrow query results from CARTO database.
bbox <- paste(st_bbox(tracts), collapse = ",")

# Download only those HOLC areas that are within the rectangle created by the extents of census geographies. 
# Uses a PostGIS query to digitalscholarshiplab CARTO API.
holc.url <- "https://digitalscholarshiplab.carto.com:443/api/v2/sql?format=GeoJSON&q="
holc.query <- paste("SELECT holc_grade, the_geom 
                    FROM digitalscholarshiplab.holc_polygons AS holc 
                    WHERE ST_Intersects(holc.the_geom, ST_Transform(ST_MakeEnvelope(", bbox, ",2249), 4326))", sep='') %>%
  URLencode(.)

# Download and project to Massachusetts Mainland (NAD 1983), EPSG: 2249.
holc <- st_read(paste(holc.url, holc.query, sep = '')) %>%
  st_transform(2249) %>%
  .[tracts,]

# Determine which counties intersect with HOLC zones
# (To constrain number of census hydrology downloads)
counties <- counties("MA", year = 2013) %>%
  st_as_sf() %>%
  st_transform(2249) %>%
  .[holc,] %>%
  rename_all(tolower) %>%
  st_set_geometry(NULL) %>%
  pull(countyfp) %>%
  as.numeric()

# Download water for counties with HOLC districts.
for (c in counties) {
  w <- area_water("MA", county = as.numeric(c), year = 2013) %>%
    st_as_sf()
  if (exists('water')) {
    water <- rbind(water, w)
  } else {
    water <- w 
  }
}

# Select only those census tracts that intersect HOLC zones.
# (To constrain size and complexity of unioned water polygon)
tracts_holc <- tracts[holc,]
# Project water polygon, limit to only those tracts with HOLC zones.
water <- water %>%
  st_transform(2249) %>%
  .[tracts_holc,] %>%
  st_union()

# 'Erase' water polygons.
# (This is by far the most computationally taxing step.)
diff_tracts <- st_difference(tracts_holc, water)

# Calculate total land area of census geographies.
diff_tracts$calc_area <- as.numeric(st_area(diff_tracts$geometry))

# Intersect holc polygons with census tracts.
int <- st_intersection(holc, diff_tracts)

# Calculate area of intersected polgons.
int$holc_area <- st_area(int$geometry)

# Calculate total area in redlined areas by census geographies and remove geometries (can't merge two spatial dataframes).
holc_by_tract <- int %>%
  group_by(geoid, holc_grade) %>%
  summarise(
    holc_area = as.numeric(sum(holc_area)),
    calc_area = as.numeric(max(calc_area))
  ) %>%
  st_set_geometry(NULL) %>%
  spread(holc_grade, holc_area, fill = 0)

# Merge redlined areas with original census geographies.
tracts <- merge(tracts, holc_by_tract, by = "geoid", all.x = TRUE) %>%
  st_transform(4326)

# Export to GeoJSON.
st_write(tracts, dsn = '<path>/<filename>.<format>', delete_dsn = TRUE)