{
 "metadata": {
  "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.8.6-final"
  },
  "orig_nbformat": 2,
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "'6.0alpha'"
      ]
     },
     "metadata": {},
     "execution_count": 1
    }
   ],
   "source": [
    "import os\n",
    "\n",
    "os.environ.get('GDS_ENV_VERSION')"
   ]
  },
  {
   "source": [
    "# Generate illustrations of tessellation\n",
    "\n",
    "This notebook contains one function `pipeline`, which for a given point (lat, lon) generates a sequence of seven images illustrating the process of creation of morphologicla tessellation within 250m buffer. The function is used to generate animations and figures in the blogpost."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "Matplotlib created a temporary config/cache directory at /tmp/matplotlib-dq2ou8_z because the default path (/home/jovyan/.config/matplotlib) is not a writable directory; it is highly recommended to set the MPLCONFIGDIR environment variable to a writable directory, in particular to speed up the import of Matplotlib and to better support multiprocessing.\n"
     ]
    }
   ],
   "source": [
    "import geopandas as gpd\n",
    "import momepy as mm\n",
    "import pygeos\n",
    "import numpy as np\n",
    "from scipy.spatial import Voronoi\n",
    "import pandas as pd\n",
    "from mapclassify import greedy\n",
    "import contextily as ctx\n",
    "import matplotlib.pyplot as plt\n",
    "from palettable.wesanderson import FantasticFox2_5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n  and should_run_async(code)\n"
     ]
    }
   ],
   "source": [
    "def pipeline(lat, lon, path, prefix, dist=250, figsize=(12, 12)):\n",
    "    point = (lat, lon)\n",
    "    gdf = ox.geometries.geometries_from_point(point, dist=dist, tags={'building':True})\n",
    "\n",
    "\n",
    "    gdf_projected = ox.projection.project_gdf(gdf)\n",
    "\n",
    "    bounds = gdf_projected.total_bounds\n",
    "    limit = Point(np.mean([bounds[0], bounds[2]]), np.mean([bounds[1], bounds[3]])).buffer(250)\n",
    "    blg = gpd.clip(gdf_projected, limit).explode()\n",
    "    bounds = limit.bounds\n",
    "\n",
    "    # figure 1 - aerial\n",
    "    fig, ax = plt.subplots(figsize=figsize)\n",
    "    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])\n",
    "    gpd.GeoSeries([limit.buffer(150).difference(limit)]).plot(ax=ax, color='white')\n",
    "    ctx.add_basemap(ax, crs=blg.crs, source=ctx.providers.Esri.WorldImagery)\n",
    "    ax.set_axis_off()\n",
    "    plt.savefig(path + prefix + \"01.png\", bbox_inches='tight')\n",
    "    plt.close()\n",
    "    print(\"Figure 1 saved to \" + path + prefix + \"01.png\")\n",
    "\n",
    "    # figure 2 - overlay\n",
    "    fig, ax = plt.subplots(figsize=figsize)\n",
    "    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])\n",
    "    gpd.GeoSeries([limit.buffer(150).difference(limit)]).plot(ax=ax, color='white')\n",
    "    ctx.add_basemap(ax, crs=blg.crs, source=ctx.providers.Esri.WorldImagery)\n",
    "    blg.plot(ax=ax, color='#0ea48f', edgecolor='k', alpha=.6)\n",
    "    ax.set_axis_off()\n",
    "    plt.savefig(path + prefix + \"02.png\", bbox_inches='tight')\n",
    "    plt.close()\n",
    "    print(\"Figure 2 saved to \" + path + prefix + \"02.png\")\n",
    "\n",
    "\n",
    "    # figure 3 - footprints\n",
    "    fig, ax = plt.subplots(figsize=figsize)\n",
    "    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])\n",
    "    blg.plot(ax=ax, color='#0ea48f', edgecolor='k').set_axis_off()\n",
    "    plt.savefig(path + prefix + \"03.png\", bbox_inches='tight')\n",
    "    plt.close()\n",
    "    print(\"Figure 3 saved to \" + path + prefix + \"03.png\")\n",
    "\n",
    "    shrinked = blg.buffer(-2)\n",
    "    shrinked = shrinked[~shrinked.is_empty]\n",
    "\n",
    "    # figure 4 - shrinked\n",
    "    fig, ax = plt.subplots(figsize=figsize)\n",
    "    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])\n",
    "    blg.plot(ax=ax, facecolor='none', linewidth=.5, edgecolor='k')\n",
    "    shrinked.plot(ax=ax, color='#0ea48f')\n",
    "    ax.set_axis_off()\n",
    "    plt.savefig(path + prefix + \"04.png\", bbox_inches='tight')\n",
    "    plt.close()\n",
    "    print(\"Figure 4 saved to \" + path + prefix + \"04.png\")\n",
    "\n",
    "    distance = 4\n",
    "    points = np.empty((0, 2))\n",
    "    ids = []\n",
    "    lines = shrinked.boundary.values.data\n",
    "    lengths = shrinked.length\n",
    "    for ix, line, length in zip(shrinked.index, lines, lengths):\n",
    "        if length > distance:\n",
    "            pts = pygeos.line_interpolate_point(\n",
    "                line,\n",
    "                np.linspace(0.1, length - 0.1, num=int((length - 0.1) // distance)),\n",
    "            )  # .1 offset to keep a gap between two segments\n",
    "            if len(pts) > 0:\n",
    "                points = np.append(points, pygeos.get_coordinates(pts), axis=0)\n",
    "                ids += [ix] * len(pts)\n",
    "\n",
    "    # figure 5 - points\n",
    "    fig, ax = plt.subplots(figsize=figsize)\n",
    "    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])\n",
    "    blg.plot(ax=ax, facecolor='none', linewidth=.5, edgecolor='k')\n",
    "    gpd.GeoSeries(pygeos.points(points)).plot(ax=ax, markersize=1, color='#0ea48f')\n",
    "    ax.set_axis_off()\n",
    "    plt.savefig(path + prefix + \"05.png\", bbox_inches='tight')\n",
    "    plt.close()\n",
    "    print(\"Figure 5 saved to \" + path + prefix + \"05.png\")\n",
    "\n",
    "    # add hull to resolve issues with infinity\n",
    "    # this is just a correction step ensuring the algorithm will work correctly\n",
    "    stop = points.shape[0]\n",
    "    series = gpd.GeoSeries(limit)\n",
    "    hull = series.geometry[[0]].buffer(500)\n",
    "    line = hull.boundary.values.data[0]\n",
    "    length = hull.length[0]\n",
    "    pts = pygeos.line_interpolate_point(\n",
    "        line,\n",
    "        np.linspace(0.1, length - 0.1, num=int((length - 0.1) // distance)),\n",
    "    )  # .1 offset to keep a gap between two segments\n",
    "    points = np.append(points, pygeos.get_coordinates(pts), axis=0)\n",
    "    ids += [-1] * len(pts)\n",
    "\n",
    "    voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "    vertices = pd.Series(voronoi_diagram.regions).take(voronoi_diagram.point_region)\n",
    "    polygons = []\n",
    "    for region in vertices:\n",
    "        if -1 not in region:\n",
    "            polygons.append(pygeos.polygons(voronoi_diagram.vertices[region]))\n",
    "        else:\n",
    "            polygons.append(None)\n",
    "\n",
    "    regions_gdf = gpd.GeoDataFrame(\n",
    "        {'unique_id': ids}, geometry=polygons\n",
    "    ).dropna()\n",
    "    regions_gdf = regions_gdf.loc[\n",
    "        regions_gdf['unique_id'] != -1\n",
    "    ]  # delete hull-based cells\n",
    "    voronoi_tessellation = gpd.clip(regions_gdf, limit)\n",
    "\n",
    "\n",
    "    # figure 6 - voronoi\n",
    "    fig, ax = plt.subplots(figsize=figsize)\n",
    "    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])\n",
    "    gpd.GeoSeries(pygeos.points(points[:stop])).plot(ax=ax, markersize=1, zorder=3, color='#0ea48f')\n",
    "    voronoi_tessellation.plot(ax=ax, facecolor='none', linewidth=.2, edgecolor='gray')\n",
    "    ax.set_axis_off()\n",
    "    plt.savefig(path + prefix + \"06.png\", bbox_inches='tight')\n",
    "    plt.close()\n",
    "    print(\"Figure 6 saved to \" + path + prefix + \"06.png\")\n",
    "\n",
    "    # figure 7 - tessellation\n",
    "    fig, ax = plt.subplots(figsize=figsize)\n",
    "    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])\n",
    "    blg = blg[blg.geom_type == 'Polygon']\n",
    "    blg = blg.reset_index(drop=True)\n",
    "    blg['uid'] = range(len(blg))\n",
    "    tessellation = mm.Tessellation(blg, 'uid', limit, verbose=False).tessellation\n",
    "    tessellation.plot(greedy(tessellation, strategy='smallest_last'), ax=ax, categorical=True, edgecolor='w', alpha=.6, cmap=FantasticFox2_5.mpl_colormap)\n",
    "    ax.set_axis_off()\n",
    "    plt.savefig(path + prefix + \"07.png\", bbox_inches='tight')\n",
    "    plt.close()\n",
    "    print(\"Figure 7 saved to \" + path + prefix + \"07.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "  and should_run_async(code)\n",
      "/opt/conda/lib/python3.8/site-packages/osmnx/footprints.py:44: UserWarning: The `footprints` module has been deprecated and will be removed in a future release. Instead, use the `geometries` module's `geometries_from_point` function, passing `tags={'building':True}`.\n",
      "  warnings.warn(msg)\n",
      "Figure 1 saved to ./la_01.png\n",
      "Figure 2 saved to ./la_02.png\n",
      "Figure 3 saved to ./la_03.png\n",
      "Figure 4 saved to ./la_04.png\n",
      "Figure 5 saved to ./la_05.png\n",
      "Figure 6 saved to ./la_06.png\n",
      "/opt/conda/lib/python3.8/site-packages/momepy/elements.py:474: UserWarning: Tessellation does not fully match buildings. 5 element(s) collapsed during generation - unique_id: {37, 133, 363, 375, 376}\n",
      "  warnings.warn(\n",
      "Figure 7 saved to ./la_07.png\n"
     ]
    }
   ],
   "source": [
    "pipeline(33.9488360, -118.2372975, path='./', prefix='la_', figsize=(15, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Figure 1 saved to ./bcn_01.png\n",
      "Figure 2 saved to ./bcn_02.png\n",
      "Figure 3 saved to ./bcn_03.png\n",
      "Figure 4 saved to ./bcn_04.png\n",
      "Figure 5 saved to ./bcn_05.png\n",
      "Figure 6 saved to ./bcn_06.png\n",
      "/opt/conda/lib/python3.8/site-packages/momepy/elements.py:474: UserWarning: Tessellation does not fully match buildings. 5 element(s) collapsed during generation - unique_id: {206, 277, 279, 345, 252}\n",
      "  warnings.warn(\n",
      "Figure 7 saved to ./bcn_07.png\n"
     ]
    }
   ],
   "source": [
    "pipeline(41.3907594, 2.1573404, path='./', prefix='bcn_', figsize=(15, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "  and should_run_async(code)\n",
      "Figure 1 saved to ./atl_01.png\n",
      "Figure 2 saved to ./atl_02.png\n",
      "Figure 3 saved to ./atl_03.png\n",
      "Figure 4 saved to ./atl_04.png\n",
      "Figure 5 saved to ./atl_05.png\n",
      "Figure 6 saved to ./atl_06.png\n",
      "/opt/conda/lib/python3.8/site-packages/momepy/elements.py:474: UserWarning: Tessellation does not fully match buildings. 1 element(s) collapsed during generation - unique_id: {21}\n",
      "  warnings.warn(\n",
      "Figure 7 saved to ./atl_07.png\n"
     ]
    }
   ],
   "source": [
    "pipeline(38.995888, -77.135073, path='./', prefix='atl_', figsize=(15, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "  and should_run_async(code)\n",
      "Figure 1 saved to ./bol_01.png\n",
      "Figure 2 saved to ./bol_02.png\n",
      "Figure 3 saved to ./bol_03.png\n",
      "Figure 4 saved to ./bol_04.png\n",
      "Figure 5 saved to ./bol_05.png\n",
      "Figure 6 saved to ./bol_06.png\n",
      "/opt/conda/lib/python3.8/site-packages/momepy/elements.py:474: UserWarning: Tessellation does not fully match buildings. 6 element(s) collapsed during generation - unique_id: {352, 808, 684, 685, 716, 175}\n",
      "  warnings.warn(\n",
      "Figure 7 saved to ./bol_07.png\n"
     ]
    }
   ],
   "source": [
    "pipeline(44.4942640, 11.3473233, path='./', prefix='bol_', figsize=(15, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "  and should_run_async(code)\n",
      "Figure 1 saved to ./bra_01.png\n",
      "Figure 2 saved to ./bra_02.png\n",
      "Figure 3 saved to ./bra_03.png\n",
      "Figure 4 saved to ./bra_04.png\n",
      "Figure 5 saved to ./bra_05.png\n",
      "Figure 6 saved to ./bra_06.png\n",
      "Figure 7 saved to ./bra_07.png\n"
     ]
    }
   ],
   "source": [
    "pipeline(-15.8038355, -47.8918796, path='./', prefix='bra_', figsize=(15, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ]
}