.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "traffic_assignment/_auto_examples/plot_forecasting.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_traffic_assignment__auto_examples_plot_forecasting.py: .. _example_usage_forecasting: Forecasting =========== In this example, we present a full forecasting workflow for the Sioux Falls example model. We start creating the skim matrices, running the assignment for the base-year, and then distributing these trips into the network. Later, we estimate a set of future demand vectors which are going to be the input of a future year assignnment with select link analysis. .. GENERATED FROM PYTHON SOURCE LINES 14-24 .. seealso:: Several functions, methods, classes and modules are used in this example: * :func:`aequilibrae.paths.Graph` * :func:`aequilibrae.paths.TrafficClass` * :func:`aequilibrae.paths.TrafficAssignment` * :func:`aequilibrae.distribution.Ipf` * :func:`aequilibrae.distribution.GravityCalibration` * :func:`aequilibrae.distribution.GravityApplication` * :func:`aequilibrae.distribution.SyntheticGravityModel` .. GENERATED FROM PYTHON SOURCE LINES 26-35 .. code-block:: Python # Imports from uuid import uuid4 from os.path import join from tempfile import gettempdir import pandas as pd from aequilibrae.utils.create_example import create_example .. GENERATED FROM PYTHON SOURCE LINES 37-44 .. code-block:: Python # We create the example project inside our temp folder fldr = join(gettempdir(), uuid4().hex) project = create_example(fldr) logger = project.logger .. GENERATED FROM PYTHON SOURCE LINES 45-49 Traffic assignment with skimming -------------------------------- In this step, we'll set the skims for the variable ``free_flow_time``, and execute the traffic assignment for the base-year. .. GENERATED FROM PYTHON SOURCE LINES 49-52 .. code-block:: Python from aequilibrae.paths import TrafficAssignment, TrafficClass .. GENERATED FROM PYTHON SOURCE LINES 53-72 .. code-block:: Python # We 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. # We grab the graph for cars graph = project.network.graphs["c"] # Let's say we want to minimize the free_flow_time graph.set_graph("free_flow_time") # And will skim time and distance while we are at it graph.set_skimming(["free_flow_time", "distance"]) # And we will allow paths to be computed going through other centroids/centroid connectors # required for the Sioux Falls network, as all nodes are centroids graph.set_blocked_centroid_flows(False) .. GENERATED FROM PYTHON SOURCE LINES 73-74 Let's get the demand matrix directly from the project record, and inspect what matrices we have in the project. .. GENERATED FROM PYTHON SOURCE LINES 74-77 .. code-block:: Python proj_matrices = project.matrices proj_matrices.list() .. raw:: html
name file_name cores procedure procedure_id timestamp description status
0 demand_omx demand.omx 1 None None 2020-11-24 08:47:18 Original data imported to OMX format
1 demand_mc demand_mc.omx 3 None None 2021-02-24 00:51:35 None
2 skims skims.omx 2 None None None Example skim
3 demand_aem demand.aem 1 None None 2020-11-24 08:46:42 Original data imported to AEM format


.. GENERATED FROM PYTHON SOURCE LINES 78-79 We get the demand matrix, and prepare it for computation .. GENERATED FROM PYTHON SOURCE LINES 79-82 .. code-block:: Python demand = proj_matrices.get_matrix("demand_omx") demand.computational_view(["matrix"]) .. GENERATED FROM PYTHON SOURCE LINES 83-84 Let's perform the traffic assignment .. GENERATED FROM PYTHON SOURCE LINES 84-113 .. code-block:: Python # Create the assignment class assigclass = TrafficClass(name="car", graph=graph, matrix=demand) assig = TrafficAssignment() # We start by adding the list of traffic classes to be assigned assig.add_class(assigclass) # Then we set these parameters, which an only be configured after adding one class to the assignment assig.set_vdf("BPR") # This is not case-sensitive # Then we set the volume delay function and its parameters assig.set_vdf_parameters({"alpha": "b", "beta": "power"}) # The capacity and free flow travel times as they exist in the graph assig.set_capacity_field("capacity") assig.set_time_field("free_flow_time") # And the algorithm we want to use to assign assig.set_algorithm("bfw") # Since we haven't checked the parameters file, let's make sure convergence criteria is good assig.max_iter = 1000 assig.rgap_target = 0.001 # we then execute the assignment assig.execute() .. rst-class:: sphx-glr-script-out .. code-block:: none car : 0%| | 0/24 [00:00
iteration rgap alpha warnings beta0 beta1 beta2
0 1 inf 1.000000 1.000000 0.000000 0.0
1 2 0.855075 0.328400 1.000000 0.000000 0.0
2 3 0.476346 0.186602 1.000000 0.000000 0.0
3 4 0.235513 0.241148 1.000000 0.000000 0.0
4 5 0.109241 0.818547 0.607382 0.392618 0.0


