generate spatial weights for blocks.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-09-24T18:46:49.023169Z",
     "start_time": "2023-09-24T18:46:49.005184Z"
    }
   },
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*- \n",
    "from __future__ import unicode_literals\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import json\n",
    "import itertools as itt\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "\n",
    "import geopandas as gpd\n",
    "import shapely\n",
    "import fiona\n",
    "import CoordinatesConverter as convert_xy  #the .py file for bd09 conversioin\n",
    "\n",
    "%matplotlib inline\n",
    "#plt.style.use('ggplot')\n",
    "plt.rcParams['figure.figsize']=[18,8]\n",
    "plt.rcParams['axes.titlesize'] = 'medium'\n",
    "plt.rcParams['font.sans-serif'] = 'SimHei'\n",
    "plt.rcParams['axes.unicode_minus'] = False\n",
    "plt.rcParams['legend.frameon'] = False\n",
    "\n",
    "\n",
    "WGS84 = 'EPSG:4326'\n",
    "#UTM_bj = ccrs.UTM(50)\n",
    "UTM_bj = 'EPSG:32650'\n",
    "UTM_nj = 'EPSG:32650'\n",
    "\n",
    "\n",
    "## set the projection (planar) crs of this area\n",
    "CRS_planar = UTM_nj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# loading data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## overall boundary\n",
    "the maximum extent for  blocks generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-09-24T18:01:32.125296Z",
     "start_time": "2023-09-24T18:01:32.080287Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "boundary = gpd.read_file(r'D:\\Research\\Some Basic Data\\Nanjing\\cityborder_2035').to_crs(CRS_planar)\n",
    "boundary = boundary.geometry.explode(index_parts=True).reset_index(drop=True).reset_index()\n",
    "boundary['geometry'] = [shapely.geometry.Polygon(i) for i in boundary['geometry']]\n",
    "boundary['geometry'] = boundary.buffer(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## blocks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-26T18:59:26.412986Z",
     "start_time": "2022-04-26T18:59:25.958204Z"
    }
   },
   "outputs": [],
   "source": [
    "blocks = gpd.read_file('../data/blocks').to_crs(CRS_planar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# voronoi diagram of block polygons (not used)\n",
    "- voronoi diagram of polygons or lines can be approximated with voronoi of points on polygon boundaries   \n",
    "but \n",
    "- how to deal with crossings?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-26T18:59:28.524385Z",
     "start_time": "2022-04-26T18:59:28.498388Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def voronoi_diagram_polygons(data,densify_max_span=5, envelope=None,tolerance=0):\n",
    "    # ============================densify the polygon boundary vertices========================\n",
    "    #### retrieve the polygons exterior as well as interior\n",
    "    geom = data.boundary.explode(index_parts=False)    \n",
    "\n",
    "    #### explode into line segments\n",
    "    def explode_line(temp):\n",
    "        t = np.array(temp.coords)\n",
    "        return [shapely.geometry.LineString([i,j]) for i,j in zip(t[:-1],t[1:])]\n",
    "    geom = [j for i in geom for j in explode_line(i.simplify(tolerance=0.5))]\n",
    "    geom = gpd.GeoSeries(geom,crs=data.crs).reset_index(name='geometry')\n",
    "\n",
    "    #### divide long line segments\n",
    "    geom['length'] = geom.length\n",
    "    geom['divide'] = (geom['length']/densify_max_span).apply(np.ceil).astype(int)\n",
    "    geom['geometry'] = [shapely.geometry.LineString([line.interpolate(p, normalized=True) for p in np.linspace(0,1,divide+1)]) \\\n",
    "                        for line,divide in zip(geom['geometry'],geom['divide'])]\n",
    "\n",
    "    #### merge into whole\n",
    "    geom = shapely.ops.linemerge(list(geom['geometry']))\n",
    "    \n",
    "    \n",
    "    \n",
    "    # ====================direct voronoi of the densified polygon boundaries=====================\n",
    "    temp = shapely.ops.voronoi_diagram(geom, tolerance=tolerance)\n",
    "    temp = gpd.GeoDataFrame({'geometry':list(temp.geoms)},crs=data.crs)\n",
    "    temp['geometry'] = temp['geometry'].buffer(0) ### the returned geometry may be invalid! so first buffer(0) to make them valid\n",
    "    \n",
    "    \n",
    "    \n",
    "    # =======================join voronoi result with original polygon and dissolve============= \n",
    "    s = temp.sjoin(data,how='inner').dissolve('index_right')\n",
    "    \n",
    "    \n",
    "    # ====================================clip by envelope======================================\n",
    "    if envelope is not None:\n",
    "        s = s.clip(envelope)\n",
    "    \n",
    "    s = s['geometry'].reset_index(drop=True).reset_index()\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 285,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-25T00:45:02.589304Z",
     "start_time": "2022-04-25T00:45:02.475439Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\envs\\pycaret\\lib\\site-packages\\geopandas\\geodataframe.py:1327: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  result = super().__getitem__(key)\n",
      "C:\\ProgramData\\Anaconda3\\envs\\pycaret\\lib\\site-packages\\geopandas\\geodataframe.py:1327: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  result = super().__getitem__(key)\n"
     ]
    }
   ],
   "source": [
    "s = blocks.sjoin(boundary)\n",
    "\n",
    "blocks1 = blocks[s['index']==0].copy()\n",
    "blocks2 = blocks[s['index']==1].copy()\n",
    "\n",
    "\n",
    "\n",
    "blocks1_vornoi = voronoi_diagram_polygons(blocks1,envelope=boundary['geometry'][0],tolerance=0.1)\n",
    "blocks2_vornoi = voronoi_diagram_polygons(blocks2,envelope=boundary['geometry'][1],tolerance=0.1)\n",
    "\n",
    "blocks_voronoi = pd.concat([blocks1_vornoi,blocks2_vornoi],ignore_index=True)\n",
    "#blocks_voronoi.to_file('../data/blocks_voronoi')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# skeletonize, e.g. medial axis transform (MAT), chordal axis\n",
    "\n",
    "[This package](https://github.com/NRCan/geo_sim_processing/blob/7894e6913d6ac51781e4656baf5bca17925efaec/chordal_axis_algorithm.py#L504) seems relevant. It employs chordal axis, and ___deals with crossing correction___. But it is in QGIS.  \n",
    "The result of this package seems __reasonable__.\n",
    "## generate the separation polygon for processing in QGIS\n",
    "Note:\n",
    "- fist check the geometry of generated blocks, or tessellation in QGIS would report error\n",
    "- better to buffer out the overall boundary, making generated chordal axis longer than the actual overall boundary, easy for cut\n",
    "    - do not buffer out too much (about 500m), or chordal axis calculation in QGIS would report error\n",
    "- should densify polygon boundary to produce finer result\n",
    "- better to export separation polygon separately into QGIS, e.g. \n",
    "    - first QGIS processes the main separation polygon (only one big polygon), then other separtion polygons (perhaps several polygons)\n",
    "    - finally combine the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-26T18:59:38.697997Z",
     "start_time": "2022-04-26T18:59:38.682006Z"
    },
    "code_folding": [
     0,
     5
    ]
   },
   "outputs": [],
   "source": [
    "def densify_polygon_boundary(polygon, densify_max_span=5):\n",
    "    densified = shapely.geometry.Polygon(shell = densify_line(polygon.exterior, densify_max_span=densify_max_span),\n",
    "                                         holes = [densify_line(i, densify_max_span=densify_max_span) for i in polygon.interiors])\n",
    "    return densified\n",
    "\n",
    "def densify_line(line, densify_max_span=5):    \n",
    "\n",
    "    line = line.simplify(tolerance=0.5)  ##simplify to remove vertices that are extremely close\n",
    "    t = np.array(line.coords)\n",
    "\n",
    "    vertices = []\n",
    "    for i,j in zip(t[:-1],t[1:]):\n",
    "        segment = shapely.geometry.LineString([i,j])\n",
    "        divide = int(np.ceil(segment.length/densify_max_span))\n",
    "        vertices += [segment.interpolate(p, normalized=True) for p in np.linspace(0,1,divide,endpoint=False)] ##middle segments including start without end\n",
    "    vertices.append(shapely.geometry.Point(t[-1]))  ##add the last point\n",
    "\n",
    "    densified = shapely.geometry.LineString(vertices)\n",
    "    return densified"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-26T19:01:08.393117Z",
     "start_time": "2022-04-26T18:59:58.251678Z"
    }
   },
   "outputs": [],
   "source": [
    "blocks['geometry'] = [densify_polygon_boundary(i) for i in blocks['geometry']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-26T19:01:08.693944Z",
     "start_time": "2022-04-26T19:01:08.397115Z"
    }
   },
   "outputs": [],
   "source": [
    "boundary1 = boundary[:1]\n",
    "boundary2 = boundary[1:]\n",
    "\n",
    "s = blocks.sjoin(boundary1)\n",
    "blocks1 = blocks[blocks.index.isin(s.index)].copy()\n",
    "\n",
    "s = blocks.sjoin(boundary2)\n",
    "blocks2 = blocks[blocks.index.isin(s.index)].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-27T22:17:55.257273Z",
     "start_time": "2022-04-27T22:17:55.237286Z"
    }
   },
   "outputs": [],
   "source": [
    "bg1 = boundary1.buffer(500,join_style=2)\n",
    "bg2 = boundary2.buffer(500,join_style=2) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-27T22:41:53.109150Z",
     "start_time": "2022-04-27T22:41:44.550843Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sep = densify_polygon_boundary(bg1.unary_union).difference(blocks1.unary_union)\n",
    "sep = gpd.GeoDataFrame({'geometry':list(sep.geoms)},crs=blocks.crs)\n",
    "sep['geometry'] = sep['geometry'].buffer(0)\n",
    "sep['area'] = sep.area\n",
    "sep = sep.sort_values('area',ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-27T22:42:37.916623Z",
     "start_time": "2022-04-27T22:42:37.862669Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\envs\\pycaret\\lib\\site-packages\\geopandas\\io\\file.py:362: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.\n",
      "  pd.Int64Index,\n"
     ]
    }
   ],
   "source": [
    "sep.iloc[1:].to_file('../data/sep2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## process QGIS result\n",
    "- polygonize the chordal axis to form polygons\n",
    "- cut polygons with boundary, as boundary is previously buffered out to elongate chordal axis\n",
    "- spatial join the polygons with blocks to get the polygons that represent a real block (not some sliver polygons generated from boundary cut), to ensure that:\n",
    "    - every block gets __one and only one__ polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-28T21:13:02.311821Z",
     "start_time": "2022-04-28T21:13:01.492214Z"
    }
   },
   "outputs": [],
   "source": [
    "upper_axis1 = gpd.read_file('../data/chordal axis/axis2.shp')\n",
    "upper_axis2 = gpd.read_file('../data/chordal axis/axis3.shp')\n",
    "upper_axis = pd.concat([upper_axis1,upper_axis2],ignore_index=True)\n",
    "\n",
    "upper = shapely.ops.polygonize(upper_axis.unary_union)\n",
    "upper = gpd.GeoDataFrame({'geometry':list(upper.geoms)},crs=upper_axis.crs)\n",
    "upper = upper.overlay(boundary2)\n",
    "\n",
    "s = upper.sjoin(blocks2,how='left')\n",
    "upper = s[['blockID','geometry']].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-28T17:31:14.344924Z",
     "start_time": "2022-04-28T17:31:13.730985Z"
    }
   },
   "outputs": [],
   "source": [
    "lower_axis1 = gpd.read_file('../data/chordal axis/lower_axis1.shp')\n",
    "lower_axis2 = gpd.read_file('../data/chordal axis/lower_axis2.shp')\n",
    "lower_axis = pd.concat([lower_axis1,lower_axis2],ignore_index=True)\n",
    "\n",
    "lower = shapely.ops.polygonize(lower_axis.unary_union)\n",
    "lower = gpd.GeoDataFrame({'geometry':list(lower)},crs=lower_axis.crs)\n",
    "lower = lower.overlay(boundary1)\n",
    "\n",
    "s = lower.sjoin(blocks1,how='left')\n",
    "lower = s[['blockID','geometry']].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-28T21:33:20.758911Z",
     "start_time": "2022-04-28T21:33:13.117029Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\envs\\pycaret\\lib\\site-packages\\geopandas\\io\\file.py:362: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.\n",
      "  pd.Int64Index,\n"
     ]
    }
   ],
   "source": [
    "adjoining_blocks = pd.concat([upper,lower],ignore_index=True).sort_values('blockID')\n",
    "adjoining_blocks.to_file('../data/adjoining_blocks')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:pycaret]",
   "language": "python",
   "name": "conda-env-pycaret-py"
  },
  "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.16"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "288px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}