.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "other_applications/_auto_examples/plot_create_zoning.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_other_applications__auto_examples_plot_create_zoning.py: .. _create_zones: Create a zone system based on Hex Bins ====================================== In this example, we show how to create hex bin zones covering an arbitrary area. We also add centroid connectors and a special generator zone to our network to make it a pretty complete example. We use the Nauru example to create roughly 100 zones covering the whole modeling area as delimited by the entire network. You are obviously welcome to create whatever zone system you would like, as long as you have the geometries for them. In that case, you can just skip the hex bin computation part of this notebook. .. GENERATED FROM PYTHON SOURCE LINES 37-40 .. admonition:: References * :ref:`Accessing project zones ` .. GENERATED FROM PYTHON SOURCE LINES 42-47 .. seealso:: Several functions, methods, classes and modules are used in this example: * :func:`aequilibrae.project.Zoning` * :func:`aequilibrae.project.network.Nodes` .. GENERATED FROM PYTHON SOURCE LINES 49-63 .. code-block:: Python # Imports from uuid import uuid4 from tempfile import gettempdir from os.path import join from math import sqrt from shapely.geometry import Point import shapely.wkb from aequilibrae.utils.create_example import create_example, list_examples from aequilibrae.utils.aeq_signal import simple_progress, SIGNAL s = SIGNAL(object) .. GENERATED FROM PYTHON SOURCE LINES 65-66 Let's print the list of examples that ship with AequilibraE .. GENERATED FROM PYTHON SOURCE LINES 66-68 .. code-block:: Python print(list_examples()) .. rst-class:: sphx-glr-script-out .. code-block:: none ['coquimbo', 'sioux_falls', 'nauru'] .. GENERATED FROM PYTHON SOURCE LINES 69-76 .. code-block:: Python # We create an empty project on an arbitrary folder fldr = join(gettempdir(), uuid4().hex) # Let's use the Nauru example project for display project = create_example(fldr, "nauru") .. GENERATED FROM PYTHON SOURCE LINES 77-78 We said we wanted 100 zones .. GENERATED FROM PYTHON SOURCE LINES 78-80 .. code-block:: Python zones = 100 .. GENERATED FROM PYTHON SOURCE LINES 81-83 Hex Bins using Spatialite ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 85-87 Spatialite requires a few things to compute hex bins. One of them is the area you want to cover. .. GENERATED FROM PYTHON SOURCE LINES 87-89 .. code-block:: Python network = project.network .. GENERATED FROM PYTHON SOURCE LINES 90-91 So we use the convenient network method ``convex_hull()`` (it may take some time for very large networks) .. GENERATED FROM PYTHON SOURCE LINES 91-93 .. code-block:: Python geo = network.convex_hull() .. GENERATED FROM PYTHON SOURCE LINES 94-96 The second thing is the side of the hex bin, which we can compute from its area. The approximate area of the desired hex bin is .. GENERATED FROM PYTHON SOURCE LINES 96-98 .. code-block:: Python zone_area = geo.area / zones .. GENERATED FROM PYTHON SOURCE LINES 99-101 Since the area of the hexagon is :math:`\frac{3\sqrt{3}}{2} * side^{2}` the side is equal to :math:`\sqrt{\frac{2\sqrt{3} * area}{9}}` .. GENERATED FROM PYTHON SOURCE LINES 101-103 .. code-block:: Python zone_side = sqrt(2 * sqrt(3) * zone_area / 9) .. GENERATED FROM PYTHON SOURCE LINES 104-108 Now we can run an SQL query to compute the hexagonal grid. There are many ways to create hex bins (including with a GUI on QGIS), but we find that using SpatiaLite is a pretty neat solution, for which we will use the entire network bounding box to make sure we cover everything. .. GENERATED FROM PYTHON SOURCE LINES 108-110 .. code-block:: Python extent = network.extent() .. GENERATED FROM PYTHON SOURCE LINES 111-120 .. code-block:: Python curr = project.conn.cursor() b = extent.bounds curr.execute( "select st_asbinary(HexagonalGrid(GeomFromWKB(?), ?, 0, GeomFromWKB(?)))", [extent.wkb, zone_side, Point(b[2], b[3]).wkb], ) grid = curr.fetchone()[0] grid = shapely.wkb.loads(grid) .. GENERATED FROM PYTHON SOURCE LINES 121-123 Since we used the bounding box, we have way more zones than we wanted, so we clean them by only keeping those that intersect the network convex hull. .. GENERATED FROM PYTHON SOURCE LINES 123-125 .. code-block:: Python grid = [p for p in grid.geoms if p.intersects(geo)] .. GENERATED FROM PYTHON SOURCE LINES 126-128 Let's re-number all nodes with IDs smaller than 300 to something bigger as to free space to our centroids to go from 1 to N. .. GENERATED FROM PYTHON SOURCE LINES 128-133 .. code-block:: Python nodes = network.nodes for i in range(1, 301): nd = nodes.get(i) nd.renumber(i + 1300) .. GENERATED FROM PYTHON SOURCE LINES 134-145 .. code-block:: Python # Now we can add them to the model and add centroids to them while we are at it. zoning = project.zoning for i, zone_geo in enumerate(simple_progress(grid, s, "Add zone centroids")): zone = zoning.new(i + 1) zone.geometry = zone_geo zone.save() # None means that the centroid will be added in the geometric point of the zone # But we could provide a Shapely point as an alternative zone.add_centroid(None) .. rst-class:: sphx-glr-script-out .. code-block:: none Add zone centroids : 0%| | 0/120 [00:00= 10: break .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( .. GENERATED FROM PYTHON SOURCE LINES 164-168 Special generator zones ----------------------- Let's add a special generator zone by adding a centroid at the airport terminal. .. GENERATED FROM PYTHON SOURCE LINES 170-171 Let's use some silly number for its ID, like 10,000, just so we can easily differentiate it .. GENERATED FROM PYTHON SOURCE LINES 171-175 .. code-block:: Python airport = nodes.new_centroid(10000) airport.geometry = Point(166.91749582, -0.54472590) airport.save() .. GENERATED FROM PYTHON SOURCE LINES 176-178 When connecting a centroid not associated with a zone, we need to tell AequilibraE what is the initial area around the centroid that needs to be considered when looking for candidate nodes. .. GENERATED FROM PYTHON SOURCE LINES 178-180 .. code-block:: Python airport.connect_mode(mode_id="c", link_types="ytrusP", connectors=1) .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/geopandas/array.py:403: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. warnings.warn( .. GENERATED FROM PYTHON SOURCE LINES 181-182 .. code-block:: Python project.close() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 4.938 seconds) .. _sphx_glr_download_other_applications__auto_examples_plot_create_zoning.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_create_zoning.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_create_zoning.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_create_zoning.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_