Areal interpolation with ACS data removing water and rights-of-way and clipping to filtered zoning area

areal_interpolation_clean.R

if(!require(pacman)){install.packages("pacman");library(pacman)}
p_load(tidyverse, tidycensus, tigris, areal, sf, mapview, geojsonsf, rmapshaper, lwgeom)
options(tigris_use_cache = T); options(tigris_class = "sf")

## Download Portland zoning
zoning <- geojson_sf(readLines("https://opendata.arcgis.com/datasets/6b26f2ccb71d431f9ce8f34fd8ec1558_16.geojson")) %>% 
st_transform(2913) %>% group_by(ZONE) %>% summarize()

## Download water featuers
water <- tigris::area_water("OR", "Multnomah", 2017) %>% summarize() %>% st_transform(2913)

## Get right-of-way here: "https://raw.githubusercontent.com/plnnr/dump/master/multco_row.rds"
multco_row <- readRDS("multco_row.rds")

## OPTIONAL STEP: Only if you want to verify unit counts with actual building footprints
## For faster load time, download footprints here (http://gis-pdx.opendata.arcgis.com/datasets/building-footprints) and 
## import shapefile; geojson file is too huge; alternatively use code below:
# building_footprints <- geojson_sf(readLines("https://opendata.arcgis.com/datasets/1078959662ce49cca059fe2f8930c194_48.geojson")) %>%
#   filter(UNITS_RES > 0,
#          YEAR_BUILT < 2016) %>% ## Look only at residential units built before 2016, so that it's more comparable to 2013-17 ACS data
#   st_centroid(.)

## Download Portland neighborhoods
neighborhoods <- geojson_sf(readLines("https://opendata.arcgis.com/datasets/9f50a605cf4945259b983fa35c993fe9_125.geojson")) %>% 
  st_transform(2913) 

## Filter for Bridgeton example, which is located in an industrial area and right next to the Columbia River
bridgeton <- neighborhoods %>%  
  filter(NAME == "BRIDGETON") %>%
  select(NAME)

## Define "livable" zones--zones that allow residential uses
res_uses <-  c("RF", "R20", "R10", "R7", "R5", "R2.5", "R3", "R2", "R1", "RH", 
               "RX", "RMP", "CR", "CM1", "CM2", "CM3", "CE", "CX", "EX", "CI2", "IR")

## Filter zoning to residential uses and erase water and rights-of-way 
livable_area <- zoning %>% 
  filter(ZONE %in% res_uses) %>%
  rmapshaper::ms_erase(., water) %>%
  rmapshaper::ms_erase(., multco_row)

## Download unit count at block-group level from ACS; keep whole block group area
units_bg_whole <- get_acs(geography = "block group", 
                            variables = "B25001_001", 
                            year = 2017, 
                            state = "OR", 
                            county = "Multnomah",
                            geometry = TRUE) %>%
  st_transform(2913) %>%
  select(GEOID, estimate)

## Clip block groups to the livable area
units_bg_livable <- units_bg_whole %>%
  rmapshaper::ms_clip(., livable_area) %>%
  lwgeom::st_make_valid(.)

## Preview the three layers in mapview
mapview(units_bg_whole%>%filter(GEOID == "410510072021")) + 
  mapview(units_bg_livable%>%filter(GEOID == "410510072021")) + 
  mapview(bridgeton)

## Preview areal weights
aw_preview_weights(bridgeton, tid = NAME, source = units_bg_livable, sid = GEOID, type = "extensive")

## Interpolate number of housing units using entire block group
aw_interpolate(bridgeton, tid = NAME, source = units_bg_whole, sid = "GEOID",
               weight = "total", output = "tibble", extensive = "estimate")

## Interpolate number of housing units using livable area of block group
aw_interpolate(bridgeton, tid = NAME, source = units_bg_livable, sid = "GEOID",
               weight = "total", output = "tibble", extensive = "estimate")