Skip to main content
`Ctrl`+`K`

Note

Go to the end to download the full example code.

# Route Choice with sub-area analysis

In this example, we show how to perform sub-area analysis using route choice assignment, for a city in La Serena Metropolitan Area in Chile.

References

- Route Choice

See also

Several functions, methods, classes and modules are used in this example:

- `aequilibrae.paths.Graph()`

- `aequilibrae.paths.RouteChoice()`

- `aequilibrae.paths.SubAreaAnalysis()`

- `aequilibrae.matrix.AequilibraeMatrix()`

    # Imports
    from uuid import uuid4
    from tempfile import gettempdir
    from os.path import join
    import itertools

    import pandas as pd
    import numpy as np
    import folium

    from aequilibrae.utils.create_example import create_example

    # We create the example project inside our temp folder
    fldr = join(gettempdir(), uuid4().hex)

    project = create_example(fldr, "coquimbo")

    import logging
    import sys

    # We the project opens, we can tell the logger to direct all messages to the terminal as well
    logger = project.logger
    stdout_handler = logging.StreamHandler(sys.stdout)
    formatter = logging.Formatter("%(asctime)s;%(levelname)s ; %(message)s")
    stdout_handler.setFormatter(formatter)
    logger.addHandler(stdout_handler)

## Model parameters

We’ll set the parameters for our route choice model. These are the parameters that will be used to calculate the utility of each path. In our example, the utility is equal to \(distance * theta\), and the path overlap factor (PSL) is equal to \(beta\).

    # Distance factor
    theta = 0.011

    # PSL parameter
    beta = 1.1

Let’s build all graphs

    project.network.build_graphs()
    # We get warnings that several fields in the project are filled with NaNs.
    # This is true, but we won't use those fields.

    2024-11-04 00:39:42,148;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
    2024-11-04 00:39:42,228;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
    2024-11-04 00:39:42,320;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
    2024-11-04 00:39:42,410;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations

We grab the graph for cars

    graph = project.network.graphs["c"]

We also see what graphs are available

    project.network.graphs.keys()

    dict_keys(['b', 'c', 't', 'w'])

Let’s say that utility is just a function of distance. So we build our utility field as the \(distance * theta\).

    graph.network = graph.network.assign(utility=graph.network.distance * theta)

Prepare the graph with all nodes of interest as centroids

    graph.prepare_graph(graph.centroids)

    2024-11-04 00:39:42,476;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations

And set the cost of the graph the as the utility field just created

    graph.set_graph("utility")

We allow flows through “centroid connectors” because our centroids are not really centroids. If we have actual centroid connectors in the network (and more than one per centroid), then we should remove them from the graph.

    graph.set_blocked_centroid_flows(False)
    graph.graph.head()

| link_id | a_node | b_node | direction | id | distance | speed | travel_time | capacity | osm_id | lanes | modes | utility | __supernet_id__ | __compressed_id__
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
0 | 34806 | 0 | 8743 | 1 | 0 | 957.271672 | NaN | NaN | 99999.0 | NaN | NaN | ct | 10.529988 | 34273 | 0
1 | 34807 | 1 | 1180 | 1 | 1 | 1446.491380 | NaN | NaN | 99999.0 | NaN | NaN | ct | 15.911405 | 34275 | 1
2 | 34808 | 2 | 1222 | 1 | 2 | 1646.723896 | NaN | NaN | 99999.0 | NaN | NaN | ct | 18.113963 | 34277 | 2
3 | 34809 | 3 | 5674 | 1 | 3 | 1196.124943 | NaN | NaN | 99999.0 | NaN | NaN | ct | 13.157374 | 34279 | 3
4 | 34810 | 4 | 6339 | 1 | 4 | 536.920082 | NaN | NaN | 99999.0 | NaN | NaN | ct | 5.906121 | 34281 | 4

## Mock demand matrix

We’ll create a mock demand matrix with demand `10` for every zone and prepare it for computation.

    from aequilibrae.matrix import AequilibraeMatrix

    names_list = ["demand"]

    mat = AequilibraeMatrix()
    mat.create_empty(zones=graph.num_zones, matrix_names=names_list, memory_only=True)
    mat.index = graph.centroids[:]
    mat.matrices[:, :, 0] = np.full((graph.num_zones, graph.num_zones), 10.0)
    mat.computational_view()

