Ascribing Neighborhood dividedness

This markdown shows the workflow and functions for allocating neighborhood (nbhd) dividedness.

This package bundles a couple functions to query nbhds within a region, and allocate each to a sub-area bounded by a type of division. For example, localities or school districts within a metro area, sides of a interstate/highway relative to one another, or other linear divisions that might be considered.

Sample call

CZ

The function can be called like:

# I use my helper library geox (see below for installation)

# define a place to get nbhd divisions for 
czsf <- geox::build.CZs('20100') # Portland, ME commuting zone

# get divisions that can be overlaid 
hwys <- geox::get.NHPN(sfx = st_transform(czsf, 4326))
## Reading layer `National_Highway_Planning_Network' from data source 
##   `C:\Users\kiram\Dropbox (Princeton)\shapefiles\nhpn\National_Highway_Planning_Network.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2055 features and 46 fields (with 10 geometries empty)
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -71.14064 ymin: 43.53229 xmax: -69.28706 ymax: 45.80578
## Geodetic CRS:  WGS 84
devtools::load_all()

# allocate census tracts to sub-area divisions, using just interstates
nb.div <- divM::gen.cross.tract.dividedness(
  region = czsf
  ,divs = filter(hwys, signt1 == 'I')
  ,nbd.query.fcn = tigris::tracts
  ,fill.nhpn.gaps = T
  ,region.id.colm = 'cz'
  ,erase.water = F
  ,year = 2019
  )

# this object represents neighborhood dividedness in the region.
nb.div
## # A tibble: 185 x 3
##    geoid       poly.id perc.in.div
##    <chr>       <fct>         <dbl>
##  1 23001044000 1             1.00 
##  2 23001046000 1             1.00 
##  3 23001046500 1             0.578
##  4 23001010100 1             1.00 
##  5 23001010300 1             1.00 
##  6 23001010400 1             1.00 
##  7 23001010500 1             1.00 
##  8 23001020100 1             1.00 
##  9 23001020200 1             1.00 
## 10 23001020300 1             1.00 
## # ... with 175 more rows

The nb.div object, which represents neighborhood dividedness in the region has a geoid column for tract or block group identifier, a poly.id column for the sub-polygon division, and a perc.in.div column for the percentage of the neighborhood within the identified sub-polygon division.

We can visualize the neighborhood allocation like:

nb.divsf <- nb.div %>% geox::attach.geos(query.fcn = tigris::tracts, year = 2019)

ggplot() +
  geom_sf(data = nb.divsf
          ,aes(fill = poly.id)
          ,color = 'white', size = .6
          ,alpha = .5) +
  geom_sf(data =  filter(hwys, signt1 == 'I')
          ,aes(color = signt1)
          ,size = 1.3) +
  #scale_fill_discrete(guide= 'none') +
  theme_void()

The pre-generated cross-tract or cross-blockgroup dividedness datasets are generated using Della scripts that just iterate through different regions.

The example, and the function itself, lean on geox, a helper library I wrote that has a lot of convenience functions working with census spatial data, querying through the tigris library, and referencing CZs and CBSAs.

To install this library, use: devtools::install_github("https://github.com/kmcd39/geox.git") Github link: https://github.com/kmcd39/geox/tree/main/R

City-level variation

Variations are easy to generate, with different regions or divisions objects. For example, here’s a separate generation block-group dividedness within the city of Portland (rather than tracts within the commuting zone). We can also use all highways instead of just interstates

# get all places in the Portland CZ
plcs <- geox::places.wrapper(x = czsf)

# filter to just portland
smplc <- plcs %>%
  filter(grepl('^Portland', name)) %>%
  st_transform(4326)

# allocate block groups to sub-area divisions, using all hwys
nb.plc.div <- divM::gen.cross.tract.dividedness(
  region = smplc
  ,divs = hwys
  ,nbd.query.fcn = tigris::block_groups
  ,fill.nhpn.gaps = T
  ,region.id.colm = 'placefp'
  ,erase.water = F
  ,year = 2019
  )