.. GENERATED FROM PYTHON SOURCE LINES 119-120 And we can also see the results of the assignment .. GENERATED FROM PYTHON SOURCE LINES 120-123 .. code-block:: Python results = assig.results() results.head() .. raw:: html
matrix_ab matrix_ba matrix_tot Preload_AB Preload_BA Preload_tot Congested_Time_AB Congested_Time_BA Congested_Time_Max Delay_factor_AB Delay_factor_BA Delay_factor_Max VOC_AB VOC_BA VOC_max PCE_AB PCE_BA PCE_tot
link_id
1 4502.545113 0.0 4502.545113 0.0 0.0 0.0 6.000822 0.0 6.000822 1.000137 0.0 1.000137 0.173842 0.0 0.173842 4502.545113 0.0 4502.545113
2 8222.240524 0.0 8222.240524 0.0 0.0 0.0 4.009141 0.0 4.009141 1.002285 0.0 1.002285 0.351326 0.0 0.351326 8222.240524 0.0 8222.240524
3 4622.925028 0.0 4622.925028 0.0 0.0 0.0 6.000913 0.0 6.000913 1.000152 0.0 1.000152 0.178490 0.0 0.178490 4622.925028 0.0 4622.925028
4 5897.692905 0.0 5897.692905 0.0 0.0 0.0 6.501414 0.0 6.501414 1.300283 0.0 1.300283 1.189487 0.0 1.189487 5897.692905 0.0 5897.692905
5 8101.860609 0.0 8101.860609 0.0 0.0 0.0 4.008617 0.0 4.008617 1.002154 0.0 1.002154 0.346182 0.0 0.346182 8101.860609 0.0 8101.860609


