Subsetting by hwy types

There are many different ways to subset hwys within a region for generating the measures:

• By route type: “Interstate,” “State route,” “non-signed other arterial,” etc.. (SIGNT1 in data)
• By “Functional classification”: — can get to “limited-access” hwys this way. (FCLASS in data)
• By spatial relationship: Interstates and all hwys (of X types) that intersect interstates

# get sample area & trim hwys
(phl <- czs[grepl("Philadelphia",czs$cz_name),] %>% divM::region.reorg("cz")) 

## Simple feature collection with 1 feature and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 1177577 ymin: 4706815 xmax: 1342598 ymax: 4902574
## CRS:            +proj=lcc +lon_0=-90 +lat_1=33 +lat_2=45
## # A tibble: 1 x 4
##   region.id region.name  region.type                                    geometry
##   <chr>     <chr>        <chr>                                     <POLYGON [m]>
## 1 19700     Philadelphia cz          ((1326289 4770204, 1326012 4769725, 132542~

hwys <- hwys %>% st_crop(phl)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

# use FCLASS to create a hwy subset
lac <- hwys %>% filter(FCLASS %in% divM::lac_codes |
                          SIGNT1 == "I") # all limited-access

# spatial group hwy data
hwys <- hwys %>% denode.lines(group.cols = c("SIGNT1", "SIGN1"))
lac <- lac %>% denode.lines(group.cols = c("SIGNT1", "SIGN1"))

# Get hwy subsets for philadelphia: ----
# all hwys
all <- subset.polys.divs(region = phl,
                               div.sf = hwys)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

# all interstates
int <- subset.polys.divs(phl, hwys,
                               div.identifier.column = "SIGNT1",
                               always.include = "I", 
                               include.intersecting = F) 

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

 # all interstates and hwys that intersect interstates
tint <- subset.polys.divs(phl,
                                hwys,
                                div.identifier.column = "SIGNT1",
                                always.include = "I",
                                include.intersecting = T)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

lac <- subset.polys.divs(phl,
                                lac)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

# visualize differences
plot(all['SIGN1'], main = "all philly hwys", lwd=2)

plot(int['SIGN1'], main = "all philly interstates", lwd=2)

plot(tint['SIGN1'], main = "all philly interstates+touching",lwd=2)

plot(lac['SIGN1'], main = "all philly LAC",lwd=2)

“Filling gaps” / The linebreak problem

Sometimes a hwy has a minuscule break. Often the hwy classification changes for an interchange or a short stretch. For example, I676 in Camden.

I run a helper fcn to fill these before generating polygons. The logic of the function is: If a line segment has an endpoint within x threshold of another segment, and they have the same hwy identifiers, I connect the segments.

Default threshold for filling gap is 200 meters. The argument threshold can be passed to any wrapper function that fills gaps to change this.

Below are the 3 interstates in Philadelphia CZ that have these gaps.

# filling gaps in Philly interstates:
viz.fix <- Fix.all.hwys(int, return.gap.map = T, verbose = T)

viz.fix$I676

viz.fix$I76

viz.fix$I95

To continue creating the polygon measures, create versions of the above hwy subsets that have their gaps filled. Suffix names with ‘.c’

int.c <- suppressWarnings(Fix.all.hwys(int, verbose = F))
tint.c <- suppressWarnings(Fix.all.hwys(tint, verbose = F))
lac.c <- suppressWarnings(Fix.all.hwys(lac, verbose = F))

Creating polygon measures

From here, creating the polygon measure is pretty straightforward. Fcn polygonal.div does the work; it takes linear division data and a polygon representing study area.

Let’s see how it looks based on the different hwy subsets and whether or not the gaps are filled. I use CZ as boundary.

The other core parameter for generating these measures is the area floor— the minimum size of each polygon subdivisions. If a polygon is formed that is below this area, it is filtered out and not counted towards the total. By default this floor is 1/2 square km.

# Summary table:
Summary.table <- 
  tibble( gaps_filled = c(FALSE, TRUE)
           ,interstates_only = c(polygonal.div(phl, int)$n.polys,
                                 polygonal.div(phl, int.c)$n.polys)
           ,interstates_intersecting = c(polygonal.div(phl, tint)$n.polys,
                                         polygonal.div(phl, tint.c)$n.polys)
           ,limited_access = c(polygonal.div(phl, lac)$n.polys,
                                 polygonal.div(phl, lac.c)$n.polys)
           )

knitr::kable(Summary.table)

# to view the polygons as map, ask the same function to return.sf to map
polygonal.div(phl, lac.c, return.sf = T) %>%
  select("id") %>%
  plot(main = "LAC polygons in philly")

Wrapping Up

To end on a pun, a single function, Polys.wrapper, can replicate this whole workflow: of subsetting the hwys, filling the gaps, and generating the polygons. It will pass on any arguments to its wrapped fcns. See full options with ?Polys.wrapper.

suppressWarnings( Polys.wrapper(phl, lac,
              fill.gaps = T, verbose = F) )

## subsetting done
## gaps filled

## # A tibble: 1 x 4
##   region.id region.type region.name  n.polys
##   <chr>     <chr>       <chr>          <int>
## 1 19700     cz          Philadelphia      34