12  Multimodal and More

Code
import contextlib
import contextily as ctx
import pandarm as pdna
import matplotlib.pyplot as plt
import urbanaccess as ua
from access import Access
from geosnap import DataStore
from geosnap import io as gio
from libpysal.graph import Graph
from mapclassify import classify
from urbanaccess.gtfs.headways import headways
from urbanaccess import gtfsfeeds

%load_ext watermark
%watermark -a 'eli knaap' -iv
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
Author: eli knaap

matplotlib : 3.10.8
pandarm    : 0.0.0
libpysal   : 4.13.0
contextily : 1.6.2
access     : 1.1.9
geosnap    : 0.15.3
mapclassify: 2.10.0
urbanaccess: 0.2.2
Warning

This notebook currently requires my personal fork of the urbanaccess package

12.1 Multimodal Travel Networks

Public transit creates immense value for cities, and dramatically expands peoples’ ability to move through the urban landscape. Creating multimodal accessibility measures that visualize this phenomenon can be challenging because doing so requires creating a routable multi-modal network that integrates, e.g. both the pedestrian network and the bus/rail network. The urbanaccess package is designed for exactly that purpose (Blanchard & Waddell, 2017). It consumes OpenStreetMap data which encodes the pedestrian network and the General Transit Feed Specification (GTFS) data specification used to encode a public transit network, both of which are global data standards, meaning these analyses are relatively easy to replicate anywhere in the world.