.. GENERATED FROM PYTHON SOURCE LINES 124-125 We can export our results to CSV or get a Pandas DataFrame, but let's put it directly into the results database .. GENERATED FROM PYTHON SOURCE LINES 125-127 .. code-block:: Python assig.save_results("base_year_assignment") .. GENERATED FROM PYTHON SOURCE LINES 128-129 And save the skims .. GENERATED FROM PYTHON SOURCE LINES 129-131 .. code-block:: Python assig.save_skims("base_year_assignment_skims", which_ones="all", format="omx") .. GENERATED FROM PYTHON SOURCE LINES 132-135 Trip distribution ----------------- First, let's have a function to plot the Trip Length Frequency Distribution. .. GENERATED FROM PYTHON SOURCE LINES 135-138 .. code-block:: Python from math import log10, floor import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 139-156 .. code-block:: Python def plot_tlfd(demand, skim, name): plt.clf() b = floor(log10(skim.shape[0]) * 10) n, bins, patches = plt.hist( np.nan_to_num(skim.flatten(), 0), bins=b, weights=np.nan_to_num(demand.flatten()), density=False, facecolor="g", alpha=0.75, ) plt.xlabel("Trip length") plt.ylabel("Probability") plt.title(f"Trip-length frequency distribution for {name}") return plt .. GENERATED FROM PYTHON SOURCE LINES 157-160 Calibration ~~~~~~~~~~~ We will calibrate synthetic gravity models using the skims for ``free_flow_time`` that we just generated .. GENERATED FROM PYTHON SOURCE LINES 160-164 .. code-block:: Python import numpy as np from aequilibrae.distribution import GravityCalibration .. GENERATED FROM PYTHON SOURCE LINES 165-166 We need the demand matrix and to prepare it for computation .. GENERATED FROM PYTHON SOURCE LINES 166-169 .. code-block:: Python demand = proj_matrices.get_matrix("demand_aem") demand.computational_view(["matrix"]) .. GENERATED FROM PYTHON SOURCE LINES 170-171 We also need the skims we just saved into our project .. GENERATED FROM PYTHON SOURCE LINES 171-176 .. code-block:: Python imped = proj_matrices.get_matrix("base_year_assignment_skims_car") # We can check which matrix cores were created for our skims to decide which one to use imped.names .. rst-class:: sphx-glr-script-out .. code-block:: none ['distance_blended', 'distance_final', 'free_flow_time_blended', 'free_flow_time_final'] .. GENERATED FROM PYTHON SOURCE LINES 177-181 Where ``free_flow_time_final`` is actually the congested time for the last iteration But before using the data, let's get some impedance for the intrazonals. Let's assume it is 75% of the closest zone. .. GENERATED FROM PYTHON SOURCE LINES 181-195 .. code-block:: Python imped_core = "free_flow_time_final" imped.computational_view([imped_core]) # If we run the code below more than once, we will be overwriting the diagonal values with non-sensical data # so let's zero it first np.fill_diagonal(imped.matrix_view, 0) # We compute it with a little bit of NumPy magic intrazonals = np.amin(imped.matrix_view, where=imped.matrix_view > 0, initial=imped.matrix_view.max(), axis=1) intrazonals *= 0.75 # Then we fill in the impedance matrix np.fill_diagonal(imped.matrix_view, intrazonals) .. GENERATED FROM PYTHON SOURCE LINES 196-198 Since we are working with an OMX file, we cannot overwrite a matrix on disk. So let's give it a new name to save. .. GENERATED FROM PYTHON SOURCE LINES 198-200 .. code-block:: Python imped.save(names=["final_time_with_intrazonals"]) .. GENERATED FROM PYTHON SOURCE LINES 201-202 This also updates these new matrices as those being used for computation .. GENERATED FROM PYTHON SOURCE LINES 202-204 .. code-block:: Python imped.view_names .. rst-class:: sphx-glr-script-out .. code-block:: none ['final_time_with_intrazonals'] .. GENERATED FROM PYTHON SOURCE LINES 205-206 Let's calibrate our Gravity Model .. GENERATED FROM PYTHON SOURCE LINES 206-221 .. code-block:: Python for function in ["power", "expo"]: gc = GravityCalibration(matrix=demand, impedance=imped, function=function, nan_as_zero=True) gc.calibrate() model = gc.model # We save the model model.save(join(fldr, f"{function}_model.mod")) _ = plot_tlfd(gc.result_matrix.matrix_view, imped.matrix_view, f"{function} model") # We can save the result of applying the model as well # We can also save the calibration report with open(join(fldr, f"{function}_convergence.log"), "w") as otp: for r in gc.report: otp.write(r + "\n") .. image-sg:: /traffic_assignment/_auto_examples/images/sphx_glr_plot_forecasting_001.png :alt: Trip-length frequency distribution for expo model :srcset: /traffic_assignment/_auto_examples/images/sphx_glr_plot_forecasting_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 222-223 And let's plot a trip length frequency distribution for the demand itself .. GENERATED FROM PYTHON SOURCE LINES 223-226 .. code-block:: Python plt = plot_tlfd(demand.matrix_view, imped.matrix_view, "demand") plt.show() .. image-sg:: /traffic_assignment/_auto_examples/images/sphx_glr_plot_forecasting_002.png :alt: Trip-length frequency distribution for demand :srcset: /traffic_assignment/_auto_examples/images/sphx_glr_plot_forecasting_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 227-232 Forecast -------- We create a set of 'future' vectors using some random growth factors. We apply the model for inverse power, as the trip frequency length distribution (TFLD) seems to be a better fit for the actual one. .. GENERATED FROM PYTHON SOURCE LINES 234-236 .. code-block:: Python from aequilibrae.distribution import Ipf, GravityApplication, SyntheticGravityModel .. GENERATED FROM PYTHON SOURCE LINES 237-241 Compute future vectors ~~~~~~~~~~~~~~~~~~~~~~ First thing to do is to compute the future vectors from our matrix. .. GENERATED FROM PYTHON SOURCE LINES 241-250 .. code-block:: Python origins = np.sum(demand.matrix_view, axis=1) destinations = np.sum(demand.matrix_view, axis=0) # Then grow them with some random growth between 0 and 10%, and balance them orig = origins * (1 + np.random.rand(origins.shape[0]) / 10) dest = destinations * (1 + np.random.rand(origins.shape[0]) / 10) dest *= orig.sum() / dest.sum() vectors = pd.DataFrame({"origins":orig, "destinations":dest}, index=demand.index[:]) .. GENERATED FROM PYTHON SOURCE LINES 251-255 IPF for the future vectors ~~~~~~~~~~~~~~~~~~~~~~~~~~ Let's balance the future vectors. The output of this step is going to be used later in the traffic assignment for future year. .. GENERATED FROM PYTHON SOURCE LINES 255-267 .. code-block:: Python args = { "matrix": demand, "vectors": vectors, "column_field": "destinations", "row_field": "origins", "nan_as_zero": True, } ipf = Ipf(**args) ipf.fit() .. GENERATED FROM PYTHON SOURCE LINES 268-269 When saving our vector into the project, we'll get an output that it was recored .. GENERATED FROM PYTHON SOURCE LINES 269-272 .. code-block:: Python ipf.save_to_project(name="demand_ipfd", file_name="demand_ipfd.aem") ipf.save_to_project(name="demand_ipfd_omx", file_name="demand_ipfd.omx") .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 273-277 Impedance ~~~~~~~~~ Let's get the base-year assignment skim for car we created before and prepare it for computation .. GENERATED FROM PYTHON SOURCE LINES 277-280 .. code-block:: Python imped = proj_matrices.get_matrix("base_year_assignment_skims_car") imped.computational_view(["final_time_with_intrazonals"]) .. GENERATED FROM PYTHON SOURCE LINES 281-282 If we wanted the main diagonal to not be considered... .. GENERATED FROM PYTHON SOURCE LINES 282-285 .. code-block:: Python # np.fill_diagonal(imped.matrix_view, np.nan) .. GENERATED FROM PYTHON SOURCE LINES 286-287 Now we apply the Synthetic Gravity model .. GENERATED FROM PYTHON SOURCE LINES 287-307 .. code-block:: Python for function in ["power", "expo"]: model = SyntheticGravityModel() model.load(join(fldr, f"{function}_model.mod")) outmatrix = join(proj_matrices.fldr, f"demand_{function}_model.aem") args = { "impedance": imped, "vectors": vectors, "row_field": "origins", "model": model, "column_field": "destinations", "nan_as_zero": True, } gravity = GravityApplication(**args) gravity.apply() # We get the output matrix and save it to OMX too, gravity.save_to_project(name=f"demand_{function}_modeled", file_name=f"demand_{function}_modeled.omx") .. GENERATED FROM PYTHON SOURCE LINES 308-309 We update the matrices table/records and verify that the new matrices are indeed there .. GENERATED FROM PYTHON SOURCE LINES 309-312 .. code-block:: Python proj_matrices.update_database() proj_matrices.list() .. raw:: html
name file_name cores procedure procedure_id timestamp description status
0 demand_omx demand.omx 1 None None 2020-11-24 08:47:18 Original data imported to OMX format
1 demand_mc demand_mc.omx 3 None None 2021-02-24 00:51:35 None
2 skims skims.omx 2 None None None Example skim
3 demand_aem demand.aem 1 None None 2020-11-24 08:46:42 Original data imported to AEM format
4 base_year_assignment_skims_car base_year_assignment_skims_car.omx 4 Traffic Assignment 440b5c8e4eb84aa48885321663ea23ae 2025-01-15 19:57:13.631339 Skimming for assignment procedure. Class car
5 demand_ipfd demand_ipfd.aem 1 Iterative Proportional fitting ad37ca792d7444448da0729e6c356925 2025-01-15 19:57:14.924919 None
6 demand_ipfd_omx demand_ipfd.omx 1 Iterative Proportional fitting ad37ca792d7444448da0729e6c356925 2025-01-15 19:57:14.924919 None
7 demand_power_modeled demand_power_modeled.omx 1 Synthetic gravity trip distribution 8f5238716b744dc18029c5fb841d6c01 2025-01-15 19:57:15.005479 Synthetic gravity trip distribution. POWER
8 demand_expo_modeled demand_expo_modeled.omx 1 Synthetic gravity trip distribution 129cfb7cae014a2a888c9125b4f35636 2025-01-15 19:57:15.098963 Synthetic gravity trip distribution. EXPO


