.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "_auto_examples/creating_models/plot_create_zoning.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download__auto_examples_creating_models_plot_create_zoning.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr__auto_examples_creating_models_plot_create_zoning.py:


.. _create_zones:

Creating a zone system based on Hex Bins
========================================

In this example, we show how to create hex bin zones covering an arbitrary area.

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.

We also add centroid connectors to our network to make it a pretty complete example

.. GENERATED FROM PYTHON SOURCE LINES 21-22

Imports

.. GENERATED FROM PYTHON SOURCE LINES 22-30

.. code-block:: Python

    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








.. GENERATED FROM PYTHON SOURCE LINES 32-33

We create an empty project on an arbitrary folder

.. GENERATED FROM PYTHON SOURCE LINES 33-38

.. code-block:: Python

    fldr = join(gettempdir(), uuid4().hex)

    # Let's use the Nauru example project for display
    project = create_example(fldr, "nauru")








.. GENERATED FROM PYTHON SOURCE LINES 39-40

We said we wanted 100 zones

.. GENERATED FROM PYTHON SOURCE LINES 40-42

.. code-block:: Python

    zones = 100








.. GENERATED FROM PYTHON SOURCE LINES 43-45

Hex Bins using Spatialite
-------------------------

.. GENERATED FROM PYTHON SOURCE LINES 47-49

Spatialite requires a few things to compute hex bins.
One of them is the area you want to cover.

.. GENERATED FROM PYTHON SOURCE LINES 49-54

.. code-block:: Python

    network = project.network

    # So we use the convenient network method ``convex_hull()`` (it may take some time for very large networks)
    geo = network.convex_hull()








.. GENERATED FROM PYTHON SOURCE LINES 55-57

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 57-62

.. code-block:: Python

    zone_area = geo.area / zones
    # Since the area of the hexagon is **3 * sqrt(3) * side^2 / 2**
    # is side is equal to  **sqrt(2 * sqrt(3) * A/9)**
    zone_side = sqrt(2 * sqrt(3) * zone_area / 9)








.. GENERATED FROM PYTHON SOURCE LINES 63-67

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 67-78

.. code-block:: Python

    extent = network.extent()

    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 79-81

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 81-90

.. code-block:: Python

    grid = [p for p in grid.geoms if p.intersects(geo)]

    # 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
    nodes = network.nodes
    for i in range(1, 301):
        nd = nodes.get(i)
        nd.renumber(i + 1300)








.. GENERATED FROM PYTHON SOURCE LINES 91-93

Now we can add them to the model
And add centroids to them while we are at it

.. GENERATED FROM PYTHON SOURCE LINES 93-103

.. code-block:: Python

    zoning = project.zoning
    for i, zone_geo in enumerate(grid):
        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)









.. GENERATED FROM PYTHON SOURCE LINES 104-106

Centroid connectors
-------------------

.. GENERATED FROM PYTHON SOURCE LINES 108-120

.. code-block:: Python

    for zone_id, zone in zoning.all_zones().items():
        # We will connect for walk, with 1 connector per zone
        zone.connect_mode(mode_id="w", connectors=1)

        # And for cars, for cars with 2 connectors per zone
        # We also specify the link types we accept to connect to (can be used to avoid connection to ramps or freeways)
        zone.connect_mode(mode_id="c", link_types="ytrusP", connectors=2)

        # This takes a few minutes to compute, so we will break after processing the first 10 zones
        if zone_id >= 10:
            break





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/scipy/cluster/vq.py:602: UserWarning: One of the clusters is empty. Re-run kmeans with a different initialization.
      warnings.warn("One of the clusters is empty. "




.. GENERATED FROM PYTHON SOURCE LINES 121-123

Let's add special generator zones
We also add a centroid at the airport terminal

.. GENERATED FROM PYTHON SOURCE LINES 123-135

.. code-block:: Python

    nodes = project.network.nodes

    # Let's use some silly number for its ID, like 10,000, just so we can easily differentiate it
    airport = nodes.new_centroid(10000)
    airport.geometry = Point(166.91749582, -0.54472590)
    airport.save()

    # 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.
    # Distance here is in degrees, so 0.01 is equivalent to roughly 1.1km
    airport.connect_mode(airport.geometry.buffer(0.01), mode_id="c", link_types="ytrusP", connectors=1)








.. GENERATED FROM PYTHON SOURCE LINES 136-137

.. code-block:: Python

    project.close()








.. rst-class:: sphx-glr-timing

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


.. _sphx_glr_download__auto_examples_creating_models_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 <plot_create_zoning.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_create_zoning.py <plot_create_zoning.py>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_