## Sub-area preparation

We need to define some polygon for out sub-area analysis, here we’ll use a section of zones and create out polygon as the union of their geometry. It’s best to choose a polygon that avoids any unnecessary intersections with links as the resource requirements of this approach grow quadratically with the number of links cut.

    zones_of_interest = [29, 30, 31, 32, 33, 34, 37, 38, 39, 40, 49, 50, 51, 52, 57, 58, 59, 60]
    zones = project.zoning.data.set_index("zone_id")
    zones = zones.loc[zones_of_interest]
    zones.head()

| ogc_fid | area | name | population | employment | geometry
--- | --- | --- | --- | --- | --- | ---
zone_id |  |  |  |  |  |
29 | 29 | None | None | 4665.233898 | None | MULTIPOLYGON (((-71.32613 -29.97973, -71.32665...
30 | 30 | None | None | 4190.454424 | None | MULTIPOLYGON (((-71.33702 -29.98829, -71.33754...
31 | 31 | None | None | 4107.854460 | None | MULTIPOLYGON (((-71.3308 -29.98206, -71.33131 ...
32 | 32 | None | None | 4815.264246 | None | MULTIPOLYGON (((-71.34687 -29.99529, -71.34791...
33 | 33 | None | None | 4556.225468 | None | MULTIPOLYGON (((-71.33909 -29.97428, -71.33857...

## Sub-area analysis

From here there are two main paths to conduct a sub-area analysis, manual or automated. AequilibraE ships with a small class that handle most of the details regarding the implementation and extract of the relevant data. It also exposes all the tools necessary to conduct this analysis yourself if you need fine grained control.

### Automated sub-area analysis

We first construct out SubAreaAnalysis object from the graph, zones, and matrix we previously constructed, then configure the route choice assignment and execute it. From there the `post_process` method is able to use the route choice assignment results to construct the desired demand matrix as a DataFrame.

    from aequilibrae.paths import SubAreaAnalysis

    subarea = SubAreaAnalysis(graph, zones, mat)
    subarea.rc.set_choice_set_generation("lp", max_routes=5, penalty=1.02, store_results=False)
    subarea.rc.execute(perform_assignment=True)
    demand = subarea.post_process()
    demand

    2024-11-04 00:39:43,337;INFO ; Created: 650 edge pairs from 26 edges
    /opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/aequilibrae/paths/route_choice.py:464: UserWarning: Two input links map to the same compressed link in the network, removing superfluous link 31425 and direction -1 with compressed id 9483
      warnings.warn(
    /opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/aequilibrae/paths/route_choice.py:464: UserWarning: Two input links map to the same compressed link in the network, removing superfluous link 31425 and direction 1 with compressed id 14421
      warnings.warn(
    /opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/aequilibrae/paths/route_choice.py:464: UserWarning: Two input links map to the same compressed link in the network, removing superfluous link 21724 and direction 1 with compressed id 14421
      warnings.warn(
    /opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/aequilibrae/paths/route_choice.py:464: UserWarning: Two input links map to the same compressed link in the network, removing superfluous link 21724 and direction -1 with compressed id 9483
      warnings.warn(

|  | demand
--- | --- | ---
origin id | destination id |
29 | 29 | 10.000000
30 | 10.000000
31 | 10.000000
32 | 10.000000
33 | 10.000000
... | ... | ...
79892 | 67891 | 47.673554
72092 | 1604.171242
72161 | 461.310580
73381 | 33.512579
73541 | 33.512579

593 rows × 1 columns

We’ll re-prepare our graph but with our new “external” ODs.

    new_centroids = np.unique(demand.reset_index()[["origin id", "destination id"]].to_numpy().reshape(-1))
    graph.prepare_graph(new_centroids)
    graph.set_graph("utility")
    new_centroids

    2024-11-04 00:40:03,038;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations

    array([   29,    30,    31,    32,    33,    34,    37,    38,    39,
              40,    49,    50,    51,    52,    57,    58,    59,    60,
           61044, 67891, 68671, 72081, 72092, 72096, 72134, 72161, 73381,
           73394, 73432, 73506, 73541, 73565, 73589, 75548, 77285, 77287,
           77289, 79297, 79892])

We can then perform an assignment using our new demand matrix on the limited graph

    from aequilibrae.paths import RouteChoice

    rc = RouteChoice(graph)
    rc.add_demand(demand)
    rc.set_choice_set_generation("link-penalisation", max_routes=5, penalty=1.02, store_results=False, seed=123)
    rc.execute(perform_assignment=True)

And plot the link loads for easy viewing

    subarea_zone = folium.Polygon(
        locations=[(x[1], x[0]) for x in zones.unary_union.boundary.coords],
        fill_color="blue",
        fill_opacity=0.5,
        fill=True,
        stroke=False,
    )

    def plot_results(link_loads):
        link_loads = link_loads[link_loads.tot > 0]
        max_load = link_loads["tot"].max()
        links = project.network.links.data
        loaded_links = links.merge(link_loads, on="link_id", how="inner")

        loads_lyr = folium.FeatureGroup("link_loads")

        # Maximum thickness we would like is probably a 10, so let's make sure we don't go over that
        factor = 10 / max_load

        # Let's create the layers
        for _, rec in loaded_links.iterrows():
            points = rec.geometry.wkt.replace("LINESTRING ", "").replace("(", "").replace(")", "").split(", ")
            points = "[[" + "],[".join([p.replace(" ", ", ") for p in points]) + "]]"
            # we need to take from x/y to lat/long
            points = [[x[1], x[0]] for x in eval(points)]
            _ = folium.vector_layers.PolyLine(
                points,
                tooltip=f"link_id: {rec.link_id}, Flow: {rec.tot:.3f}",
                color="red",
                weight=factor * rec.tot,
            ).add_to(loads_lyr)
        long, lat = project.conn.execute("select avg(xmin), avg(ymin) from idx_links_geometry").fetchone()

        map_osm = folium.Map(location=[lat, long], tiles="Cartodb Positron", zoom_start=12)
        loads_lyr.add_to(map_osm)
        folium.LayerControl().add_to(map_osm)
        return map_osm

    map = plot_results(rc.get_load_results()["demand"])
    subarea_zone.add_to(map)
    map

    /home/runner/work/aequilibrae/aequilibrae/docs/source/examples/assignment_workflows/plot_subarea_analysis.py:169: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.
      locations=[(x[1], x[0]) for x in zones.unary_union.boundary.coords],

Make this Notebook Trusted to load map: File -> Trust Notebook

#### Sub-area further preparation

We take the union of this GeoDataFrame as our polygon.

    poly = zones.union_all()
    poly

    <POLYGON ((-71.348 -29.993, -71.349 -29.993, -71.35 -29.992, -71.349 -29.991...>

It’s useful later on to know which links from the network cross our polygon.

    links = project.network.links.data
    inner_links = links[links.crosses(poly.boundary)].sort_index()
    inner_links.head()

| ogc_fid | link_id | a_node | b_node | direction | distance | modes | link_type | name | speed_ab | speed_ba | travel_time_ab | travel_time_ba | capacity_ab | capacity_ba | osm_id | lanes_ab | lanes_ba | geometry
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
560 | 734 | 734 | 75548 | 75321 | 1 | 118.672175 | ct | primary | Avenida Costanera | 50.0 | NaN | None | None | NaN | NaN | 1.087246e+09 | 1.0 | NaN | LINESTRING (-71.31995 -29.96046, -71.32004 -29...
599 | 777 | 777 | 68671 | 77289 | 1 | 98.456777 | ct | trunk | None | NaN | NaN | None | None | NaN | NaN | 1.087084e+09 | 1.0 | NaN | LINESTRING (-71.31844 -29.96439, -71.31855 -29...
1821 | 3799 | 3799 | 73589 | 72898 | 1 | 38.022984 | ct | residential | Las Violetas | NaN | NaN | None | None | NaN | NaN | 7.925576e+08 | NaN | NaN | LINESTRING (-71.3339 -29.98877, -71.33423 -29....
9462 | 20823 | 20823 | 66687 | 73506 | 1 | 36.529160 | ct | residential | Las Orquídeas | NaN | NaN | None | None | NaN | NaN | 2.322910e+08 | NaN | NaN | LINESTRING (-71.33567 -29.98926, -71.33536 -29...
10112 | 21724 | 21724 | 77269 | 77266 | 0 | 105.722347 | ct | residential | Suiza | NaN | NaN | None | None | NaN | NaN | 2.253009e+08 | NaN | NaN | LINESTRING (-71.31896 -29.96747, -71.31838 -29...

As well as which nodes are interior.

    nodes = project.network.nodes.data.set_index("node_id")
    inside_nodes = nodes.sjoin(zones, how="inner").sort_index()
    inside_nodes.head()

| ogc_fid_left | is_centroid | modes | link_types | osm_id | geometry | zone_id | ogc_fid_right | area | name | population | employment
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
node_id |  |  |  |  |  |  |  |  |  |  |  |
29 | 70084 | 1 | ct | z | NaN | POINT (-71.3234 -29.97278) | 29 | 29 | None | None | 4665.233898 | None
30 | 70085 | 1 | ct | z | NaN | POINT (-71.33567 -29.98311) | 30 | 30 | None | None | 4190.454424 | None
31 | 70086 | 1 | ct | z | NaN | POINT (-71.33088 -29.97783) | 31 | 31 | None | None | 4107.854460 | None
32 | 70087 | 1 | ct | z | NaN | POINT (-71.34347 -29.98923) | 32 | 32 | None | None | 4815.264246 | None
33 | 70088 | 1 | ct | z | NaN | POINT (-71.33627 -29.96964) | 33 | 33 | None | None | 4556.225468 | None

Here we filter those network links to graph links, dropping any dead ends and creating a link_id, dir multi-index.

    g = (
        graph.graph.set_index("link_id")
        .loc[inner_links.link_id]
        .drop(graph.dead_end_links, errors="ignore")
        .reset_index()
        .set_index(["link_id", "direction"])
    )
    g.head()

|  | a_node | b_node | id | distance | speed | travel_time | capacity | osm_id | lanes | modes | utility | __supernet_id__ | __compressed_id__
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
link_id | direction |  |  |  |  |  |  |  |  |  |  |  |  |
734 | 1 | 33 | 12522 | 51 | 118.672175 | 50.0 | NaN | NaN | 1.087246e+09 | 1.0 | ct | 1.305394 | 703 | 49
777 | 1 | 20 | 36 | 25 | 98.456777 | NaN | NaN | NaN | 1.087084e+09 | 1.0 | ct | 1.083025 | 760 | 23
3799 | 1 | 32 | 11071 | 50 | 38.022984 | NaN | NaN | NaN | 7.925576e+08 | NaN | ct | 0.418253 | 2775 | 48
20823 | 1 | 8070 | 29 | 16220 | 36.529160 | NaN | NaN | NaN | 2.322910e+08 | NaN | ct | 0.401821 | 15745 | 10070
21724 | -1 | 13709 | 13710 | 29612 | 105.722347 | NaN | NaN | NaN | 2.253009e+08 | NaN | ct | 1.162946 | 16781 | 9234

#### Sub-area visualisation

Here we’ll quickly visualise what out sub-area is looking like. We’ll plot the polygon from our zoning system and the links that it cuts.

    points = [(link_id, list(x.coords)) for link_id, x in zip(inner_links.link_id, inner_links.geometry)]
    subarea_layer = folium.FeatureGroup("Cut links")

    for link_id, line in points:
        _ = folium.vector_layers.PolyLine(
            [(x[1], x[0]) for x in line],
            tooltip=f"link_id: {link_id}",
            color="red",
        ).add_to(subarea_layer)

    long, lat = project.conn.execute("select avg(xmin), avg(ymin) from idx_links_geometry").fetchone()

    map_osm = folium.Map(location=[lat, long], tiles="Cartodb Positron", zoom_start=12)

    subarea_zone.add_to(map_osm)

    subarea_layer.add_to(map_osm)
    _ = folium.LayerControl().add_to(map_osm)
    map_osm

Make this Notebook Trusted to load map: File -> Trust Notebook

### Manual sub-area analysis

In order to perform out analysis we need to know what OD pairs have flow that enters and/or exists our polygon. To do so we perform a select link analysis on all links and pairs of links that cross the boundary. We create them as tuples of tuples to make represent the select link AND sets.

    edge_pairs = {x: (x,) for x in itertools.permutations(g.index, r=2)}
    single_edges = {x: ((x,),) for x in g.index}
    f"Created: {len(edge_pairs)} edge pairs from {len(single_edges)} edges"

    'Created: 650 edge pairs from 26 edges'

Here we’ll construct and use the Route Choice class to generate our route sets

    from aequilibrae.paths import RouteChoice

We’ll re-prepare out graph quickly

    project.network.build_graphs()
    graph = project.network.graphs["c"]
    graph.network = graph.network.assign(utility=graph.network.distance * theta)
    graph.prepare_graph(graph.centroids)
    graph.set_graph("utility")
    graph.set_blocked_centroid_flows(False)

    2024-11-04 00:40:08,372;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
    2024-11-04 00:40:08,447;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
    2024-11-04 00:40:08,539;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
    2024-11-04 00:40:08,628;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
    2024-11-04 00:40:08,695;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations

This object construction might take a minute depending on the size of the graph due to the construction of the compressed link to network link mapping that’s required. This is a one time operation per graph and is cached. We need to supply a Graph and an AequilibraeMatrix or DataFrame via the `add_demand` method, if demand is not provided link loading cannot be preformed.

    rc = RouteChoice(graph)
    rc.add_demand(mat)

Here we add the union of edges as select link sets.

    rc.set_select_links(single_edges | edge_pairs)

    /opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/aequilibrae/paths/route_choice.py:464: UserWarning: Two input links map to the same compressed link in the network, removing superfluous link 31425 and direction -1 with compressed id 9483
      warnings.warn(
    /opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/aequilibrae/paths/route_choice.py:464: UserWarning: Two input links map to the same compressed link in the network, removing superfluous link 31425 and direction 1 with compressed id 14421
      warnings.warn(
    /opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/aequilibrae/paths/route_choice.py:464: UserWarning: Two input links map to the same compressed link in the network, removing superfluous link 21724 and direction 1 with compressed id 14421
      warnings.warn(
    /opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/aequilibrae/paths/route_choice.py:464: UserWarning: Two input links map to the same compressed link in the network, removing superfluous link 21724 and direction -1 with compressed id 9483
      warnings.warn(

For the sake of demonstration we limit out demand matrix to a few OD pairs. This filter is also possible with the automated approach, just edit the `subarea.rc.demand.df` DataFrame, however make sure the index remains intact.

    ods_pairs_of_interest = [
        (4, 39),
        (92, 37),
        (31, 58),
        (4, 19),
        (39, 34),
    ]
    ods_pairs_of_interest = ods_pairs_of_interest + [(x[1], x[0]) for x in ods_pairs_of_interest]
    rc.demand.df = rc.demand.df.loc[ods_pairs_of_interest].sort_index().astype(np.float32)
    rc.demand.df

|  | demand
--- | --- | ---
origin id | destination id |
4 | 19 | 10.0
39 | 10.0
19 | 4 | 10.0
31 | 58 | 10.0
34 | 39 | 10.0
37 | 92 | 10.0
39 | 4 | 10.0
34 | 10.0
58 | 31 | 10.0
92 | 37 | 10.0

Perform the assignment

    rc.set_choice_set_generation("link-penalisation", max_routes=5, penalty=1.02, store_results=False, seed=123)
    rc.execute(perform_assignment=True)

We can visualise the current links loads

    map = plot_results(rc.get_load_results()["demand"])
    subarea_zone.add_to(map)
    map

Make this Notebook Trusted to load map: File -> Trust Notebook

We’ll pull out just OD matrix results as well we need it for the post-processing, we’ll also convert the sparse matrices to SciPy COO matrices.

    sl_od = rc.get_select_link_od_matrix_results()
    edge_totals = {k: sl_od[k]["demand"].to_scipy() for k in single_edges}
    edge_pair_values = {k: sl_od[k]["demand"].to_scipy() for k in edge_pairs}

For the post processing, we are interested in the demand of OD pairs that enter or exit the sub-area, or do both. For the single enters and exists we can extract that information from the single link select link results. We also need to map the links that cross the boundary to the origin/destination node and the node that appears on the outside of the sub-area.

    from collections import defaultdict

    entered = defaultdict(float)
    exited = defaultdict(float)
    for (link_id, dir), v in edge_totals.items():
        link = g.loc[link_id, dir]
        for (o, d), load in v.todok().items():
            o = graph.all_nodes[o]
            d = graph.all_nodes[d]

            o_inside = o in inside_nodes.index
            d_inside = d in inside_nodes.index

            if o_inside and not d_inside:
                exited[o, graph.all_nodes[link.b_node]] += load
            elif not o_inside and d_inside:
                entered[graph.all_nodes[link.a_node], d] += load
            elif not o_inside and not d_inside:
                pass

Here he have the load that entered the sub-area

    entered

    defaultdict(<class 'float'>, {(34, 37): 10.0, (20, 39): 9.88913345336914, (23, 39): 0.1108664870262146})

and the load that exited the sub-area

    exited

    defaultdict(<class 'float'>, {(39, 20): 9.873265266418457, (37, 36): 0.07007550448179245, (39, 23): 0.12673552334308624, (37, 19): 9.929924011230469})

To find the load that both entered and exited we can look at the edge pair select link results.

    through = defaultdict(float)
    for (l1, l2), v in edge_pair_values.items():
        link1 = g.loc[l1]
        link2 = g.loc[l2]

        for (o, d), load in v.todok().items():
            o_inside = o in inside_nodes.index
            d_inside = d in inside_nodes.index

            if not o_inside and not d_inside:
                through[graph.all_nodes[link1.a_node], graph.all_nodes[link2.b_node]] += load

    through

    defaultdict(<class 'float'>, {(21, 23): 4.1274213790893555, (35, 22): 0.2188785821199417, (35, 25): 9.781121253967285, (22, 37): 0.435827374458313, (26, 36): 0.2188785821199417, (23, 36): 5.337530136108398, (23, 38): 4.443591117858887, (25, 37): 9.564172744750977, (39, 26): 0.435827374458313, (39, 23): 5.436751365661621})

With these results we can construct a new demand matrix. Usually this would be now transplanted onto another network, however for demonstration purposes we’ll reuse the same network.

    demand = pd.DataFrame(
        list(entered.values()) + list(exited.values()) + list(through.values()),
        index=pd.MultiIndex.from_tuples(
            list(entered.keys()) + list(exited.keys()) + list(through.keys()), names=["origin id", "destination id"]
        ),
        columns=["demand"],
    ).sort_index()
    demand.head()

|  | demand
--- | --- | ---
origin id | destination id |
20 | 39 | 9.889133
21 | 23 | 4.127421
22 | 37 | 0.435827
23 | 36 | 5.337530
38 | 4.443591

We’ll re-prepare our graph but with our new “external” ODs.

    new_centroids = np.unique(demand.reset_index()[["origin id", "destination id"]].to_numpy().reshape(-1))
    graph.prepare_graph(new_centroids)
    graph.set_graph("utility")
    new_centroids

    2024-11-04 00:40:13,047;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations

    array([19, 20, 21, 22, 23, 25, 26, 34, 35, 36, 37, 38, 39])

Re-perform our assignment

    rc = RouteChoice(graph)
    rc.add_demand(demand)
    rc.set_choice_set_generation("link-penalisation", max_routes=5, penalty=1.02, store_results=False, seed=123)
    rc.execute(perform_assignment=True)

And plot the link loads for easy viewing

    map = plot_results(rc.get_load_results()["demand"])
    subarea_zone.add_to(map)
    map

Make this Notebook Trusted to load map: File -> Trust Notebook

Total running time of the script: (0 minutes 34.603 seconds)

Gallery generated by Sphinx-Gallery