.. GENERATED FROM PYTHON SOURCE LINES 313-316 Traffic assignment with Select Link Analysis -------------------------------------------- We'll perform traffic assignment for the future year. .. GENERATED FROM PYTHON SOURCE LINES 316-318 .. code-block:: Python logger.info("\n\n\n TRAFFIC ASSIGNMENT FOR FUTURE YEAR WITH SELECT LINK ANALYSIS") .. GENERATED FROM PYTHON SOURCE LINES 319-321 Let's get our future demand matrix, which corresponds to the IPF result we just saved, and see what is the core we ended up getting. It should be ``matrix``. .. GENERATED FROM PYTHON SOURCE LINES 321-324 .. code-block:: Python demand = proj_matrices.get_matrix("demand_ipfd") demand.names .. rst-class:: sphx-glr-script-out .. code-block:: none ['matrix'] .. GENERATED FROM PYTHON SOURCE LINES 325-326 Let's prepare our data for computation .. GENERATED FROM PYTHON SOURCE LINES 326-328 .. code-block:: Python demand.computational_view("matrix") .. GENERATED FROM PYTHON SOURCE LINES 329-330 The future year assignment is quite similar to the one we did for the base-year. .. GENERATED FROM PYTHON SOURCE LINES 330-355 .. code-block:: Python # So, let's create the assignment class assigclass = TrafficClass(name="car", graph=graph, matrix=demand) assig = TrafficAssignment() # Add at a list of traffic classes to be assigned assig.add_class(assigclass) assig.set_vdf("BPR") # Set the volume delay function and its parameters assig.set_vdf_parameters({"alpha": "b", "beta": "power"}) # Set the capacity and free flow travel times as they exist in the graph assig.set_capacity_field("capacity") assig.set_time_field("free_flow_time") # And the algorithm we want to use to assign assig.set_algorithm("bfw") # Once again we haven't checked the parameters file, so let's make sure convergence criteria is good assig.max_iter = 500 assig.rgap_target = 0.00001 .. GENERATED FROM PYTHON SOURCE LINES 356-357 Now we select two sets of links to execute select link analysis. .. GENERATED FROM PYTHON SOURCE LINES 357-362 .. code-block:: Python select_links = { "Leaving node 1": [(1, 1), (2, 1)], "Random nodes": [(3, 1), (5, 1)], } .. GENERATED FROM PYTHON SOURCE LINES 363-375 .. note:: As we are executing the select link analysis on a particular ``TrafficClass``, we should set the links we want to analyze. The input is a dictionary with string as keys and a list of tuples as values, so that each entry represents a separate set of selected links to compute. ``select_link_dict = {"set_name": [(link_id1, direction1), ..., (link_id, direction)]}`` The string name will name the set of links, and the list of tuples is the list of selected links in the form ``(link_id, direction)``, as it occurs in the :ref:`Graph `. Direction can be one of 0, 1, or -1, where 0 denotes bi-directionality. .. GENERATED FROM PYTHON SOURCE LINES 377-384 .. code-block:: Python # We call this command on the class we are analyzing with our dictionary of values assigclass.set_select_links(select_links) # we then execute the assignment assig.execute() .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/aequilibrae/paths/traffic_class.py:164: UserWarning: Input string name has a space in it. Replacing with _ warnings.warn("Input string name has a space in it. Replacing with _") car : 0%| | 0/24 [00:00` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_forecasting.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_forecasting.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_