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

See also

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

# 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\).

theta = 0.011  # Distance factor

beta = 1.1  # PSL parameter

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.
2025-06-17 01:06:44,519;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/graph.py:248: UserWarning: Found centroids not present in the graph!
[  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108
 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
 127 128 129 130 131 132 133]
  warnings.warn("Found centroids not present in the graph!\n" + str(centroids[~present_centroids]))
2025-06-17 01:06:44,592;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
2025-06-17 01:06:44,679;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
2025-06-17 01:06:44,767;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/graph.py:248: UserWarning: Found centroids not present in the graph!
[  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108
 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
 127 128 129 130 131 132 133]
  warnings.warn("Found centroids not present in the graph!\n" + str(centroids[~present_centroids]))

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)
2025-06-17 01:06:44,833;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")

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. If we were interested in the original origin and destination IDs for each entry we could use subarea.post_process(keep_original_ods=True) instead. This will attach the true ODs from the select link OD matrix as part of the index. However, this will create a significantly larger, but more flexible matrix.

from aequilibrae.paths import SubAreaAnalysis

subarea = SubAreaAnalysis(graph, zones, mat)
subarea.rc.set_choice_set_generation("lp", max_routes=3, penalty=1.02, store_results=False)
subarea.rc.execute(perform_assignment=True)
demand = subarea.post_process()
demand
2025-06-17 01:06:45,323;INFO ; Created: 650 edge pairs from 26 edges
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/route_choice.py:581: 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(
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/route_choice.py:581: 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(
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/route_choice.py:581: 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(
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/route_choice.py:581: 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(
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/route_choice.py:309: UserWarning: found unreachable OD pairs, no choice sets generated for: [(64, 1), (64, 2), (64, 3), (64, 4), (64, 5), (64, 6), (64, 7), (64, 8), (64, 9), (64, 10), (64, 11), (64, 12), (64, 13), (64, 14), (64, 15), (64, 16), (64, 17), (64, 18), (64, 19), (64, 20), (64, 21), (64, 22), (64, 23), (64, 24), (64, 25), (64, 26), (64, 27), (64, 28), (64, 29), (64, 30), (64, 31), (64, 32), (64, 33), (64, 34), (64, 35), (64, 36), (64, 37), (64, 38), (64, 39), (64, 40), (64, 41), (64, 42), (64, 43), (64, 44), (64, 45), (64, 46), (64, 47), (64, 48), (64, 49), (64, 50), (64, 51), (64, 52), (64, 53), (64, 54), (64, 55), (64, 56), (64, 57), (64, 58), (64, 59), (64, 60), (64, 61), (64, 62), (64, 63), (64, 65), (64, 66), (64, 67), (64, 68), (64, 69), (64, 70), (64, 71), (64, 72), (64, 73), (64, 74), (64, 75), (64, 76), (64, 77), (64, 78), (64, 79), (64, 80), (64, 81), (64, 82), (64, 83), (64, 84), (64, 85), (64, 86), (64, 87), (64, 88), (64, 89), (64, 90), (64, 91), (64, 92), (64, 93), (64, 94), (64, 95), (64, 96), (64, 97), (64, 98), (64, 99), (64, 100), (64, 101), (64, 102), (64, 103), (64, 104), (64, 105), (64, 106), (64, 107), (64, 108), (64, 109), (64, 110), (64, 111), (64, 112), (64, 113), (64, 114), (64, 115), (64, 116), (64, 117), (64, 118), (64, 119), (64, 120), (64, 121), (64, 122), (64, 123), (64, 124), (64, 125), (64, 126), (64, 127), (64, 128), (64, 129), (64, 130), (64, 131), (64, 132), (64, 133)]
  self.__rc.batched(
demand
origin id destination id
29 29 10.000000
30 10.000000
31 10.000000
32 10.000000
33 10.000000
... ... ...
79892 67891 44.203742
72092 1780.272598
72161 472.216805
73381 44.203742
73541 44.203742

558 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
2025-06-17 01:06:57,660;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("lp", max_routes=3, penalty=1.02, store_results=False, seed=123)
rc.execute(perform_assignment=True)
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/route_choice.py:309: UserWarning: found unreachable OD pairs, no choice sets generated for: [(34, 67891), (52, 67891), (52, 73381), (52, 73432), (52, 77287), (52, 79297), (34, 73381), (34, 77287), (34, 79297), (49, 67891), (49, 72092), (49, 72161), (49, 73381), (49, 77287), (49, 79297), (29, 67891), (29, 72092), (29, 72161), (29, 73381), (29, 77287), (29, 79297), (57, 67891), (57, 72092), (57, 77287), (37, 67891), (37, 77287), (50, 67891), (50, 72092), (50, 72161), (50, 73381), (50, 77287), (50, 79297), (30, 67891), (30, 73381), (30, 77287), (30, 79297), (58, 67891), (58, 77287), (58, 79297), (38, 67891), (38, 77287), (51, 67891), (51, 73381), (31, 67891), (31, 72092), (31, 72161), (31, 73381), (51, 77287), (51, 79297), (31, 77287), (31, 79297), (59, 67891), (59, 77287), (59, 79297), (39, 67891), (39, 72092), (39, 77287), (67891, 29), (67891, 30), (67891, 31), (32, 67891), (67891, 32), (67891, 33), (32, 73381), (67891, 34), (67891, 37), (67891, 38), (32, 73432), (67891, 39), (67891, 40), (67891, 49), (67891, 50), (67891, 51), (67891, 52), (67891, 57), (67891, 58), (32, 77287), (67891, 59), (67891, 60), (32, 79297), (67891, 77287), (67891, 79297), (68671, 29), (68671, 30), (68671, 31), (68671, 32), (68671, 33), (68671, 34), (68671, 37), (68671, 38), (68671, 40), (68671, 49), (68671, 50), (68671, 51), (68671, 52), (68671, 57), (68671, 58), (68671, 59), (68671, 60), (68671, 67891), (68671, 72092), (68671, 72161), (68671, 73381), (68671, 73541), (72081, 77289), (72092, 29), (72092, 31), (72092, 33), (72092, 34), (72092, 37), (72092, 38), (72092, 39), (72092, 40), (72092, 49), (72092, 50), (72092, 57), (72092, 58), (72092, 59), (72092, 60), (72092, 77287), (72092, 79297), (40, 67891), (40, 72092), (72134, 77289), (72161, 29), (72161, 31), (72161, 33), (72161, 37), (72161, 38), (72161, 49), (72161, 50), (72161, 60), (40, 77287), (72161, 77287), (72161, 79297), (73381, 32), (73381, 52), (73394, 73381), (60, 67891), (60, 72092), (73394, 77289), (73432, 32), (73432, 34), (73432, 52), (60, 77287), (77285, 72081), (77285, 72134), (77285, 73394), (79892, 29), (79892, 30), (79892, 31), (79892, 32), (79892, 33), (79892, 34), (79892, 38), (79892, 40), (79892, 49), (79892, 50), (79892, 51), (79892, 52), (79892, 59), (79892, 60), (79892, 67891), (79892, 72092), (79892, 72161), (79892, 73381), (79892, 73541), (73432, 77289), (33, 67891), (33, 72092), (33, 77287), (33, 79297), (73565, 67891), (73565, 77289)]
  self.__rc.batched(

Let’s take the union of the zones GeoDataFrame as a polygon

poly = zones.union_all()
poly
<POLYGON ((-71.348 -29.993, -71.349 -29.993, -71.35 -29.992, -71.349 -29.991...>

And prepare the sub-area to plot.

subarea_zone = folium.Polygon(
    locations=[(x[1], x[0]) for x in poly.boundary.coords],
    fill_color="blue",
    fill_opacity=0.1,
    fill=True,
    weight=1,
)

We create a function to plot out link loads data more easily

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

    return loaded_links.explore(
        color="red",
        style_kwds={
            "style_function": lambda x: {
                "weight": x["properties"]["demand_tot"] * factor,
            }
        },
    )

And plot our data!

map = plot_results(rc.get_load_results())
subarea_zone.add_to(map)
map
<folium.folium.Map object at 0x7fb305585db0>

Sub-area further preparation#

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


Let’s 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


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

map = inner_links.explore(color="red", style_kwds={"weight": 4})
subarea_zone.add_to(map)
map
Make this Notebook Trusted to load map: File -> Trust Notebook


Manual sub-area analysis#

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

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'

Let’s prepare our graph once again

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)
2025-06-17 01:07:00,635;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/graph.py:248: UserWarning: Found centroids not present in the graph!
[  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108
 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
 127 128 129 130 131 132 133]
  warnings.warn("Found centroids not present in the graph!\n" + str(centroids[~present_centroids]))
2025-06-17 01:07:00,712;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
2025-06-17 01:07:00,802;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
2025-06-17 01:07:00,892;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/graph.py:248: UserWarning: Found centroids not present in the graph!
[  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108
 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
 127 128 129 130 131 132 133]
  warnings.warn("Found centroids not present in the graph!\n" + str(centroids[~present_centroids]))
2025-06-17 01:07:00,957;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)
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/route_choice.py:581: 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(
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/route_choice.py:581: 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(
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/route_choice.py:581: 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(
/home/runner/work/aequilibrae/aequilibrae/aequilibrae/paths/route_choice.py:581: 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("lp", max_routes=3, 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())
subarea_zone.add_to(map)
map
<folium.folium.Map object at 0x7fb2ed69e260>

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): 10.0})

and the load that exited the sub-area

exited
defaultdict(<class 'float'>, {(39, 20): 10.0, (37, 19): 10.0})

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): 3.5468435287475586, (35, 25): 10.0, (23, 36): 6.613396167755127, (23, 38): 3.386603832244873, (25, 37): 10.0, (39, 23): 6.453156471252441})

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 10.000000
21 23 3.546844
23 36 6.613396
38 3.386604
25 37 10.000000


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
2025-06-17 01:07:04,216;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations

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

Re-perform our assignment

rc = RouteChoice(graph)
rc.add_demand(demand)
rc.set_choice_set_generation("lp", max_routes=3, 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())
subarea_zone.add_to(map)
map
<folium.folium.Map object at 0x7fb2ed4895a0>
project.close()

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

Gallery generated by Sphinx-Gallery