nb.plc.div
## # A tibble: 54 x 3
##    geoid        poly.id perc.in.div
##    <chr>        <fct>         <dbl>
##  1 230050024003 1             0.893
##  2 230050021022 2             0.515
##  3 230050020012 6             0.907
##  4 230050020022 6             0.418
##  5 230050015002 6             0.438
##  6 230050017002 6             0.993
##  7 230050017003 6             0.925
##  8 230050021013 7             0.912
##  9 230050021014 7             0.994
## 10 230050021021 7             0.977
## # ... with 44 more rows
# get background tiles for vis this time
sttm <- visaux::get.stamen.bkg(smplc, zoom = 12)

nb.plc.divsf <- nb.plc.div %>% geox::attach.geos(query.fcn = tigris::block_groups)

ggmap(sttm) +
  geom_sf(
    data = nb.plc.divsf
    ,aes(fill = poly.id)
    ,color = 'white'
    ,alpha = .6
    ,inherit.aes = F
  ) +
  geom_sf(data = hwys
          ,aes(color = factor(signt1))
          ,size =  1.3
          ,inherit.aes = F
  ) +
  scale_fill_discrete(guide = 'none') +
  visaux::bbox2ggcrop(nb.plc.divsf) +
  theme_void() +
  theme(legend.position = 'bottom')

See the R script “sample polydiv gen.R” to see a longer walkthrough of the workflow that generates this.

Using Census Roads

I switched from using NHPN data to census roads for polygons. The data’s neater and more complete, and the measures it generates will be a little more accurate. Here’s a sample call, loading and ascribing dividedness based on census urban arterials.

smplc
## Simple feature collection with 1 feature and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -70.3474 ymin: 43.54435 xmax: -69.97585 ymax: 43.72769
## Geodetic CRS:  WGS 84
##   statefp placefp  placens   geoid     name      namelsad lsad classfp pcicbsa
## 1      23   60545 00582683 2360545 Portland Portland city   25      C5       Y
##   pcinecta mtfcc funcstat    aland    awater    intptlat     intptlon
## 1        Y G4110        A 55863362 123988333 +43.6331574 -070.1853051
##                         geometry
## 1 MULTIPOLYGON (((-70.13284 4...
cosf <- geox::county.subset(st_bbox(smplc))

rds <- tigris::roads( state = cosf$statefp
                     ,county = cosf$countyfp
                     ,year = 2019) %>% 
  rename_with(tolower) %>% 
  st_transform(4326)

# just interstates
ints <- rds %>% filter(rttyp == 'I') 

# arterials
arts <- rds %>% filter(mtfcc %in%
           Hmisc::Cs(S1100, S1200, S1630))
# -> Census codes for road features are referenced here:
# https://www2.census.gov/geo/pdfs/reference/mtfccs2019.pdf

# allocate block groups to sub-area divisions, using arterials
nb.plc.div <- divM::gen.cross.tract.dividedness(
  region = smplc
  ,divs = arts
  ,nbd.query.fcn = tigris::block_groups
  ,region.id.colm = 'placefp'
  ,erase.water = F
  ,year = 2019
  )

nb.plc.div
## # A tibble: 54 x 3
##    geoid        poly.id perc.in.div
##    <chr>        <fct>         <dbl>
##  1 230050024003 1             0.893
##  2 230050020012 23            0.910
##  3 230050020022 23            0.415
##  4 230050017002 23            0.591
##  5 230050017003 23            0.925
##  6 230050021022 26            0.367
##  7 230050019003 26            0.996
##  8 230050020011 26            0.808
##  9 230050019004 26            0.996
## 10 230050021023 26            0.992
## # ... with 44 more rows
# get background tiles for vis this time
sttm <- visaux::get.stamen.bkg(smplc, zoom = 12)

nb.plc.divsf <- nb.plc.div %>% geox::attach.geos(query.fcn = tigris::block_groups)

ggmap(sttm) +
  geom_sf(
    data = nb.plc.divsf
    ,aes(fill = poly.id)
    ,color = 'white'
    ,alpha = .6
    ,inherit.aes = F
  ) +
  geom_sf(data = arts
          ,aes(color = factor(rttyp))
          ,size =  1.3
          ,inherit.aes = F
  ) +
  scale_fill_discrete(guide = 'none') +
  visaux::bbox2ggcrop(nb.plc.divsf) +
  theme_void() +
  theme(legend.position = 'bottom')