{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "For this project I create a voronoi diagram on the map based on data points (or point of interests), voronoi diagram have applications in almost all areas of science and engineering. For geospatial use case, it is useful to tell us the closest point of interest (POI) by representing each POI with a dot inside a polygon shape. So within a polygon, the closest POI is definitely the dot inside the polygon.\n", "\n", "Ok, let's start the project. As usual, for the first step let's import all necessary packages" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import folium\n", "import numpy as np\n", "import geopandas as gpd\n", "import osmnx as ox\n", "import matplotlib.pyplot as plt\n", "from shapely.ops import cascaded_union\n", "from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area\n", "from geovoronoi import voronoi_regions_from_coords, points_to_coords\n", "from shapely.geometry import Point, LineString, Polygon" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Insert data points/POI into geodaframe" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
geometry
0POINT (106.64409 -6.30529)
1POINT (106.65326 -6.30131)
2POINT (106.63775 -6.28477)
3POINT (106.66506 -6.28460)
4POINT (106.62758 -6.28352)
5POINT (106.64136 -6.27659)
6POINT (106.62597 -6.30364)
\n", "
" ], "text/plain": [ " geometry\n", "0 POINT (106.64409 -6.30529)\n", "1 POINT (106.65326 -6.30131)\n", "2 POINT (106.63775 -6.28477)\n", "3 POINT (106.66506 -6.28460)\n", "4 POINT (106.62758 -6.28352)\n", "5 POINT (106.64136 -6.27659)\n", "6 POINT (106.62597 -6.30364)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gdf = gpd.GeoDataFrame()\n", "gdf = gdf.append({'geometry': Point(106.644085,-6.305286)}, ignore_index=True)\n", "\n", "gdf = gdf.append({'geometry': Point(106.653261,-6.301309)}, ignore_index=True)\n", "gdf = gdf.append({'geometry': Point(106.637751,-6.284774)}, ignore_index=True)\n", "gdf = gdf.append({'geometry': Point(106.665062,-6.284598)}, ignore_index=True)\n", "gdf = gdf.append({'geometry': Point(106.627582,-6.283521)}, ignore_index=True)\n", "gdf = gdf.append({'geometry': Point(106.641365,-6.276593)}, ignore_index=True)\n", "gdf = gdf.append({'geometry': Point(106.625972,-6.303643)}, ignore_index=True)\n", "gdf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Determine the coverage area of the voronoi diagram & save it to geodaframe" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "area_max_lon = 106.670929\n", "area_min_lon = 106.619602\n", "area_max_lat = -6.275227\n", "area_min_lat = -6.309795\n", "\n", "lat_point_list = [area_min_lat, area_max_lat,area_max_lat,area_min_lat]\n", "lon_point_list = [area_min_lon, area_min_lon, area_max_lon, area_max_lon]\n", "\n", "polygon_geom = Polygon(zip(lon_point_list, lat_point_list))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
geometry
0POLYGON ((106.61960 -6.30980, 106.61960 -6.275...
\n", "
" ], "text/plain": [ " geometry\n", "0 POLYGON ((106.61960 -6.30980, 106.61960 -6.275..." ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "boundary = gpd.GeoDataFrame()\n", "boundary = boundary.append({'geometry': polygon_geom}, ignore_index=True)\n", "boundary" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Convert both dataframes to Web Mercator projection" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/anaconda3/lib/python3.7/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n", " return _prepare_from_string(\" \".join(pjargs))\n" ] } ], "source": [ "gdf.crs = {'init' :'epsg:3395'}\n", "boundary.crs = {'init' :'epsg:3395'}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Convert the boundary geometry into a union of the polygon and POI data to NumPy array of coordinates." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "boundary_shape = cascaded_union(boundary.geometry)\n", "coords = points_to_coords(gdf.geometry)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "boundary_shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate voronoi regions" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, boundary_shape)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEYCAYAAABV8iGRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhUUlEQVR4nO3daYxk63kX8P9zau+9p6e7Z+np6Zm5s/XM3Hvt62vPNTZxosgWYBRzpYAUB+PEhiAQIAgCJVgBBHxIPqAAIoIvIEwSIwQC4YU4Il7kLTg2iZc7+9zp7up9rb3q1FlePlR1T92Znpnq7qp6zznv/ye1pqeruurp6anzr/Oe931eUUqBiIjMY+kugIiI9GAAEBEZigFARGQoBgARkaEYAEREhmIAEBEZigFAxhGRj4vI70ftuYgOSrgOgDpFRH4PwHeVUr/2xNd/BsC/BzCllHK1FEdET+EZAHXSfwLw8yIiT3z9LwP4nYMc/EUk3tHKDikodRB1AwOAOul/AhgD8MHdL4jIKICPAvisiKRE5DdFZLn58Zsikmre70Misigi/1BEVgH8xzbv/8sisi4iKyLyCy3POywinxWRDRGZF5HPiIjVvO2TIvLN/X4AEZkRESUinxKRBQBfaX79F0XktojsiMiXReRsy/d8WETuikheRH5LRL4uIp/e77lE5P0i8kfN+/6RiLy/5bavicg/E5FviUhRRH5fRI43b0uLyG+LyJaI5JrfO3nE3xcZjgFAHaOUqgL4rwA+0fLlvwjgjlLqBwD+EYCbAF4F8AqA9wL4TMt9TwA4BuAsgL/W5v2HAZwG8CkA/7YZOADwb5q3nQfwE82afgHt+wkAVwF8pDmE9asA3gQwDuAbAD4HAM0D9H8D8CtohN9dAO/f7wFF5BiALwL41837/ksAXxSRsZa7/VyzzgkASQB/v/n1v9L8ec40v/evA6ge4OcheppSih/86NgHgA8AyAFIN//+LQB/t/n5QwB/tuW+HwEw1/z8QwDqu9/X5v2rAOItt6+jERix5mPNttz2SwC+1vz8kwC++Yz6ZwAoAOdbvva/AXyq5e8WgAoaQfUJAN9puU0AZAF8+snnQmMo7LtPPN93AHyy+fnXAHym5ba/AeD3mp//IoBvA3hZ9++YH9H54BkAdZRS6psANgF8TEQuoPGu/XebN58CMN9y9/nm13ZtKKVqLX9/0f231DuvK1QADAA4DiCxz/eePsCPkm35/CyAf9UceskB2EbjQH+6Wc/efZVSCsDiMx7zyZ9nv7pWWz7f/XkA4D8D+DKA/9IcDvsNEUkc4OchegoDgLrhs2i8M/55AF9WSq01v76MxsF013Tza7uenJL2ovs/yyYAZ5/vXWrje/erJQvgl5RSIy0fGaXUtwGsAJjavWPzAvgU9vfkz9N2XUopRyn1T5VSs2gMMX0U7xxqIzowBgB1w2cB/DSAv4rGzKBdnwPwGREZb46d/xqA337O4xz0/gAApZSHxrWIfyEig80Ltn+vne99hn8H4FdE5Bqwd4H5Z5u3fRHADRH5WHPG0N9E49rEfr4E4JKI/JyIxEXkLwGYBfCFFxUgIj8pIjdEJAaggEbA+Yf8eYgAMACoC5RSc2iMV/cD+F8tN/1zAN8D8EMAPwLw/5pfe5aD3r/V3wJQBvA2gG+iMQz1H9r9GVoppf4HgF9HY/ilAODHAP5M87ZNAD8L4DcAbKFxQP8eAHufx9lC4537Lzfv+w8AfLT5GC9yAo2LzQUAtwF8HY1hIaJD40Iwog5qTjVdBPBxpdRXdddD9Dw8AyA6IhH5iIiMNNco/CoaF4j/UHNZRC/EACA6ujfQmLK6CeDPA/iYaqyJIAo0DgERERmKZwBERIY6UKOr48ePq5mZmS6VQkREnfb9739/Uyk1vt9tBwqAmZkZfO973+tMVURE1HUi8uTq8z0cAiIiMhQDgIjIUAwAIiJDMQCIiAzFACAiMhQDgIjIUAwAIiJDMQCIiAzFACAiMtSBVgIfWXUHcNgk8Sjmq2kU3ZjuMihCpq11DKGiuwzaT2YUGJnu2sP3NAA+8Rd+Gucn+nr5lJHzo03gDwrd+w9BZhH4+LR8Hv0xV3cptI8frnn42togtre3u/L4PQ2AUrmEf/KPf72XTxk5PoDfWZ7CZo2jd3R0l1Jb+HPey7rLoGeZvAa59jNde3geRULGAvCnx0u6y6CIuBpb1l0CacQACKEziRyOp33dZVDICXxMuc9sFEkGYACEUOMsoKi7DAq5K6ltJH1bdxmkEQMgpM4k8hhPe7rLoBC7ElvSXQJpxgAIKQvABw24FlC2XXg+963utBh8nHazussgzXq7DoA6qnEWMICNWjTXBXzlzjreWs6jPxnHx29OIxWP5s+pw5X0FhIuh39MxzOAEIv6jKA7qwX4Cqi5HtYKPFh10mWLwz/EAAi9qUQek5loXguYPTkEESCTiGFyKKW7nMiIi49THP4hcAgo9CwAHzhexH/PjugupeM+dHkC7zs/hlTMgmWJ7nIi40pqCwm3rrsMCgCeAUTAVKIQ2bOATCLGg3+HXbYWdZdAAcEAiAALwAePR/daAHVOXDycdDj8Qw0MgIg4ncjjRETPAqhzZlObSChHdxkUEAyAiDBlXQAdDWf/UCsGQIScivMsgJ4tCQ8nOPxDLRgAEdI4C2CPINrfbHoTccW+//QYAyBiTsULuDaRQYwzZ+gJlzj7h57AdQARIwCmjw3j1OgINss13FkvoFrnsJDpUuJiksM/9AQGQMR4VgoKFmICTA5kMDGQRqHm4P5mAZslLv4x1Wx6A3GHbwTonRgAEeNZ72yZIBAMp5N4bWoMNcfD/E4Jc9vcANw0l4Wzf+hpDICIcWX/njkCQSYRx5WJEVw4PoTlQhX31wtw2Wo58tKWiwmH4//0NAZAxLhWCnjBMT1hWTg70o8zw33YrNRwb72Aks3hgai6nlpHjMM/tA8GQMTU0X7XTEsEE/0ZjJ9Lo2A7eLhZxHqRbZej5qLw3T/tjwEQMQcJgF0CwXAqiXedPgbb9ZvXCcpQHB0KvbTlYrzO8X/aHwMgYurq8H3zBYJ0PIbL48O4MDaI1WIVd9cLcDwmQVi9nFpDzPF1l0EBxQCIEAUFG4mOPFbcsjA13I9Tw33Yrti4t15AocZVpGHzEod/6DkYABHiN9cAdJIFwfG+NMZm0ijVG9cJVgu1jj4HdUefVcfx+rLuMijA2AoiQlyre9smCoDBZAKvnBrFT16cxEvjAxB2mwi0G+kNxMDhH3o2ngFESDcDYJdAkIrF8NLYEM6NDmK11LhOUHd5oAmal8DWD/R8DIAIcZHu6fPFLMHpoT6cHMpgp1LHvY0C8lVuNhIEAxz+oTYwACLEkeQLF4F1gwXBWF8KN8+Oo1J38PZWCUv5au8LoT3XU2uwHM7eoudjAETIYdYAdJIA6E8mcOPkKC5PDGMxX8aDjRJ8LijouZfA2T/0YgyACLE1B0CrZMzC+WODODvaj/WSjTvredicj94Tg7E6xpwV3WVQCDAAIkJBwVadWQPQSTGxcHIwg8mBNAp2HffWC9iu8DpBN91IrnL4h9rCAIgI30pCIaa7jGeyRDCSTuH16eOoOB7mtkvI7rAtdTdw+IfaxQCICO8ZbaCDRiDoT8RxbXIEF8cHsZyv4v5GER7bUnfEkGVjlMM/1CYGQES4Vm+ngHZC0ophZnQAZ0b6uX1lh7ycXoXFjd+oTQyAiGhnH4Cgiok83r7SdnB/g9tXHtYFxcVf1D4GQEQ4SOou4ch221Jz+8rDGYnVMOKs6S6DQoQBEBG61wB00pPbV64UqrjH7Stf6OUUh3/oYBgAERGkNQCdlLAsTI/0Y2q4D1uVGu5y+8pnOs/eP3RADIAIaOwDEP4hoOexRDDen8Hxc2kUbRcPNgvcvrLFsXgVI/U1NNZjE7WHARABviThq+CuAegkgWAoleD2lU+4kVyF1Hnwp4NhAESAZ0X73f9+uH3lO53n7B86BAZABIRxDUAntW5fuVOxcdew7SvHYlUMO+vg8A8dFAMgAtyQrALutkZb6jTemEmjXHfwwJDtK19OrXD4hw6FARABjoR3EVg3CICB5vaVVyd9ZHNlPNwsRfY6wTl/QXcJFFIMgAioR3wG0GG1bl85c2wQa8XobV85Hq9gyN0Eh3/oMBgAEVBXHAJ6kbg83r4yV220pc5FYPvKlzn7h46AARByjTUADIB2WRAcy6TwvrPjqDgO3t4M9/aVM/687hIoxCzdBYSW7wGe/oVIShLwArwPQFAJgP5EY/vKn7o0iUsTg7AkXO+kJ+NlDLpbusugEOMZwGE4FWDpu4DygWMvAcPT2krxLL77P6qkFXvH9pV31wuoOcFvN3Ejydk/dDQMgMOobqMx7UYBxWWtAeAyADqmdfvKfK2O+xvB3r5yhrN/6IgYAIeRGQPkbUApYGhKaylcA9B5lghGM43tK6uOh0cB3L7yZKKEAXcbnP1DR8EAOIxEBpj+QGMIyNL7T8g1AN0jEPQFdPvK64llDv/QkTEADkusxodmUdoHIMie3L7y3noBZY3bV3L4hzqBARByXATWW09uX/lgo4CNHm9feTpRxICb7+lzUjQxAELO5iIwLXa3r3z31Bhqrof57d5tX3k9uQLon4EcHUoBdgFIDQbirL6XGAAh5ksMrorxOqBGAkEm3sPtK5XCWY+LvzpGKeCHnwPyS0DfGPDaJ40KAQZAiHlWCtC4eMn1FOY3qnA8hZOjSQz3JbTVEgRPbl95b72Iot3ZttRTqRL63UJHH9NovgPkmtdTKltArQBkRrSW1EvmRF0E6V4DUKq5cJsbsGwWgjtfvtd2t698/7kJvP/cOCYHO/d7uh5f6thjEQArAYxdBCDA4EkgPaS7op7iGUCI6V4DkEk2hp9EAQMZtqN4kgAYSiXwaqe2r1QKZz3O/ukoEeDam422LjG9Z9Q6MABCzNE8BTSVsHBhsg+ur5CKm/XCOYj9tq+8t1E8cFvqs6kC+rxil6o0mAgQN3NXPQZAiNUDsAgsHhPEYzz4t+so21deSyxz9g91FAMgxLgGILye3L7y4WYRK8/bvlIpTLsc/qHOYgCEGPcBCL/H21cew5VJ75nbV55LF5DxSlpqpOhiAISUD4Hrx7kGIEKet33ltfgiEPwO1RQyDICQ8mJp42YsmOKp7SvX8phys7rLogjiOoCQ4kYw0be7feWfmh5ExivrLociiAEQUq6YOW3NTDzTo+5gAISUI5wBZA4GAHUHAyCkuA+AQXith7qEARBSXANgDsUzAOoSBkBIcR8Ak+jfgpKiiQEQQr4IHGV262WT8AyAuoUBEEK+mNe10Gz8XVN3MABCSPc+ANRrDADqDgZACLkW1wAYhWd71CUMgBDSvREM9ZZKFrFx6hIvBVPHMQBCyOYUUHOkd1BOzmF+4gQeXnkfnCTP/qhzGAAhxEVghkhvoRSf3/trLp3CrSvvQf74lMaiKEoYACHERWDRppSCyqyjFM9CnrgA7FgW7k+dw8KFd8Oz2MyXjoYBEDIKQF0xAKJKKQXVv4pKbOWpg/9jgvXBAdyefR8qg2M9rY+ihQEQMo020JwVEkUKCn7/EqrWelv3r8VjuH1+FqtnrvICMR0KAyBkPE4BjSgffn8WNWvrQN+lRLA4No57V2/CTg90qTaKKgZAyHARWPQoeHD651GTnUM/RjGVxK3L78L25LkOVkZRxwAIGYczgCJFiQe3fw51KRz5sTwRvH1yCo8uvgdejBeI6cUYACFT50Yw0WG5cPrfRl1KHXxQwVZ/H27Nvg/F0RMdfFyKIgZAyPAMICIsB3bfQziodOXh7VgMd89exNK5l+ELX+a0P/7PCBmuAo6AWB21vgdwUevyEwlWhkdwd/Yman0jXX4uCiMGQIhwDUAExGuoZu7DQ71nT1lOxHHr4g32E6KnMABCxLdSUPyVhVeigkr6AXy4PX9qX6SlnxCHEamBR5MQ4RTQEEuUUUk9hIKntYxGP6HX2U+IADAAQoVtoEMqVUA59RAKvu5KALT2E3oX+wkZjgEQIjwDCKFUDuXEIwRvY3fB+uAg7rCfkNEYACHCNtDhotLbKCXmdJfxXFX2EzIaAyBEHE4BDYXdds7l+MJzOnoGR2s/oXq6X3c51EMMgBCxFc8Agq69ds7BVEwlcevSu9lPyCAMgJBQUFwEFnAKCqp/ue12zkHkWuwnZBIGQEhwDUDQNdo5V61N3YV0wON+QqWRSd3FUBfxiBISnAEUXAoe3P6FI7VzDiI7FsOdmUtYmrnBfkIRxd9qSHANQECJD7d/DrbkdVfSJYKVkdFmP6Fh3cVQhzEAQoIBEEDiot7/sMPtnIOp0U/oZWycusjpohHCAAgJhwEQLJYDu7977ZyDqNFP6CQeXn4v+wlFBAMgJLgILEBiddT6HvagnXMw5TLpRj+hMfYTCjsGQEhwCmgwSLyGauYBPNi6S9HKsSzcP9PoJ+RbMd3l0CExAEJAQcHmPgD6JSoopx/Ah6O7koBo9BO6PXsTlYFjuouhQ2AAhIBvJaHAd1laBaSdcxBV4zHcvnANq2dmeYE4ZBgAIeBxDYBeAWvnHESNfkLHcZ/9hEKFARACrqR1l2CuwLZzDqZCs5/QzsSM7lKoDQyAEHCF4/86qPQ2yol53WWEjmsJHp46w35CIcAACAGuAeg9ldlAOb6gu4wQYz+hMGAAhADXAPSOUgp+3yoqseXQtXMOIjsWw92ZS1ieuQEl/PcMGgZACNgMgJ5QSkENLKNqrekuJVIUBMsjo7gz+wb7CQVMTwfoCj9axfa3FmAlE7CScUgqASsZg5WIQZIxWEkLkrBgJQArbkFifMegoFDnGoAe8OEPLEauo2eQ7PYTml5fxdjKfZ5fBUBPA0A8H16+0uZMagWJx2GlEpBkAlaqNTDikJQFKxGDlbAgCWn8GQckYqeZviTgcQ1AVyl48PqzEe7oGRy+COYmTyI3NIqzj36ARN3sFdW6BfgSvUC5HjzXA8pt9lwRwEolIal4ZM4yuAagy8SH2zePuhR1V2KUXCaN8pXXMbM0h+GtRd3lGCvAAXAICvBrdaBWP/xZRjIOKxUPzFmGa3ENQNeIi3r/I6M6egbJbj+hiZFxTD36ISyfq6x7LVoBcGDBP8vgPgBdYjmw+942tqNncDT6CRVnb+Lc3G30lbZ1F2QUwwPgEA5zlhGLwUonD3WW4bALaMd5VgVO34LxHT2DpBqP4c6Fazi9vYWJ7C1eIO4RBkDXCZTnwyvXDnWW4Y3UEJsdgpfq626ZEaaUQkVtY9tZxLq9gJpfxit9VzjPP2B8EWTHjiM/cBMzb/8QSZtDc93GAAiilrMMyc9hdHkJ9vXXUJy+BHBz7rb48FDy1rHpZLFeW0DVe+fBpGY7yKR4dhVEhVQSty6/hrOrCxhdZyuObmIAhIHjIPXHf4jkwkMUX3kD9SH2Xt+Pq+rIe8uNg351Ea56dt/+xdwKXpqc5llAQDX6CU3j+PBxnHn7B4h5ru6SIokBECKytYGhr34ezuVrKFx8BSqW0F2SdrYqI+c2hna27FWoNrt27tRy8LwpxGNcYxFcgs3+/sYF4uw9DOTWdRcUOQyAsFEKiTs/xrHsI1RfvYnK+BndFfWUUgpVlcO2k8V6PYucvXnoabnrxS2cGpnocIXUaXbMwt2ZyziZm8TJ+R9DFFtzdwoDIKSkXEbmW/8HqalzKFx/HV6EN+FojOdvYMvJYt3OouKW9m47ypqMpfwKTgyPw4rY6vEo2u0nVOi7iZm5W0hXuGq7ExgAISYQxBbnMLoWvYvEHhzkvWVs1bNYqy3C8esdfw4FIF8uYnRgqOOPTd1RSibYT6iDGABRsHeR+AFKr74Be3BMd0WH4qgKdrwlbNgL2LBXoFT3t2Ccyy1iZGCWB5IQ2e0nlB8axTT7CR0JAyBCZGsTg1/5PNKXrqFw6dXAXyRWSqGm8th2s9iws9i2N3rezM/xHFRqNfSn2XIjbHYyaZSuvI5zS48wtLWku5xQYgBEjQISd9/CscW5QF4k9uGj7G1gy81irZZFxX3chE1XJ9dsbgmXT5znlNAQciwL986cx+TwOE7P/Yj9hA6IARBRexeJz5xD4Zrei8QeHBTcVWw6C1ivLaLuB+uUvWCX4Lo+EnFOCQ0nwdrQEAqzN3F+7jYy7CfUNgZAhAkEsewcRlYWUb/xnp5eJK6rKnLuEjbrWWzUluG32TlJl9XCBs4cO6G7DDqCajyG2xeu4fT2Jiayt3k+1wYGgAHEdR9fJH7lDdhD3blIXFMFbDvN8fx6uBbtrBRXcXpkEpbFw0aYNfoJjSM/MMx+Qm1gABhEtjYx+NXPI32xeZE4frSLxD58lP0tbDtZrNUWUHYLHapUB8FOpYCxAe5ZGwXsJ9QeBoBpFJC49xaOLTZXEk9MH+jbPeWi6K82++1kYfvR6ac/t5PFsf6hyG0rair2E3oxBoChpFJB5tt/gNTUDArX3/vci8QOai3j+UvwVLDH8w/L8z2U7RoG0hndpVDHsJ/Q8zAADNZYSTyPkdUlODdeQ2H68t5F4ppfxI7XaL2wba9prrR3FnaWcPXkBU4JjZi9fkI7Ezi58Bb7CTUxAAjiukj+8f/F2PwD/Pi1KSx5WRTqO0YOhZTqZTiOh2SCL42oURAsjx5Dof8mzj16C6lqmK9ZdUY0GsdQR8j2Fh6V30LRyRl58N+1nDfnjMdEpWQCty69gs2TL7XZPDy6GAC0x7cA24rm+P5BrJc34PmmHxqizRPB3OQpvH35vXAT5u4MxwCgPV46ARj8zv8xwVZpR3cR1AM7mTRuXX0vCmOndZeiBQOA9jhJtkLYtZBbhOKFQiPUm/2EsudfhW+Z9RpgANAeJ8X/Drt8pVCschWpORr9hG7P3kR1wJw9t/mKpz12kv8dWi3kltreY5iiYbef0NrUFSN+83zF0x47yfH/VhWnCtvh6lHT+CLIHp/A/as3UU/16S6nqxgAtKcW7P1jtFjKreougTTZ7Se0M3FWdyldwwCgPdVE97dgDJutyjZcj/8uptrtJzR38TV4segtDmQA0J5q3IRRz4Pb5AYjhmv0E7o1exPl4XHdxXQUA4D2lBJcBLafbH6JU0IJdszCnXNXsHz2OlRE1sswAGhPOc4A2I9SQL5a1l0GBcBuP6G7V2/CzgzpLufIGAAEgG0gXmR+Z5FTQmnPbj+hrRMXQv2/ggFAANgG4kVs10bNrusugwLEE8GjE6dD3U+IAUAA2AaiHYu5VZ4F0FP2+gkdO6W7lANjABAAtoFox04tB9fllFB6Wt2ycG/6Quj6CfFVTwDYBqJdG6Ut3SVQYDX6Cd2ZvYlq/6juYtrCVz0BAOwEx//bsZRfgc8pofQclXgMt1+6jvWpq4EfMGQAEACgFs5rWD2nAOTLRd1lUMD5Ilg4Po77V2/CSWV0l/NMDAACANQSQX+vEhxzOU4JpfYUUkm8dfk9yI0Hs58QA4AAAJU4L262y/EcVGq27jIoJFxL8OD0NOZeenfg+gkxAAgA20AcVJZ7BdCBCDYHBnA7YP2EGAAEAKjE2Pf+IAp2iVNC6cBqAesnxAAg+BZQi/FgdlCrhQ3dJVAIBamfEAOA2AbikFaKq/B9DgPR4QShnxADgNgG4tAE2+W87iIoxHb7CT26/LqWfkIMAGIbiCOYzy1yrwA6su1MRks/Ib7yCXW2gTg0z/dQtjkllI5ut5/Q4rlXetZPiK98Qo1tII5kIUB7Bbi+g/nqQyxV5+ErXtgPH8Hq8HDP+gkFa1UCaWGzDcSRlOplOI6HZEL/yylbe4RtZwsCIC4JTKbD16KYHvcTOut0t40EzwAIVbaBOLLl/JruEgAAgsdnc8KZXaHmiyA/MNjV59D/loW0YxuIo1svb+DMsVOIWXoPumcy5xCXBOISx3jyhNZaKPgYAIQKN4PvAMFWaQcTQ8e0VhGTOKYyM1proPDgEBChFGcbiE5Y4JRQChkGgOHYBqJzfKVQrFZ0l0HUNgaA4bxUnG0gOmg+QFNCiV6EAWA4J8XLQJ1UdWuw647uMojawgAwnMNVwB23lF/jWQCFAl/9hquzD1DHbVW24XsMAAo+vvoNxzYQ3bFR2tZdAtELMQAMZyd0VxBN2fwSp4RS4DEADFdN8iDVDUoB+WpZdxlEz8UAMFyVbSC6hlNCKegYAC38Wg3VP/kT1LNZ3aX0TJltILrGdm3U7LruMoieiQHQovCFL6D8ne+g8KUvwd0wY8NvtoHormxuhWcBFFgMgBZ+pQL4jSERv1rVXE33sQ1E9+Vqebgu/40pmBgALQY//GHET55E+vp1JM6c0V1O17ENRG+sF7d0l0C0L/YBaJE4cQIjb76pu4yecVNxAByj7rblwgpOjozDYthSwPAMwGDcDL43FIB8uaS7DKKn8AhgMLaB6J25XJYXgylweAQwmM02ED3jeA4qNVt3GUTvwAAwWI1tIHoqm1viWQAFCgPAYDW2geipgl3ilFAKFAaAwSpsA9FzqwUzFhhSODAADGZCGwilFL6x+A187s7n8Cj/SHc5WCmuwvd55kXBwAAwmAkBsFpZxZ2dO8jZOXxl4Su6ywEg2C7ndRdBBIABYCzfAqqx6AdAf6IfUEBc4hhMDuouBwAwn1vkXgEUCFwJbKhGG4joXwMYSg7hzYtvYr2yjvPD53WXAwDwfA9lu4aBdEZ3KWQ4ngEYyk2ak/1jmTFcHbuKVDylu5Q989vcK4D0YwAYiquA9So7FdSd6A/BUbDxKGAoh32AtFvJr+kugQzHo4Chakm2gdBtvbwBz4/+dRgKLgaAodgGIggEW6Wc7iLIYAwAQ9USvAAZBAucEkoaMQAMVUlw6CEIfKVQrFZ0l0GGYgAYqswACIz5HU4JJT0YAIYqx1zdJVBT1a3Brju6yyADMQAMpMSMNhBhspRf5VkA9RwDwEBuKg5wg/JA2arswPcYANRbDAADuSlz2kCEyXppW3cJZBgGgIGcdEx3CbSPxfwSp4RSTzEADFTnZvCBpBSQr5Z1l0EGYQAYiG0ggotTQqmXGAAGstkGIrBs10bNrusugwzBADBQlW0gAi2bW+FZAPUEA8BAFQZAoOVqebguV2pT9zEADFQxYDP4sFsvbukugQzAADBQKc42EEG3XFiBzymh1GUMAMOwDUQ4KAC5clF3GRRxDADDsA1EeMznOCWUuosBYBi2gQgPx3NQqdm6y6AIYwAYxknxVx4m2dwSzwKoa3g0MEw9yV95mBTsEhyH12yoO3g0MIzNPkChs1bc1F0CRRQDwDC1pO4K6KBWiqvwfQ4DUef19Irg18tlTP+dv93Lp6QnuDHAY+yHTsyyIHy/Zhzl+xgdHe3a4/c0AOo+l7cTEQUF31IQERmKAUBEZCgGABGRoRgARESGYgAQERmKAUBEZCgGABGRoRgARESGYgAQERlK1AG2nRORDQDz3SuHiIg67KxSany/Gw4UAEREFB0cAiIiMhQDgIjIUAwAIiJDMQCIiAzFACAiMhQDgIjIUAwAIiJDMQCIiAzFACAiMtT/B90Xac9Z7L8mAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = subplot_for_map()\n", "plot_voronoi_polys_with_points_in_area(ax, boundary_shape, poly_shapes, pts, poly_to_pt_assignments)\n", "ax.set_title('Voronoi regions')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a graph from OSM within the boundaries of coverage area. Save the graph to geodataframe" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "G = ox.graph_from_polygon(boundary_shape, network_type='all_private')\n", "gdf_all_streets = ox.graph_to_gdfs(G, nodes=False, edges=True,node_geometry=False, fill_edge_geometry=True)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
osmidhighwayonewaylengthgeometryserviceaccessnamelanesest_widthbridgejunctiontunnelmaxspeedwidthuvkey
0816728678serviceFalse76.900LINESTRING (106.63699 -6.28336, 106.63630 -6.2...NaNNaNNaNNaNNaNNaNNaNNaNNaNNaN762796441878479709090
1841101956serviceTrue42.409LINESTRING (106.63699 -6.28336, 106.63684 -6.2...NaNNaNNaNNaNNaNNaNNaNNaNNaNNaN762796441878479404350
2841101955serviceFalse70.538LINESTRING (106.63485 -6.28342, 106.63491 -6.2...NaNNaNNaNNaNNaNNaNNaNNaNNaNNaN762796441978479709040
3[591697100, 695004245]serviceFalse305.493LINESTRING (106.63485 -6.28342, 106.63484 -6.2...NaNNaNNaNNaNNaNNaNNaNNaNNaNNaN762796441965283946350
4816728678serviceFalse39.691LINESTRING (106.63485 -6.28342, 106.63502 -6.2...NaNNaNNaNNaNNaNNaNNaNNaNNaNNaN762796441976279644130
\n", "
" ], "text/plain": [ " osmid highway oneway length \\\n", "0 816728678 service False 76.900 \n", "1 841101956 service True 42.409 \n", "2 841101955 service False 70.538 \n", "3 [591697100, 695004245] service False 305.493 \n", "4 816728678 service False 39.691 \n", "\n", " geometry service access name \\\n", "0 LINESTRING (106.63699 -6.28336, 106.63630 -6.2... NaN NaN NaN \n", "1 LINESTRING (106.63699 -6.28336, 106.63684 -6.2... NaN NaN NaN \n", "2 LINESTRING (106.63485 -6.28342, 106.63491 -6.2... NaN NaN NaN \n", "3 LINESTRING (106.63485 -6.28342, 106.63484 -6.2... NaN NaN NaN \n", "4 LINESTRING (106.63485 -6.28342, 106.63502 -6.2... NaN NaN NaN \n", "\n", " lanes est_width bridge junction tunnel maxspeed width u \\\n", "0 NaN NaN NaN NaN NaN NaN NaN 7627964418 \n", "1 NaN NaN NaN NaN NaN NaN NaN 7627964418 \n", "2 NaN NaN NaN NaN NaN NaN NaN 7627964419 \n", "3 NaN NaN NaN NaN NaN NaN NaN 7627964419 \n", "4 NaN NaN NaN NaN NaN NaN NaN 7627964419 \n", "\n", " v key \n", "0 7847970909 0 \n", "1 7847940435 0 \n", "2 7847970904 0 \n", "3 6528394635 0 \n", "4 7627964413 0 " ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gdf_all_streets.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create new dataframe to collect street networks within each voronoi regions" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "gdf_streets_by_region = gpd.GeoDataFrame()\n", "for x in range(len(poly_shapes)):\n", " gdf_streets = gpd.GeoDataFrame()\n", " gdf_streets['geometry'] = gdf_all_streets.intersection(poly_shapes[x])\n", " gdf_streets['voronoi_region'] = x\n", " gdf_streets = gdf_streets[gdf_streets['geometry'].astype(str) != 'LINESTRING EMPTY']\n", " gdf_streets_by_region = gdf_streets_by_region.append(gdf_streets)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
geometryvoronoi_region
34LINESTRING (106.65206 -6.28108, 106.65205 -6.2...0
1126LINESTRING (106.64790 -6.28385, 106.64835 -6.2...0
1127LINESTRING (106.64790 -6.28385, 106.64761 -6.2...0
1128LINESTRING (106.64790 -6.28385, 106.64811 -6.2...0
1132LINESTRING (106.64854 -6.28395, 106.64838 -6.2...0
\n", "
" ], "text/plain": [ " geometry voronoi_region\n", "34 LINESTRING (106.65206 -6.28108, 106.65205 -6.2... 0\n", "1126 LINESTRING (106.64790 -6.28385, 106.64835 -6.2... 0\n", "1127 LINESTRING (106.64790 -6.28385, 106.64761 -6.2... 0\n", "1128 LINESTRING (106.64790 -6.28385, 106.64811 -6.2... 0\n", "1132 LINESTRING (106.64854 -6.28395, 106.64838 -6.2... 0" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gdf_streets_by_region.head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visualize the voronoi diagram using folium" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = folium.Map([-6.304029, 106.645874], zoom_start=13, tiles='cartodbpositron')\n", "#draw the voronoi diagram within coverage area\n", "for x in range(len(poly_shapes)):\n", " folium.GeoJson(poly_shapes[x]).add_to(m)\n", " \n", "#draw the data points\n", "points = np.array([[-6.305286, 106.644085], [-6.301309, 106.653261], [-6.284774, 106.637751], [-6.284598, 106.665062], [-6.283521, 106.627582],[-6.276593, 106.641365], [-6.303643, 106.625972]])\n", "locs = points\n", "for location in locs:\n", " folium.CircleMarker(location=location, \n", " color = \"#4925a2\", radius=0.01).add_to(m)\n", "\n", "\n", "#draw street networks for each voronoi region\n", "folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 0], style_function=lambda x: {'color': '#30e861', 'weight':'1'}).add_to(m)\n", "folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 1], style_function=lambda x: {'color': '#ca1890', 'weight':'1'}).add_to(m)\n", "folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 2], style_function=lambda x: {'color': '#d5e849', 'weight':'1'}).add_to(m)\n", "folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 3], style_function=lambda x: {'color': '#1b12e4', 'weight':'1'}).add_to(m)\n", "folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 4], style_function=lambda x: {'color': '#ee3d07', 'weight':'1'}).add_to(m)\n", "folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 5], style_function=lambda x: {'color': '#eead07', 'weight':'1'}).add_to(m)\n", "folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 6], style_function=lambda x: {'color': '#059b3a', 'weight':'1'}).add_to(m)\n", "\n", "\n", "\n", " \n", "m" ] }, { "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.7.4" } }, "nbformat": 4, "nbformat_minor": 2 }