Test for building PySAL weights objrect from shapely polygons.

weights_from_shapely.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing the construction of PySAL weights object from list of shapely objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Author: David C. Folch <[email protected]>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on Serge Rey's weights_from_geojson.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pysal as ps\n",
    "from shapely.geometry import MultiPoint\n",
    "import numpy as np\n",
    "from scipy.spatial import Voronoi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class PolygonCollection:\n",
    "    def __init__(self, polygons, bbox=None):\n",
    "        \"\"\"\n",
    "\n",
    "        Parameters\n",
    "        ==========\n",
    "        polygons: dict\n",
    "                  key is polygon Id, value is PySAL Polygon object\n",
    "        bbox: list (optional)\n",
    "              [left, lower, right, upper]\n",
    "\n",
    "        Notes\n",
    "        =====\n",
    "        bbox is supported in geojson specification at both the feature and feature collection level. However, not all geojson writers generate the bbox at the feature collection level. \n",
    "        In those cases, the bbox property will be set on initial access.\n",
    "\n",
    "        \"\"\"\n",
    "              \n",
    "        self.type=ps.cg.Polygon\n",
    "        self.n = len(polygons)\n",
    "        self.polygons = polygons\n",
    "        if bbox is None:\n",
    "            self._bbox = None\n",
    "        else:\n",
    "            self._bbox = bbox\n",
    "            \n",
    "    @property\n",
    "    def bbox(self):\n",
    "        bboxes = np.array([self.polygons[p].bbox for p in self.polygons])\n",
    "        mins = bboxes.min(axis=0)\n",
    "        maxs = bboxes.max(axis=0)\n",
    "        self._bbox = [ mins[0], mins[1], maxs[2], maxs[3] ]\n",
    "        return self._bbox\n",
    "        \n",
    "    \n",
    "    def __getitem__(self, index):\n",
    "        return self.polygons[index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some code to get finite voronoi polygons (ignores the vertices at infinity, which works for my use case). Basically it's a way to gin up some random connected shapely objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_voronoi(points):\n",
    "    vor = Voronoi(points)\n",
    "    vor_polys = []\n",
    "    for region in vor.regions:\n",
    "        if region:\n",
    "            if -1 in region:\n",
    "                region = region[:] # so we don't change the original voronoi output\n",
    "                region.remove(-1)  # there can only be one -1 in the index list\n",
    "            vor_polys.append(MultiPoint(vor.vertices[region]).convex_hull)\n",
    "    return vor_polys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def shapely2w(polys, wttype):\n",
    "    '''\n",
    "    Parameters:\n",
    "    ----------\n",
    "    polys:  list\n",
    "            List of shapely polygons\n",
    "    wttype: string\n",
    "            \"rook\" or \"queen\"\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    pysal weights object\n",
    "\n",
    "    '''\n",
    "    if wttype.upper() == 'QUEEN':\n",
    "        wttype = 1\n",
    "    elif wttype.upper() == 'ROOK':\n",
    "        wttype = 2\n",
    "    else:\n",
    "        raise Exception, \"unknown weight type\"\n",
    "\n",
    "    polys_pysal = []\n",
    "    poly_ids = []\n",
    "    poly_index = 0\n",
    "    for feature in polys:\n",
    "        polys_pysal.append(ps.cg.asShape(feature))\n",
    "        poly_ids.append(poly_index)\n",
    "        poly_index +=1\n",
    "    pcollection = PolygonCollection(dict(zip(poly_ids, polys_pysal)))\n",
    "    neighbors = ps.weights.Contiguity.ContiguityWeightsPolygons(pcollection, wttype=wttype).w\n",
    "    return ps.W(neighbors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.random.seed(999)\n",
    "pnts = np.random.random((100, 2))\n",
    "# get some random shapely polygons\n",
    "vor_polys_shapely = get_voronoi(pnts)\n",
    "# build a W object from the polygons\n",
    "w = shapely2w(vor_polys_shapely, 'queen')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}