{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Create a zone system based on Hex Bins\n\nIn this example, we show how to create hex bin zones covering an arbitrary area.\n\nWe also add centroid connectors and a special generator zone to our network to make \nit a pretty complete example.\n\nWe use the Nauru example to create roughly 100 zones covering the whole modeling\narea as delimited by the entire network.\n\nYou are obviously welcome to create whatever zone system you would like, as long as\nyou have the geometries for them. In that case, you can just skip the hex bin computation\npart of this notebook.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. admonition:: References\n\n  * `Accessing project zones <project_zoning>`\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. seealso::\n    Several functions, methods, classes and modules are used in this example:\n\n    * :func:`aequilibrae.project.Zoning`\n    * :func:`aequilibrae.project.network.Nodes` \n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Imports\nfrom uuid import uuid4\nfrom tempfile import gettempdir\nfrom os.path import join\nfrom math import sqrt\nfrom shapely.geometry import Point\nimport shapely.wkb\n\nfrom aequilibrae.utils.create_example import create_example, list_examples\nfrom aequilibrae.utils.aeq_signal import simple_progress, SIGNAL\ns = SIGNAL(object)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's print the list of examples that ship with AequilibraE\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(list_examples())"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# We create an empty project on an arbitrary folder\nfldr = join(gettempdir(), uuid4().hex)\n\n# Let's use the Nauru example project for display\nproject = create_example(fldr, \"nauru\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We said we wanted 100 zones\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "zones = 100"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Hex Bins using Spatialite\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Spatialite requires a few things to compute hex bins.\nOne of them is the area you want to cover.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "network = project.network"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "So we use the convenient network method ``convex_hull()`` (it may take some time for very large networks)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "geo = network.convex_hull()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The second thing is the side of the hex bin, which we can compute from its area.\nThe approximate area of the desired hex bin is\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "zone_area = geo.area / zones"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Since the area of the hexagon is $\\frac{3\\sqrt{3}}{2} * side^{2}$\nthe side is equal to $\\sqrt{\\frac{2\\sqrt{3} * area}{9}}$\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "zone_side = sqrt(2 * sqrt(3) * zone_area / 9)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we can run an SQL query to compute the hexagonal grid.\nThere are many ways to create hex bins (including with a GUI on QGIS), but we find that\nusing SpatiaLite is a pretty neat solution, \nfor which we will use the entire network bounding box to make sure we cover everything.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "extent = network.extent()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "b = extent.bounds\nsql = \"select st_asbinary(HexagonalGrid(GeomFromWKB(?), ?, 0, GeomFromWKB(?)))\"\nwith project.db_connection as conn:\n    grid = conn.execute(sql, [extent.wkb, zone_side, Point(b[2], b[3]).wkb]).fetchone()[0]\n    grid = shapely.wkb.loads(grid)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Since we used the bounding box, we have way more zones than we wanted, so we clean them\nby only keeping those that intersect the network convex hull.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "grid = [p for p in grid.geoms if p.intersects(geo)]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's re-number all nodes with IDs smaller than 300 to something bigger as to free space to our\ncentroids to go from 1 to N.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nodes = network.nodes\nfor i in range(1, 301):\n    nd = nodes.get(i)\n    nd.renumber(i + 1300)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Now we can add them to the model and add centroids to them while we are at it.\nzoning = project.zoning\nfor i, zone_geo in enumerate(simple_progress(grid, s, \"Add zone centroids\")):\n    zone = zoning.new(i + 1)\n    zone.geometry = zone_geo\n    zone.save()\n    # None means that the centroid will be added in the geometric point of the zone\n    # But we could provide a Shapely point as an alternative\n    zone.add_centroid(None)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Centroid connectors\nLet's connect our zone centroids to the network.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for zone_id, zone in zoning.all_zones().items():\n    # We will connect for walk, with 1 connector per zone\n    zone.connect_mode(mode_id=\"w\", connectors=1)\n\n    # And for cars, for cars with 2 connectors per zone\n    # We also specify the link types we accept to connect to (can be used to avoid connection to ramps or freeways)\n    zone.connect_mode(mode_id=\"c\", link_types=\"ytrusP\", connectors=2)\n\n    # This takes a few minutes to compute, so we will break after processing the first 10 zones\n    if zone_id >= 10:\n        break"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Special generator zones\n\nLet's add a special generator zone by adding a centroid at the airport terminal.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's use some silly number for its ID, like 10,000, just so we can easily differentiate it\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "airport = nodes.new_centroid(10000)\nairport.geometry = Point(166.91749582, -0.54472590)\nairport.save()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "When connecting a centroid not associated with a zone, we need to tell AequilibraE what is the initial area around\nthe centroid that needs to be considered when looking for candidate nodes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "airport.connect_mode(mode_id=\"c\", link_types=\"ytrusP\", connectors=1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "project.close()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.10.18"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}