Code
d = DataStore()
sd = gio.get_acs(d, county_fips="06073", years=[2021])
sd = sd[sd.n_total_pop > 0]
sdlodes = gio.get_lodes(d, county_fips="06073", years=[2021])
/Users/knaaptime/miniforge3/envs/urban_analysis/lib/python3.12/site-packages/geosnap/io/util.py:273: UserWarning: Unable to find local adjustment year for 2021. Attempting from online data
  warn(
/Users/knaaptime/miniforge3/envs/urban_analysis/lib/python3.12/site-packages/geosnap/io/constructors.py:218: UserWarning: Currency columns unavailable at this resolution; not adjusting for inflation
  warn(

First we will download a walking network from OpenStreetMap and lightly reformat it for the urbanaccess package by changing the impedance name from “length” to “distance” and dropping the geometry column that came along for the ride. Then, we import the walking network into the urbanaccess global network. Importantly, the urbanaccess global network always uses time as its network impedance factor, so we need to define a (constant/average) travel speed for the walk network that helps make the conversion.

Code
sdnet = gio.get_network_from_gdf(sd)
sdnet.edges_df = sdnet.edges_df.drop(columns=["geometry"])
sdnet.edges_df = sdnet.edges_df.rename(columns={"length": "distance"})
sdnet.impedance_names[0] = "distance"
ua_osm = ua.create_osm_net(
    osm_edges=sdnet.edges_df, osm_nodes=sdnet.nodes_df, travel_speed_mph=3
)
Generating contraction hierarchies with 16 threads.
Setting CH node vector of size 341930
Setting CH edge vector of size 911766
Range graph removed 937944 edges of 1823532
. 10% . 20% . 30% . 40% . 50% . 60% . 70% . 80% . 90% . 100%
Created OSM network with travel time impedance using a travel speed of 3 MPH. Took 0.00 seconds.

Now it is time to get a transit network stored in GTFS format. There are many places you can do this, but transitland is one good option where you can go to search feeds relevant to a study area. You can combine multiple GTFS feeds into a single network, but here we will just use a single provider. You can download the zip file yourself or allow urbanaccess to handle the input/output. Here we will do the latter, then read the downloaded GTFS data into a set of DataFrames (stored inside the loaded_feeds object). Then, we create a transit network object from the dataframes (which also adds it to the global network). Transit schedules are designed to help channel traffic in certain directions during different times of day to optimize the infrastructure (e.g. to help pipe people downtown during the morning peak period), so we need to give the transit network a time threshold defining when our hypothetical trip occurs (again, in this examine near the “AM Peak”).

The urbanaccess package prints a lot of useful diagnostic information to the console, but it can be overwhelming for the book, so we suppress it with a context manager below.

Code
with contextlib.redirect_stdout(None):

    gtfsfeeds.download(
        feed_name="mts",
        feed_url="http://www.sdmts.com/google_transit_files/google_transit.zip"
        
    )

    loaded_feeds = ua.gtfsfeed_to_df(
        "data/gtfsfeed_text/",
        validation=True,
        bbox=tuple(sd.total_bounds),
        remove_stops_outsidebbox=True,
        append_definitions=True,
    )

    tnet = ua.create_transit_net(
        gtfsfeeds_dfs=loaded_feeds,
        day="monday",
        timerange=["07:00:00", "10:00:00"],
        calendar_dates_lookup=None,
        time_aware=False,
    )
/Users/knaaptime/Dropbox/projects/urbanaccess/urbanaccess/gtfs/utils_format.py:1159: DtypeWarning: Columns (5) have mixed types. Specify dtype option on import or set low_memory=False.
  df = pd.read_csv(file)

Finally, we can integrate the walking and transit networks into a single routable system. For measuring transit accessibility, one important consideration is the level of service provided by each route, often measured by train or bus headways, which is the frequency that vehicles arrive at the stop for each route. In a multimodal network, estimating the anticipated headways is a way to help capture the amount of time passengers expect to wait at the station. Headways are often different depending on the direction of travel to optimize the fleet and help accommodate demand moving in different directions during morning and afternoon commuting periods. Thus, it is important to remember that a multimodal network is directed and often asymmetric so when creating a multimodal pandana network it is critical to set twoway=False1. This can seem counterintuitive, but ‘twoway’ in this context means each edge can be used in both directions (i.e. an undirected edge). Instead, we need a directed network, which means the parameter is False.

Code
# compute mean headways during the AM commute
heads = headways(loaded_feeds, ["07:00:00", "10:00:00"])

# generate intermediate connections between the pedestrian and transit networks
cnet = ua.integrate_network(
    urbanaccess_network=ua.ua_network, headways=True, urbanaccess_gtfsfeeds_df=heads
)

combined_net = pdna.Network(
    cnet.net_nodes["x"],
    cnet.net_nodes["y"],
    cnet.net_edges["from_int"],
    cnet.net_edges["to_int"],
    cnet.net_edges[["weight"]],
    twoway=False,
)

# project into a local coordinate system for better snapping
sd = sd.to_crs(sd.estimate_utm_crs())
sdlodes = sdlodes.to_crs(sd.crs)
combined_net = gio.project_network(combined_net, output_crs=sd.crs)
sdnet = gio.project_network(sdnet, output_crs=sd.crs)

In the combined network, the travel impedance is now measured in time, which means the network aggregation is performed in minutes. Here we will compute a ten-minute job-accessibility query, which is a relatively small timeframe. While the query itself is still relatively fast thanks to pandana, it requires a lot of memory to store a travel matrix this large which means I cannot compute an example larger than this one on my laptop without crashing the Python kernel. If you need to run a larger example, you can spin up a cloud machine with lots of memory. Also remember that access queries computed using pandana’s aggregate method are run against every node in the network, so if you need to compute a smaller subset of distances (i.e. using less memory), there is another example of generating cost matrices in Chapter 33.

Code
# network impedance is measured in minutes now
combined_net.precompute(10)

Calling precompute on a pandana.Network object computes the shortest path between all observations within the impedance threshold. This is the most computationally-intensive portion of the analysis, where we convert an extremely sparse travel network into a (much) denser point-to-point cost matrix, which can be expensive to store in memory. As an alternative, you could use something like r5py which provides an interface to the Conveyal routing product, and allows you to conduct your analysis in the cloud.

I prefer to avoid hosted solutions whenever possible because (1) I prefer control of the full analysis without sending data to a commercial server, (2) it is much easier to edit and manipulate the transport network when you only consume simple data formats like GTFS and OSM, (3) the modeling assumptions (like network impedance, etc) are explicit when you build your own network, and (4) while you can stand up an r5 instance on your local computer, doing so is a huge pain–and there is nothing worse than Java :P.

All the same, if you run into memory issues or you want a fully-hosted solution, r5py is another option. Since we’re done using the multimodal network for now, we will delete the object from the namespace below, which will release all that memory it currently occupies.

Code
# get the nearest network node for both geometries
nodes = combined_net.get_node_ids(
    sdlodes.centroid.geometry.x, sdlodes.centroid.geometry.y
)
tract_nodes = combined_net.get_node_ids(sd.centroid.geometry.x, sd.centroid.geometry.y)

sdlodes["node_ids"] = nodes
sd["node_ids"] = tract_nodes

# the original walk-only network
walk_nodes = sdnet.get_node_ids(sd.centroid.geometry.x, sd.centroid.geometry.y)
sd["walk_nodes"] = walk_nodes
block_walk_nodes = sdnet.get_node_ids(
    sdlodes.centroid.geometry.x, sdlodes.centroid.geometry.y
)

# set the job variables on walk and multimodal networks
sdnet.set(block_walk_nodes, sdlodes.total_employees.values, name="walk_jobs")
combined_net.set(nodes, sdlodes.total_employees.values, name="jobs")

# calculate 10 min job access on both networks
job_access = combined_net.aggregate(10, name="jobs")  # network impedance in time
walk_job_access = sdnet.aggregate(800, name="walk_jobs") # impedance in meters 

# make all measures available on the dataframe
sd = sd.merge(
    job_access.rename("access"), left_on="node_ids", right_index=True, how="left"
)
sd = sd.merge(
    walk_job_access.rename("walk_access"),
    left_on="walk_nodes",
    right_index=True,
    how="left",
)
sd["diff"] = (sd.access - sd.walk_access).values

# done using the multimodal network, so reclaim our memory
del combined_net
Removed 12657 rows because they contain missing values
Removed 12657 rows because they contain missing values

Using both the pedestrian network and the multimodal network we can now take a a look at the impact of transit. We will plot the accessibility measures for both networks and the difference between the two. For ease of comparison, we will use multimodal-based quantiles in the maps.

Code
transit_quantiles = classify(sd.access.values, scheme="quantiles", k=8)
sd.explore(
    "walk_access",
    scheme="user_defined",
    classification_kwds={"bins": transit_quantiles.bins},
    cmap="YlOrBr",
    tiles="CartoDB Positron",
    tooltip=["access", "walk_access"],
    style_kwds={"weight": 0.5},
)
Make this Notebook Trusted to load map: File -> Trust Notebook