create_simplified_geometry.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating simplified geometries for administrative boundaries\n",
    "\n",
    "This Jupyter notebook can be used to generate simplified geometries of administrative areas.\n",
    "\n",
    "**Click on the cell below and click on \"Run\" to import the necessary libraries.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "import math\n",
    "from descartes.patch import PolygonPatch\n",
    "import pyproj\n",
    "from pyproj import CRS\n",
    "from shapely.geometry.geo import box, mapping\n",
    "from shapely.ops import cascaded_union, transform\n",
    "import json\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import osmnx as ox\n",
    "from ipywidgets import interact, interactive, FloatSlider\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "BLUE = '#6699cc'\n",
    "GRAY = '#999999'\n",
    "\n",
    "proj_epsg4326 = CRS(\"EPSG:4326\")\n",
    "proj_epsg3857 = CRS(\"EPSG:3857\")\n",
    "\n",
    "project_fwd = partial(\n",
    "    pyproj.transform,\n",
    "    proj_epsg4326,\n",
    "    proj_epsg3857)\n",
    "\n",
    "project_bwd = partial(\n",
    "    pyproj.transform,\n",
    "    proj_epsg3857,\n",
    "    proj_epsg4326)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 1: Selection of administrative areas\n",
    "\n",
    "The input field below the next cell expects Nominatim queries as input. The query can contain multiple places, separated by a semicolon. If no input field is shown, click on the next cell and click on \"Run\".\n",
    "You can test your query here: [https://nominatim.openstreetmap.org](https://nominatim.openstreetmap.org).\n",
    "\n",
    "If the results match your expectations, you can scroll down to Part 2:\n",
    "\n",
    "Example queries:\n",
    "\n",
    "* `Solothurn, Switzerland`\n",
    "* `Kanton Bern, Switzerland`\n",
    "* `Kanton Bern, Switzerland; Solothurn, Switzerland`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "923246cc8b93494f8942f2334bc9e10e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Text(value='Solothurn, Switzerland', description='query'), Button(description='Run Inter…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def get_boundaries(query):\n",
    "    print(\"Retrieved query: {}\".format(query))\n",
    "\n",
    "    queries = query.split(\";\")\n",
    "    results = ox.gdf_from_places(queries)\n",
    "\n",
    "    geoms = []\n",
    "    for index, row in results.iterrows():    \n",
    "        geoms.append(row['geometry'])\n",
    "    geom = cascaded_union(geoms)\n",
    "\n",
    "    # Plot query result\n",
    "    fig, ax = plt.subplots()\n",
    "    patch_original = PolygonPatch(geom, fc=GRAY, ec=GRAY, alpha=0.5, zorder=1)\n",
    "    ax.add_patch(patch_original)\n",
    "    for g in geoms:\n",
    "        patch_geom = PolygonPatch(g, fc=\"None\", ec='k', alpha=1.0, zorder=2)\n",
    "        ax.add_patch(patch_geom)\n",
    "        \n",
    "    minx, miny, maxx, maxy = geom.bounds\n",
    "    plt.xlim(minx, maxx)\n",
    "    plt.ylim(miny, maxy)\n",
    "    plt.show()\n",
    "    return geom\n",
    " \n",
    "w = interactive(get_boundaries, {'manual': True},\n",
    "                query='Solothurn, Switzerland')\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 2: Selecting parameters for simplification\n",
    "\n",
    "Using the sliders below the next cell, you can adjust the parameters of `box_size` and `buffer_distance`. If no sliders are shown, or the wrong area is displayed, select the next cell and click on \"Run\".  A simplified geometry is generated by creating rectangular boxes. Only the boxes that touch the original geometry with an extra buffer distance are kept. This method is not particularly efficent. Especially for larger geometries, a slow processing time should be expected.\n",
    "\n",
    "* `buffer_distance` sets the size of a buffer around an administrative boundary. [Units in  (WGS 84 / Pseudo-Mercator projection)]\n",
    "* `box_size` sets the size of the \"boxes\" that are used to generate the simplified geometries. [degree longitude/latitude]\n",
    "\n",
    "Following the plot of the original, respectively simplified geometry, the simplified geometry is printed in text form: If you are satisfied with results, you can copy this text. \n",
    "```\n",
    "{'coordinates': (((7.48, 47.1),\n",
    "                  (7.48, 47.08),\n",
    "                  ...\n",
    "                  (7.48, 47.1)),),\n",
    " 'type': 'Polygon'}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "422acaffa2374dcc892ea39e6ab44760",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.03, continuous_update=False, description='box_size', max=0.15, min=0…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "geom = w.result\n",
    "\n",
    "bounds = geom.bounds\n",
    "max_diff = round(max(bounds[2] - bounds[0], bounds[3] - bounds[1]), 2)\n",
    "default_box = max(round(max_diff * 0.05, 2), 0.01)\n",
    "max_box = max(0.1, default_box * 5)\n",
    "\n",
    "def create_simple_geometry(box_size, buffer_distance):\n",
    "\n",
    "    print(\"Processing\")\n",
    "    # Project to EPSG:3857 to apply buffer\n",
    "    if buffer_distance > 0:\n",
    "        geom_epsg3857 = transform(project_fwd, geom)\n",
    "        geom_epsg3857_buffered = geom_epsg3857.buffer(buffer_distance)\n",
    "        geom_buffered = transform(project_bwd, geom_epsg3857_buffered)\n",
    "    else:\n",
    "        geom_buffered = geom\n",
    "\n",
    "    # Create boxes\n",
    "    bounds = geom_buffered.bounds\n",
    "    # To save space, round bounding box to generate coordinates with less digits\n",
    "    minx = math.floor(bounds[0] / box_size) * box_size\n",
    "    miny = math.floor(bounds[1] / box_size) * box_size\n",
    "    maxx = math.ceil(bounds[2] / box_size) * box_size\n",
    "    maxy = math.ceil(bounds[3] / box_size) * box_size\n",
    "\n",
    "    boxes = []\n",
    "    for x in np.arange(minx, maxx, box_size):\n",
    "        for y in np.arange(miny, maxy, box_size):\n",
    "\n",
    "            pts = [x, y, x + box_size, y + box_size]\n",
    "            pts = map(lambda x: round(x, 6), pts)\n",
    "            b = box(*pts)\n",
    "            if b.intersects(geom_buffered):\n",
    "                boxes.append(b)\n",
    "\n",
    "    new_geom = cascaded_union(boxes)\n",
    "\n",
    "    # Remove unnecessary points\n",
    "    new_geom = new_geom.simplify(0.0001, preserve_topology=False)\n",
    "\n",
    "    geojson = mapping(new_geom)\n",
    "\n",
    "    # Plot result\n",
    "    patch_original = PolygonPatch(geom, fc=GRAY, ec=GRAY, alpha=0.3, zorder=1)\n",
    "    patch_buffer = PolygonPatch(geom_buffered, fc=GRAY, ec=GRAY, alpha=0.3, zorder=1)\n",
    "    patch_box = PolygonPatch(new_geom, fc=BLUE, ec=BLUE, alpha=0.5, zorder=2)\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.add_patch(patch_original)\n",
    "    ax.add_patch(patch_buffer)\n",
    "    ax.add_patch(patch_box)\n",
    "\n",
    "    plt.xlim(minx, maxx)\n",
    "    plt.ylim(miny, maxy)\n",
    "    plt.show()\n",
    "\n",
    "    print(json.dumps(geojson, indent=4))\n",
    "\n",
    "interactive(create_simple_geometry,\n",
    "         box_size=FloatSlider(value=default_box, min=0.01, max=max_box, step=0.01, continuous_update=False),\n",
    "         buffer_distance=FloatSlider(value=500.0, min=0.0, max=5000.0, step=100.0, continuous_update=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}