notebooks/window.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Point Pattern Windows\n",
    "\n",
    "**Author: Serge Rey <sjsrey@gmail.com>**\n",
    "\n",
    "## Introduction\n",
    "Windows play several important roles in the analysis of planar point patterns. As we saw in the [introductory notebook](pointpattern.ipynb), the area of the window can be used to develop estimates of the intensity of the point pattern. A window also defines the domain for the point pattern and can support corrections for so-called edge effects in the statistical analysis of point patterns. However, there are different ways to define a window for a point pattern.\n",
    "\n",
    "This notebook provides an overview of how to work with windows and covers the following:\n",
    "\n",
    "* [Creating a window](#Creating-a-Window)\n",
    "* [Window attributes](#Window-Attributes)\n",
    "* [Window methods](#Window-Methods)\n",
    "* [Multi-part windows](#Multi-part-Windows)\n",
    "* [Windows and point pattern intensity revisited](#Windows-and-point-pattern-intensity-revisited)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a Window\n",
    "\n",
    "We will first continue on with an example from the [introductory notebook](pointpattern.ipynb). Recall this uses  200 randomly distributed points within the counties of Virginia. Coordinates are for UTM zone 17 N."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import libpysal as ps\n",
    "import numpy as np\n",
    "from pointpats import PointPattern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Point Pattern\n",
      "200 points\n",
      "Bounding rectangle [(273959.664381352,4049220.903414295), (972595.9895779632,4359604.85977962)]\n",
      "Area of window: 216845506675.0557\n",
      "Intensity estimate for window: 9.223156295311261e-10\n",
      "               x             y\n",
      "0  865322.486181  4.150317e+06\n",
      "1  774479.213103  4.258993e+06\n",
      "2  308048.692232  4.054700e+06\n",
      "3  670711.529980  4.258864e+06\n",
      "4  666254.475614  4.256514e+06\n"
     ]
    }
   ],
   "source": [
    "f = ps.examples.get_path('vautm17n_points.shp')\n",
    "fo = ps.io.open(f)\n",
    "pp_va = PointPattern(np.asarray([pnt for pnt in fo]))\n",
    "fo.close()\n",
    "pp_va.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the summary method we see that the **Bounding Rectangle** is reported along with the **Area of the window** for the point pattern. Two things to note here. \n",
    "\n",
    "First, the only argument we passed in to the `PointPattern`s constructor was the array of coordinates for the 200 points. In this case PySAL finds the [minimum bounding box](https://en.wikipedia.org/wiki/Minimum_bounding_rectangle) for the point pattern and uses this as the window.\n",
    "\n",
    "The second thing to note is that the area of the window in this case is simply the area of the bounding rectangle. Because we are using projected coordinates (UTM) the unit of measure for the area is in square meters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Window Attributes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The window is an attribute of the `PointPattern`. It is also an object with its own attributes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "216845506675.0557"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_va.window.area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[273959.664381352, 4049220.903414295, 972595.9895779632, 4359604.85977962]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_va.window.bbox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bounding box is given in left, bottom, right, top ordering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(623277.8269796579, 4204412.881596957)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_va.window.centroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[(273959.664381352, 4049220.903414295),\n",
       "  (273959.664381352, 4359604.85977962),\n",
       "  (972595.9895779632, 4359604.85977962),\n",
       "  (972595.9895779632, 4049220.903414295),\n",
       "  (273959.664381352, 4049220.903414295)]]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_va.window.parts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `parts` attribute for the `window` is a list of polygons. In this case the window has only a single part and it is a rectangular polygon with vertices listed clockwise in closed cartographic form."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Window Methods\n",
    "\n",
    "A window has several basic geometric operations that are heavily used in some of the other modules in the the `Point` package. Most of this is done under the hood and the user typically doesn't see this. However, there can be times when direct access to these method can be handy. Let's explore.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The window supports basic point containment checks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_va.window.contains_point((623277.82697965798, 4204412.8815969583))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This also applies to sequences of points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pnts = ((-623277.82697965798, 4204412.8815969583),\n",
    "        (623277.82697965798, 4204412.8815969583),\n",
    "        (1000.01, 200.9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([ 623277.82697966, 4204412.88159696])]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pnts_in = pp_va.window.filter_contained(pnts)\n",
    "pnts_in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Multi-part Windows\n",
    "\n",
    "Thus far our window was a simple bounding box. There many instances when the relevant containing geometry for a point pattern is more complex. Examples include multi-part polygons and polygons with holes.\n",
    "\n",
    "Here we construct such a window, one with two parts and one hole."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "parts = [[(0.0, 0.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0)],\n",
    "         [(11.,11.), (11.,20.), (20.,20.), (20.,11.)]]\n",
    "holes = [[(3.0,3.0), (6.0, 3.0), (6.0, 6.0), (3.0, 6.0)]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will plot this using matplotlib to get a better understanding of the challenges that this type of window presents for statistical analysis of the associated point pattern."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADKlJREFUeJzt3H+o3fV9x/Hna6bdH9aylKYutYpObF1aWNZd7IqjWLq2KozUgUwHJWOF+IdChf1R1/6h/wxkzDoYnTTFoIyqK7TO0Lm2Kh3ZoOt6I0GjQQzO2WhIrrihf7Wo7/1xv2G38SY3N+d876/38wGHc873/Ph8vnw5z5x87/d8U1VIknr5tdWegCRp5Rl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIYmjn+SC5P8OMmhJM8k+fKw/H1JHkvy/HC9efLpSpKmIZMe559kK7C1qp5Mch6wH/gC8GfAa1V1Z5LbgM1V9ZVJJyxJmtzE3/yr6mhVPTncfgM4BFwA7ADuH552P/P/IEiS1oCJv/n/ypslFwP7gI8BL1XVbyx47H+q6h27fpLsAnYBnHvuub93+eWXT20+ktTB/v37X62qLct5zaZpDZ7kPcB3gVur6vUkZ/S6qtoN7AaYmZmp2dnZaU1JklpI8t/Lfc1UjvZJ8i7mw//tqvresPjY8PeAE38XOD6NsSRJk5vG0T4B7gUOVdXXFzy0F9g53N4JPDLpWJKk6ZjGbp8rgS8CTyc5MCz7KnAn8J0kXwJeAq6fwliSpCmYOP5V9e/AqXbwf2bS95ckTZ+/8JWkhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqSHjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpoanEP8meJMeTHFyw7I4kLyc5MFyuncZYkqTJTeub/33A1Yssv7uqtg+XR6c0liRpQlOJf1XtA16bxntJksY39j7/W5I8NewW2jzyWJKkMzRm/O8BLgW2A0eBuxZ7UpJdSWaTzM7NzY04HUnSCaPFv6qOVdVbVfU28C3gilM8b3dVzVTVzJYtW8aajiRpgdHin2TrgrvXAQdP9VxJ0sraNI03SfIgcBXw/iRHgNuBq5JsBwp4EbhpGmNJkiY3lfhX1Y2LLL53Gu8tSZo+f+ErSQ0Zf0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqSHjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkPGX5Iamkr8k+xJcjzJwQXL3pfksSTPD9ebpzGWJGly0/rmfx9w9UnLbgOeqKrLgCeG+5KkNWDTNN6kqvYlufikxTuAq4bb9wP/CnxlGuNp7Xngpy/xyIGXV3saOks7tl/An37iotWehlbQmPv8z6+qowDD9QcWe1KSXUlmk8zOzc2NOB2N6ZEDL/Ps0ddXexo6C88efd1/uBuayjf/SVTVbmA3wMzMTK3ydDSBbVvfyz/e9MnVnoaW6U+++ZPVnoJWwZjf/I8l2QowXB8fcSxJ0jKMGf+9wM7h9k7gkRHHkiQtw7QO9XwQ+AnwkSRHknwJuBP4bJLngc8O9yVJa8C0jva58RQPfWYa7y9Jmi5/4StJDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqSHjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ5vGHiDJi8AbwFvAm1U1M/aYkqTTGz3+g09X1asrNJYkaQnu9pGkhlYi/gX8KMn+JLtOfjDJriSzSWbn5uZWYDqSpJWI/5VV9XHgGuDmJJ9a+GBV7a6qmaqa2bJlywpMR5I0evyr6pXh+jjwMHDF2GNKkk5v1PgnOTfJeSduA58DDo45piRpaWMf7XM+8HCSE2M9UFU/GHlMSdISRo1/Vb0A/M6YY0iSls9DPSWpIeMvSQ0Zf0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqSHjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqaPT4J7k6yXNJDie5bezxJElLGzX+Sc4BvgFcA2wDbkyybcwxJUlLG/ub/xXA4ap6oap+CTwE7Bh5TEnSEjaN/P4XAD9fcP8I8ImFT0iyC9gFcNFFF408HY1l2wffu9pT0Fly2/U0dvyzyLL6lTtVu4HdADMzM7XI87UO3P5HH13tKegsue16Gnu3zxHgwgX3PwS8MvKYkqQljB3/nwGXJbkkybuBG4C9I48pSVrCqLt9qurNJLcAPwTOAfZU1TNjjilJWtrY+/ypqkeBR8ceR5J05vyFryQ1ZPwlqSHjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqSHjL0kNjRb/JHckeTnJgeFy7VhjSZKWZ9PI7393Vf3NyGNIkpbJ3T6S1NDY8b8lyVNJ9iTZPPJYkqQzNFH8kzye5OAilx3APcClwHbgKHDXKd5jV5LZJLNzc3OTTEeSdIZSVeMPklwMfL+qPna6583MzNTs7Ozo85GkjSTJ/qqaWc5rxjzaZ+uCu9cBB8caS5K0PGMe7fPXSbYDBbwI3DTiWJKkZRgt/lX1xbHeW5I0GQ/1lKSGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqSHjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNWT8Jakh4y9JDU0U/yTXJ3kmydtJZk567C+THE7yXJLPTzZNSdI0bZrw9QeBPwa+uXBhkm3ADcBHgQ8Cjyf5cFW9NeF4kqQpmOibf1UdqqrnFnloB/BQVf2iqv4LOAxcMclYkqTpmfSb/6lcAPzHgvtHhmXvkGQXsGu4+4skB0ea01rwfuDV1Z7EiFy/9Wsjrxts/PX7yHJfsGT8kzwO/OYiD32tqh451csWWVaLPbGqdgO7h7Fmq2pmsedtBK7f+raR128jrxv0WL/lvmbJ+FfVH57FXI4AFy64/yHglbN4H0nSCMY61HMvcEOSX09yCXAZ8J8jjSVJWqZJD/W8LskR4JPAPyf5IUBVPQN8B3gW+AFw8xke6bN7kvmsA67f+raR128jrxu4fu+QqkV3xUuSNjB/4StJDRl/SWpoTcS/02kiktyR5OUkB4bLtas9p0kluXrYPoeT3Lba85m2JC8meXrYXss+pG6tSbInyfGFv6lJ8r4kjyV5frjevJpznMQp1m9DfO6SXJjkx0kODc388rB82dtvTcSf/z9NxL6FC086TcTVwN8nOWflpzd1d1fV9uHy6GpPZhLD9vgGcA2wDbhx2G4bzaeH7bURjhW/j/nP00K3AU9U1WXAE8P99eo+3rl+sDE+d28Cf1FVvw38PnDz8Hlb9vZbE/H3NBHr2hXA4ap6oap+CTzE/HbTGlVV+4DXTlq8A7h/uH0/8IUVndQUnWL9NoSqOlpVTw633wAOMX/2hGVvvzUR/9O4APj5gvunPE3EOnNLkqeG/56u2/9eDzbqNlqogB8l2T+cjmQjOr+qjsJ8YIAPrPJ8xrCRPnckuRj4XeCnnMX2W7H4J3k8ycFFLqf7lnjGp4lYS5ZY13uAS4HtwFHgrlWd7OTW5TZapiur6uPM79q6OcmnVntCWrYN9blL8h7gu8CtVfX62bzHWCd2e4dOp4k403VN8i3g+yNPZ2zrchstR1W9MlwfT/Iw87u69p3+VevOsSRbq+pokq3A8dWe0DRV1bETt9f75y7Ju5gP/7er6nvD4mVvv7W+22fDnSZi2DAnXMf8H7vXs58BlyW5JMm7mf8D/d5VntPUJDk3yXknbgOfY/1vs8XsBXYOt3cCpzpp47q0UT53SQLcCxyqqq8veGjZ229N/MI3yXXA3wFbgP8FDlTV54fHvgb8OfN/5b61qv5l1SY6BUn+gfn/ehbwInDTiX1169Vw2NzfAucAe6rqr1Z5SlOT5LeAh4e7m4AH1vv6JXkQuIr50xwfA24H/on5U7JcBLwEXF9V6/KPpqdYv6vYAJ+7JH8A/BvwNPD2sPirzO/3X9b2WxPxlyStrLW+20eSNALjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhv4PqNMWcLBcafgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "p0 = np.asarray(parts[0])\n",
    "plt.plot(p0[:,0], p0[:,1])\n",
    "plt.xlim(-10,20)\n",
    "t = plt.ylim(-10,20) # silence the output of ylim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not, quite what we wanted, as the first part of our multi-part polygon is a ring, but it was not encoded in closed cartographic form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0.],\n",
       "       [ 0., 10.],\n",
       "       [10., 10.],\n",
       "       [10.,  0.]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can fix this with a helper function from the `window` module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0.0, 0.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0)]\n",
      "[(0.0, 0.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0), (0.0, 0.0)]\n"
     ]
    }
   ],
   "source": [
    "from pointpats.window import to_ccf\n",
    "print(parts[0])\n",
    "print(to_ccf(parts[0])) #get closed ring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADL9JREFUeJzt3F2oZfV9h/HnWyfphVE6komdGoNWbOw00Gl6sA2WYEiTqL2YWJBoIUxpYLxQiNCL2OQi3hSk1FgoqTjBQSlRG0isktokKim2YG3OyKCjgyhmakYH54gtepWi/npx1rQn45mXPXuv8/Z7PrDZe6/98v8vFvuZPeusvVJVSJJ6+aXVnoAkaeUZf0lqyPhLUkPGX5IaMv6S1JDxl6SGpo5/kvOT/DjJgSTPJvnysPycJI8keWG43jz9dCVJs5Bpj/NPshXYWlVPJTkL2At8HvhT4I2qujXJzcDmqvrKtBOWJE1v6m/+VXW4qp4abr8FHADOA3YA9wxPu4fFfxAkSWvA1N/8f+HNkguAx4GPAS9X1a8seey/quo9u36S7AJ2AZx55pm/e8kll8xsPpLUwd69e1+vqi2TvGbTrAZP8gHgu8BNVfVmklN6XVXtBnYDzM3N1fz8/KymJEktJPnPSV8zk6N9kryPxfB/u6q+Nyx+bfh7wNG/CxyZxViSpOnN4mifAHcBB6rqG0seegjYOdzeCTw47ViSpNmYxW6fy4AvAs8k2Tcs+ypwK/CdJF8CXgaumcFYkqQZmDr+VfVvwPF28H962veXJM2ev/CVpIaMvyQ1ZPwlqSHjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqaGZxD/JniRHkuxfsuyWJK8k2TdcrprFWJKk6c3qm//dwBXLLL+9qrYPl4dnNJYkaUoziX9VPQ68MYv3kiSNb+x9/jcmeXrYLbR55LEkSadozPjfAVwEbAcOA7ct96Qku5LMJ5lfWFgYcTqSpKNGi39VvVZV71TVu8C3gEuP87zdVTVXVXNbtmwZazqSpCVGi3+SrUvuXg3sP95zJUkra9Ms3iTJfcDlwAeTHAK+DlyeZDtQwEHg+lmMJUma3kziX1XXLbP4rlm8tyRp9vyFryQ1ZPwlqSHjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqSHjL0kNGX9Jamgm8U+yJ8mRJPuXLDsnySNJXhiuN89iLEnS9Gb1zf9u4Ipjlt0MPFZVFwOPDfclSWvAplm8SVU9nuSCYxbvAC4fbt8D/AvwlVmMp7Xn3idf5sF9r6z2NHSadmw/jz/5vY+s9jS0gsbc539uVR0GGK4/tNyTkuxKMp9kfmFhYcTpaEwP7nuF5w6/udrT0Gl47vCb/sPd0Ey++U+jqnYDuwHm5uZqlaejKWzbejb/cP0nVnsamtAX7nxitaegVTDmN//XkmwFGK6PjDiWJGkCY8b/IWDncHsn8OCIY0mSJjCrQz3vA54APprkUJIvAbcCn0nyAvCZ4b4kaQ2Y1dE+1x3noU/P4v0lSbPlL3wlqSHjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqSHjL0kNGX9JamjT2AMkOQi8BbwDvF1Vc2OPKUk6sdHjP/hUVb2+QmNJkk7C3T6S1NBKxL+AHyXZm2TXsQ8m2ZVkPsn8wsLCCkxHkrQS8b+sqj4OXAnckOSTSx+sqt1VNVdVc1u2bFmB6UiSRo9/Vb06XB8BHgAuHXtMSdKJjRr/JGcmOevobeCzwP4xx5QkndzYR/ucCzyQ5OhY91bVD0YeU5J0EqPGv6peAn57zDEkSZPzUE9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqSHjL0kNGX9Jasj4S1JDxl+SGho9/kmuSPJ8kheT3Dz2eJKkkxs1/knOAL4JXAlsA65Lsm3MMSVJJ7dp5Pe/FHixql4CSHI/sAN4buRxtcKe/OkbAHzhzidWeSaa1HOH32Tb1rNXexpaYWPv9jkP+NmS+4eGZf8nya4k80nmFxYWRp6OpGNt23o2O7afd/InakMZ+5t/lllWv3CnajewG2Bubq6Web7WgYO3/tFqT0HSBMb+5n8IOH/J/Q8Dr448piTpJMaO/0+Ai5NcmOT9wLXAQyOPKUk6iVF3+1TV20luBH4InAHsqapnxxxTknRyY+/zp6oeBh4eexxJ0qnzF76S1JDxl6SGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqSHjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNTRa/JPckuSVJPuGy1VjjSVJmsymkd//9qr665HHkCRNyN0+ktTQ2PG/McnTSfYk2TzyWJKkUzRV/JM8mmT/MpcdwB3ARcB24DBw23HeY1eS+STzCwsL00xHknSKUlXjD5JcAHy/qj52oufNzc3V/Pz86PORpI0kyd6qmpvkNWMe7bN1yd2rgf1jjSVJmsyYR/v8VZLtQAEHgetHHEuSNIHR4l9VXxzrvSVJ0/FQT0lqyPhLUkPGX5IaMv6S1JDxl6SGjL8kNWT8Jakh4y9JDRl/SWrI+EtSQ8Zfkhoy/pLUkPGXpIaMvyQ1ZPwlqSHjL0kNGX9Jasj4S1JDxl+SGjL+ktSQ8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkPGX5IaMv6S1NBU8U9yTZJnk7ybZO6Yx/4iyYtJnk/yuemmKUmapU1Tvn4/8MfAnUsXJtkGXAv8FvBrwKNJfqOq3plyPEnSDEz1zb+qDlTV88s8tAO4v6p+XlU/BV4ELp1mLEnS7Ez7zf94zgP+fcn9Q8Oy90iyC9g13P15kv0jzWkt+CDw+mpPYkSu3/q1kdcNNv76fXTSF5w0/kkeBX51mYe+VlUPHu9lyyyr5Z5YVbuB3cNY81U1t9zzNgLXb33byOu3kdcNeqzfpK85afyr6g9PYy6HgPOX3P8w8OppvI8kaQRjHer5EHBtkl9OciFwMfAfI40lSZrQtId6Xp3kEPAJ4J+S/BCgqp4FvgM8B/wAuOEUj/TZPc181gHXb33byOu3kdcNXL/3SNWyu+IlSRuYv/CVpIaMvyQ1tCbi3+k0EUluSfJKkn3D5arVntO0klwxbJ8Xk9y82vOZtSQHkzwzbK+JD6lba5LsSXJk6W9qkpyT5JEkLwzXm1dzjtM4zvptiM9dkvOT/DjJgaGZXx6WT7z91kT8+f/TRDy+dOExp4m4Avi7JGes/PRm7vaq2j5cHl7tyUxj2B7fBK4EtgHXDdtto/nUsL02wrHid7P4eVrqZuCxqroYeGy4v17dzXvXDzbG5+5t4M+r6jeB3wduGD5vE2+/NRF/TxOxrl0KvFhVL1XV/wD3s7jdtEZV1ePAG8cs3gHcM9y+B/j8ik5qho6zfhtCVR2uqqeG228BB1g8e8LE229NxP8EzgN+tuT+cU8Tsc7cmOTp4b+n6/a/14ONuo2WKuBHSfYOpyPZiM6tqsOwGBjgQ6s8nzFspM8dSS4Afgd4ktPYfisW/ySPJtm/zOVE3xJP+TQRa8lJ1vUO4CJgO3AYuG1VJzu9dbmNJnRZVX2cxV1bNyT55GpPSBPbUJ+7JB8AvgvcVFVvns57jHVit/fodJqIU13XJN8Cvj/ydMa2LrfRJKrq1eH6SJIHWNzV9fiJX7XuvJZka1UdTrIVOLLaE5qlqnrt6O31/rlL8j4Ww//tqvresHji7bfWd/tsuNNEDBvmqKtZ/GP3evYT4OIkFyZ5P4t/oH9olec0M0nOTHLW0dvAZ1n/22w5DwE7h9s7geOdtHFd2iifuyQB7gIOVNU3ljw08fZbE7/wTXI18LfAFuC/gX1V9bnhsa8Bf8biX7lvqqp/XrWJzkCSv2fxv54FHASuP7qvbr0aDpv7G+AMYE9V/eUqT2lmkvw68MBwdxNw73pfvyT3AZezeJrj14CvA//I4ilZPgK8DFxTVevyj6bHWb/L2QCfuyR/APwr8Azw7rD4qyzu959o+62J+EuSVtZa3+0jSRqB8Zekhoy/JDVk/CWpIeMvSQ0Zf0lqyPhLUkP/Cy0GHi63zq5KAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pointpats.window import to_ccf\n",
    "p0 = np.asarray(to_ccf(parts[0]))\n",
    "plt.plot(p0[:,0], p0[:,1])\n",
    "plt.xlim(-10,20)\n",
    "t=plt.ylim(-10,20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can print all the rings composing our window: two exterior rings, and one hole:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAD7dJREFUeJzt3W+MZXV9x/H3R/6oARKgrnSLJCClKm3qQCcbGhqz9V+BJ0ijCWxieEAypoFEE5tINWkx6QM1Io8M7RiIpHFFqxKIEhUICzFp0Vkdll03lNXSurDZHWOJ+IQW/PbBPdtMYWZnZu85c+/u7/1Kbu45v3vO+X1zcudzf3PuueekqpAkteV1ky5AkrT5DH9JapDhL0kNMvwlqUGGvyQ1yPCXpAaNHf5J3pDkh0meTLIvyae79ouSPJHkmSRfS3L6+OVKkvrQx8j/JeDdVfVOYAa4KskVwGeBO6rqEuC/gJt66EuS1IOxw79GftPNntY9Cng38I2u/R7gA+P2JUnqx6l9bCTJKcBu4PeBLwI/A16oqpe7RQ4C56+y7hwwB3DGGWf8ydvf/vY+SpKkZuzevfuXVbVlI+v0Ev5V9Qowk+Rs4D7gHSsttsq688A8wOzsbC0sLPRRkiQ1I8l/bHSdXs/2qaoXgF3AFcDZSY5+uLwFeL7PviRJx6+Ps322dCN+krwReC+wH3gU+GC32I3A/eP2JUnqRx+HfbYC93TH/V8HfL2qvp3kp8C9Sf4e+AlwVw99SZJ6MHb4V9Ue4LIV2n8ObBt3+5Kk/vkLX0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDerjBu4XJHk0yf4k+5J8tGu/LclzSRa7xzXjlytJ6kMfN3B/Gfh4Vf04yVnA7iQPda/dUVWf76EPSVKP+riB+yHgUDf9YpL9wPnjbleSNJxej/knuRC4DHiia7olyZ4kdyc5p8++JEnHr7fwT3Im8E3gY1X1a+BO4GJghtF/Brevst5ckoUkC0tLS32VI0k6hl7CP8lpjIL/K1X1LYCqOlxVr1TVb4EvAdtWWreq5qtqtqpmt2zZ0kc5kqQ19HG2T4C7gP1V9YVl7VuXLXYdsHfcviRJ/ejjbJ8rgQ8DTyVZ7No+CdyQZAYo4FngIz30JUnqQR9n+/wAyAovPTjutiVJw/AXvpLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDxg7/JBckeTTJ/iT7kny0az83yUNJnumezxm/XElSH/oY+b8MfLyq3gFcAdyc5FLgVuCRqroEeKSblyRNgVPH3UBVHQIOddMvJtkPnA9cC2zvFrsH2AV8Ytz+pGk0Pw87d066ihPHjh0wNzfpKtrW6zH/JBcClwFPAOd1HwxHPyDevMo6c0kWkiwsLS31WY60aXbuhMXFSVdxYlhc9INyGow98j8qyZnAN4GPVdWvk6xrvaqaB+YBZmdnq696pM02MwO7dk26ium3ffukKxD0NPJPchqj4P9KVX2raz6cZGv3+lbgSB99SZLG18fZPgHuAvZX1ReWvfQAcGM3fSNw/7h9SZL60cdhnyuBDwNPJTl61POTwGeArye5CfhP4EM99CVJ6kEfZ/v8AFjtAP97xt2+JKl//sJXkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDegn/JHcnOZJk77K225I8l2Sxe1zTR1+SpPH1NfL/MnDVCu13VNVM93iwp74kSWPqJfyr6nHgV31sS5I0vKGP+d+SZE93WOiclRZIMpdkIcnC0tLSwOVIkmDY8L8TuBiYAQ4Bt6+0UFXNV9VsVc1u2bJlwHIkSUcNFv5VdbiqXqmq3wJfArYN1ZckaWMGC/8kW5fNXgfsXW1ZSdLmOrWPjST5KrAdeFOSg8DfAduTzAAFPAt8pI++JEnj6yX8q+qGFZrv6mPbkqT++QtfSWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGtTL9fyl1j322Oh5+/aJlnFCWFyEmZlJVyFH/pI21cwM7Ngx6SrkyF/qQdWkK5A2ppeRf5K7kxxJsndZ27lJHkryTPd8Th99SZLG19dhny8DV72q7Vbgkaq6BHikm5ckTYG+buD+eJILX9V8LbC9m74H2AV8oo/+tLb5edi5c9JVnFx27IC5uUlXIfVjyC98z6uqQwDd85tXWijJXJKFJAtLS0sDltOWnTtHZ1WoH4uLfpjq5DLxL3yrah6YB5idnfVrsx7NzMCuXZOu4uTgKZw62Qw58j+cZCtA93xkwL4kSRswZPg/ANzYTd8I3D9gX5KkDejrVM+vAv8CvC3JwSQ3AZ8B3pfkGeB93bwkaQr0dbbPDau89J4+ti9J6peXd5CkBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBE7+Hr05g8/Obc1fzHTtgbm74fqSGOPLX8du5ExYXh+1jcXFzPmCkxjjy13hmZmDXruG2v337cNuWGjZ4+Cd5FngReAV4uapmh+5TknRsmzXy//Oq+uUm9SVJWoPH/CWpQZsR/gV8P8nuJK85ZSPJXJKFJAtLS0ubUI4kaTPC/8qquhy4Grg5ybuWv1hV81U1W1WzW7Zs2YRyJEmDh39VPd89HwHuA7YN3ack6dgGDf8kZyQ56+g08H5g75B9SpLWNvTZPucB9yU52tfOqvruwH1KktYwaPhX1c+Bdw7ZhyRp4zzVU5IaZPhLUoMMf0lqkOEvSQ3yqp46fo89Nnoe8sqbi4ujK4dK6pUjf023mZnRzVwk9cqRv45f1aQrkHScHPlLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNGjz8k1yV5OkkB5LcOnR/kqS1DRr+SU4BvghcDVwK3JDk0iH7lCStbeirem4DDnQ3cifJvcC1wE8H7rd5m3Gp/ZZ4WwGdbIY+7HM+8Itl8we7tv+TZC7JQpKFpaWlgcuRjo+3FdDJZuiRf1Zo+38Xga+qeWAeYHZ21gvE98RL7Us6lqFH/geBC5bNvwV4fuA+JUlrGDr8fwRckuSiJKcD1wMPDNynJGkNgx72qaqXk9wCfA84Bbi7qvYN2ackaW2D38O3qh4EHhy6H0nS+vkLX0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDRos/JPcluS5JIvd45qh+pIkbczQ9/C9o6o+P3AfkqQN8rCPJDVo6PC/JcmeJHcnOWfgviRJ6zRW+Cd5OMneFR7XAncCFwMzwCHg9lW2MZdkIcnC0tLSOOVIktYpVTV8J8mFwLer6o+Otdzs7GwtLCwMXo8knUyS7K6q2Y2sM+TZPluXzV4H7B2qL0nSxgx5ts/nkswABTwLfGTAviRJGzBY+FfVh4fatiRpPJ7qKUkNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDVorPBP8qEk+5L8Nsnsq177myQHkjyd5C/GK1OS1Kdxb+C+F/hL4B+XNya5FLge+EPg94CHk/xBVb0yZn+SpB6MNfKvqv1V9fQKL10L3FtVL1XVvwMHgG3j9CVJ6s+4I//VnA/867L5g13baySZA+a62ZeS7B2opj69CfjlpItYB+vs14lQ54lQI1hn39620RXWDP8kDwO/u8JLn6qq+1dbbYW2WmnBqpoH5ru+FqpqdqXlpol19ss6+3Mi1AjW2bckCxtdZ83wr6r3HkctB4ELls2/BXj+OLYjSRrAUKd6PgBcn+T1SS4CLgF+OFBfkqQNGvdUz+uSHAT+FPhOku8BVNU+4OvAT4HvAjev80yf+XHq2UTW2S/r7M+JUCNYZ982XGeqVjwUL0k6ifkLX0lqkOEvSQ2aivA/ES8TkeS2JM8lWewe10y6pqOSXNXtrwNJbp10PatJ8mySp7r9t+FT1YaS5O4kR5b/5iTJuUkeSvJM93zOJGvsalqpzql7Xya5IMmjSfZ3f+cf7dqnap8eo86p2adJ3pDkh0me7Gr8dNd+UZInun35tSSnr7mxqpr4A3gHox8p7AJml7VfCjwJvB64CPgZcMqk6+1quw3460nXsUJdp3T76a3A6d3+u3TSda1S67PAmyZdxwp1vQu4HNi7rO1zwK3d9K3AZ6e0zql7XwJbgcu76bOAf+v+tqdqnx6jzqnZp4x+Q3VmN30a8ARwBaMTbK7v2v8B+Ku1tjUVI//yMhF92gYcqKqfV9V/A/cy2o9ap6p6HPjVq5qvBe7ppu8BPrCpRa1glTqnTlUdqqofd9MvAvsZ/eJ/qvbpMeqcGjXym272tO5RwLuBb3Tt69qXUxH+x3A+8Itl86teJmJCbkmyp/v3e+KHATrTvs+WK+D7SXZ3l/mYZudV1SEYhQTw5gnXcyzT+L4EIMmFwGWMRqxTu09fVSdM0T5NckqSReAI8BCj//RfqKqXu0XW9Te/aeGf5OEke1d4HGtUuu7LRAxhjZrvBC4GZoBDwO2bVdcaJrrPNujKqrocuBq4Ocm7Jl3QSWBa35ckORP4JvCxqvr1pOtZzQp1TtU+rapXqmqG0ZUTtjE6bP6axdbazlAXdnuNOgEvE7HempN8Cfj2wOWs1wlzaY2qer57PpLkPkZv5McnW9WqDifZWlWHkmxlNOqaOlV1+Oj0NL0vk5zGKFC/UlXf6pqnbp+uVOe07tOqeiHJLkbH/M9Ocmo3+l/X3/y0H/aZ2stEdG/Wo65jdG+DafAj4JLu2//TGd1X4YEJ1/QaSc5IctbRaeD9TM8+XMkDwI3d9I3Aahc1nKhpfF8mCXAXsL+qvrDspanap6vVOU37NMmWJGd3028E3svou4lHgQ92i61vX0762+vu2+nrGI1YXwIOA99b9tqnGB3Tehq4etK1Lqvrn4CngD2M3sRbJ13TstquYXSmws8YXX114jWtUONbGZ2J9CSwb5rqBL7K6N/7/+nelzcBvwM8AjzTPZ87pXVO3fsS+DNGhyH2AIvd45pp26fHqHNq9inwx8BPulr2An/btb+V0cD4APDPwOvX2paXd5CkBk37YR9J0gAMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktSg/wXw/Hl3XOzMNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for part in parts:\n",
    "    part = np.asarray(to_ccf(part))\n",
    "    plt.plot(part[:,0], part[:,1], 'b')\n",
    "for hole in holes:\n",
    "    hole = np.asarray(to_ccf(hole))\n",
    "    plt.plot(hole[:,0], hole[:,1], 'r')\n",
    "plt.xlim(-10,30)\n",
    "t = plt.ylim(-10,30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The red hole is associated with the first exterior ring.\n",
    "\n",
    "With this visual representation, consider the problem of testing whether or not this multi-part window contains one or more points in a sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAETpJREFUeJzt3X9sZWd95/H3pyamVUAilIHOhpEm0JQfu9oa1oqwWCEvgTbknzRVkSASyR9IRisigdRKTVupTaVKoVUB7R8Vu2YTNawYKC0gIopKgxUDldyA05owYcQmsNPOkFFi1CLoP3hjvv3jnqncxJ4Zzz3H946f90u6uvc+99zzfHV0/fHjx889J1WFJKktPzXpAiRJB8/wl6QGGf6S1CDDX5IaZPhLUoMMf0lq0Njhn+Snk3wtyTeSPJbk97v265I8nOTxJH+WZHb8ciVJfehj5P9j4M1V9YvAHHBTkjcAfwh8uKquB/4ZeHcPfUmSejB2+NfIv3RPr+puBbwZ+Iuu/X7gV8btS5LUj+f1sZMkM8AjwM8DfwJ8B/hBVT3TbXIWuHaP9y4BSwBXX331f3n1q1/dR0mS1IxHHnnk+1V1ZD/v6SX8q2obmEvyIuCzwGt222yP9y4DywDz8/O1vr7eR0mS1Iwk/7Df9/S62qeqfgCsAm8AXpTk/C+XlwNP9tmXJOny9bHa50g34ifJzwBvAU4BDwG/1m12B/C5cfuSJPWjj2mfo8D93bz/TwGfqqrPJ/kW8MkkfwD8PXBvD31JknowdvhX1aPA63Zp/y5ww7j7lyT1z2/4SlKDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqUB8XcD+W5KEkp5I8luR9XfvdSb6XZKO73Tx+uZKkPvRxAfdngF+vqr9L8kLgkSQPdq99uKr+uIc+JEk96uMC7ueAc93jHyU5BVw77n4lScPpdc4/yXHgdcDDXdOdSR5Ncl+Sa/rsS5J0+XoL/yQvAD4NvL+qfgh8BHglMMfoL4MP7vG+pSTrSdY3Nzf7KkeSdAG9hH+SqxgF/8er6jMAVfVUVW1X1U+AjwI37Pbeqlquqvmqmj9y5Egf5UiSLqKP1T4B7gVOVdWHdrQf3bHZrcDJcfuSJPWjj9U+bwTeBXwzyUbX9tvAO5PMAQWcBt7TQ1+SpB70sdrnb4Ds8tIXxt23JGkYfsNXkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JatDY4Z/kWJKHkpxK8liS93XtL07yYJLHu/trxi9XktSHPkb+zwC/XlWvAd4AvDfJa4G7gJWquh5Y6Z5L/87amTXu+eo9rJ1Zm3QpUlOeN+4OquoccK57/KMkp4BrgVuAxW6z+4FV4DfH7U+Hx9qZNW782I1sbW8xOzPLyu0rLBxbmHRZl2V5GU6cmHQVV47bboOlpUlX0bZe5/yTHAdeBzwMvKz7xXD+F8RL93jPUpL1JOubm5t9lqMpt3p6la3tLbZrm63tLVZPr066pMt24gRsbEy6iivDxoa/KKfB2CP/85K8APg08P6q+mGSS3pfVS0DywDz8/PVVz2afovHF5mdmf23kf/i8cVJlzSWuTlYXZ10FdNvcXHSFQh6Cv8kVzEK/o9X1We65qeSHK2qc0mOAk/30ZcOj4VjC6zcvsLq6VUWjy9esVM+0pVo7PDPaIh/L3Cqqj6046UHgDuAD3T3nxu3Lx0+C8cWDH1pAvoY+b8ReBfwzSTnZz1/m1HofyrJu4F/BN7eQ1+SpB70sdrnb4C9JvhvHHf/kqT++Q1fSWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kN6iX8k9yX5OkkJ3e03Z3ke0k2utvNffQlSRpfXyP/PwVu2qX9w1U1192+0FNfkqQx9RL+VfUV4J/62JckaXhDz/nfmeTRblromt02SLKUZD3J+ubm5sDlSJJg2PD/CPBKYA44B3xwt42qarmq5qtq/siRIwOWI0k6b7Dwr6qnqmq7qn4CfBS4Yai+JEn7M1j4Jzm64+mtwMm9tpUkHazn9bGTJJ8AFoGXJDkL/B6wmGQOKOA08J4++pIkja+X8K+qd+7SfG8f+5Yk9c9v+EpSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLB2DtzBr3fPUe1s6sTboUCejpG76S9rZ2Zo0bP3YjW9tbzM7MsnL7CgvHFiZdlhrnyF8a2OrpVba2t9iubba2t1g9vTrpkiTDXxra4vFFZmdmmckMszOzLB5fnHRJktM+0tAWji2wcvsKq6dXWTy+6JSPpoLhLx2AhWMLhr6mitM+ktQgR/5SD7785dH94uJEy7gibGzA3Nykq5Ajf0kHam4Obrtt0lXIkb/Ug6pJVyDtTy8j/yT3JXk6yckdbS9O8mCSx7v7a/roS5I0vr6mff4UuOlZbXcBK1V1PbDSPZckTYG+LuD+lSTHn9V8C7DYPb4fWAV+s4/+dHHLy3DixKSrOFxuuw2WliZdhdSPIf/h+7KqOgfQ3b90t42SLCVZT7K+ubk5YDltOXFitKpC/djY8JepDpeJ/8O3qpaBZYD5+Xn/bdajuTlYXZ10FYeDSzh12Aw58n8qyVGA7v7pAfuSJO3DkOH/AHBH9/gO4HMD9iVJ2oe+lnp+AlgDXpXkbJJ3Ax8A3prkceCt3XNJ0hToa7XPO/d46cY+9i9J6pend5CkBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBE7+Gr65gy8sHc1Xz226DpaXh+5Ea4shfl+/ECdjYGLaPjY2D+QUjNcaRv8YzNwerq8Ptf3FxuH1LDRs8/JOcBn4EbAPPVNX80H1Kki7soEb+/62qvn9AfUmSLsI5f0lq0EGEfwF/neSRJM9ZspFkKcl6kvXNzc0DKEeSdBDh/8aqej3wNuC9Sd6088WqWq6q+aqaP3LkyAGUI0kaPPyr6snu/mngs8ANQ/cpSbqwQcM/ydVJXnj+MfBLwMkh+5QkXdzQq31eBnw2yfm+TlTVXw3cp6bE2pk1Vk+vsnh8kYVjC5MuR9IOg4Z/VX0X+MUh+9B0Wjuzxo0fu5Gt7S1mZ2ZZuX3FXwDSFHGppwaxenqVre0ttmubre0tVk+vTrokSTsY/hrE4vFFZmdmmckMszOzLB5fnHRJknbw3D4axMKxBVZuX3HOX5pShr8Gs3BswdCXppThr8v35S+P7oc88+bGxujMoZJ65Zy/ptvc3OhiLpJ65chfl69q0hVIukyO/CWpQYa/JDXI8Nee1s6scc9X72HtzNqkS5HUM+f8tStPzyAdbo78tStPzyAdboa/duXpGaTDzWkf7crTM0iHm+GvPXl6BrXssF+PwvCXpGdpYcGDc/6S9CwtLHgw/CXpWVpY8DD4tE+Sm4D/AcwA/7uqPjB0n5I0jhYWPAwa/klmgD8B3gqcBb6e5IGq+taQ/UrSuA77goehR/43AE90F3InySeBWwDDf2AHcar9lnhZAR02Q8/5Xwuc2fH8bNf2b5IsJVlPsr65uTlwOdLl8bICOmyGHvlnl7Z/dxL4qloGlgHm5+c9QXxPPNW+pAsZeuR/Fji24/nLgScH7lOSdBFDh//XgeuTXJdkFngH8MDAfUqSLmLQaZ+qeibJncAXGS31vK+qHhuyT0nSxQ2+zr+qvgB8Yeh+JEmXzm/4SlKDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lq0GDhn+TuJN9LstHdbh6qL0nS/gx9Dd8PV9UfD9yHJGmfnPaRpAYNHf53Jnk0yX1Jrhm4L0nSJRor/JN8KcnJXW63AB8BXgnMAeeAD+6xj6Uk60nWNzc3xylHknSJUlXDd5IcBz5fVf/pQtvNz8/X+vr64PVI0mGS5JGqmt/Pe4Zc7XN0x9NbgZND9SVJ2p8hV/v8UZI5oIDTwHsG7EuStA+DhX9VvWuofUuSxuNST0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBY4V/krcneSzJT5LMP+u130ryRJJvJ/nl8cqUJPVp3Au4nwR+FfhfOxuTvBZ4B/Afgf8AfCnJL1TV9pj9SZJ6MNbIv6pOVdW3d3npFuCTVfXjqvp/wBPADeP0JUnqz7gj/71cC/ztjudnu7bnSLIELHVPf5zk5EA19eklwPcnXcQlsM5+XQl1Xgk1gnX27VX7fcNFwz/Jl4Cf2+Wl36mqz+31tl3aarcNq2oZWO76Wq+q+d22mybW2S/r7M+VUCNYZ9+SrO/3PRcN/6p6y2XUchY4tuP5y4EnL2M/kqQBDLXU8wHgHUmen+Q64HrgawP1JUnap3GXet6a5CywAPxlki8CVNVjwKeAbwF/Bbz3Elf6LI9TzwGyzn5ZZ3+uhBrBOvu27zpTtetUvCTpEPMbvpLUIMNfkho0FeF/JZ4mIsndSb6XZKO73Tzpms5LclN3vJ5Ictek69lLktNJvtkdv30vVRtKkvuSPL3zOydJXpzkwSSPd/fXTLLGrqbd6py6z2WSY0keSnKq+zl/X9c+Vcf0AnVOzTFN8tNJvpbkG12Nv9+1X5fk4e5Y/lmS2YvurKomfgNew+hLCqvA/I721wLfAJ4PXAd8B5iZdL1dbXcDvzHpOnapa6Y7Tq8AZrvj99pJ17VHraeBl0y6jl3qehPweuDkjrY/Au7qHt8F/OGU1jl1n0vgKPD67vELgf/b/WxP1TG9QJ1Tc0wZfYfqBd3jq4CHgTcwWmDzjq79fwL//WL7moqRf3maiD7dADxRVd+tqi3gk4yOoy5RVX0F+KdnNd8C3N89vh/4lQMtahd71Dl1qupcVf1d9/hHwClG3/ifqmN6gTqnRo38S/f0qu5WwJuBv+jaL+lYTkX4X8C1wJkdz/c8TcSE3Jnk0e7P74lPA3Sm/ZjtVMBfJ3mkO83HNHtZVZ2DUUgAL51wPRcyjZ9LAJIcB17HaMQ6tcf0WXXCFB3TJDNJNoCngQcZ/aX/g6p6ptvkkn7mDyz8k3wpycldbhcalV7yaSKGcJGaPwK8EpgDzgEfPKi6LmKix2yf3lhVrwfeBrw3yZsmXdAhMK2fS5K8APg08P6q+uGk69nLLnVO1TGtqu2qmmN05oQbGE2bP2ezi+1nqBO7PUddgaeJuNSak3wU+PzA5VyqK+bUGlX1ZHf/dJLPMvogf2WyVe3pqSRHq+pckqOMRl1Tp6qeOv94mj6XSa5iFKgfr6rPdM1Td0x3q3Naj2lV/SDJKqM5/xcleV43+r+kn/lpn/aZ2tNEdB/W825ldG2DafB14Pruv/+zjK6r8MCEa3qOJFcneeH5x8AvMT3HcDcPAHd0j+8A9jqp4URN4+cySYB7gVNV9aEdL03VMd2rzmk6pkmOJHlR9/hngLcw+t/EQ8CvdZtd2rGc9H+vu/9O38poxPpj4Cngizte+x1Gc1rfBt426Vp31PV/gG8CjzL6EB+ddE07aruZ0UqF7zA6++rEa9qlxlcwWon0DeCxaaoT+ASjP+//f/e5fDfws8AK8Hh3/+IprXPqPpfAf2U0DfEosNHdbp62Y3qBOqfmmAL/Gfj7rpaTwO927a9gNDB+Avhz4PkX25end5CkBk37tI8kaQCGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWrQvwLm5kR1teEIvAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pnts = [(12,12), (4,4), (2,2), (25,1), (5,20)]\n",
    "for pnt in pnts:\n",
    "    plt.plot(pnt[0], pnt[1], 'g.')\n",
    "\n",
    "for part in parts:\n",
    "    part = np.asarray(to_ccf(part))\n",
    "    plt.plot(part[:,0], part[:,1], 'b')\n",
    "for hole in holes:\n",
    "    hole = np.asarray(to_ccf(hole))\n",
    "    plt.plot(hole[:,0], hole[:,1], 'r')\n",
    "plt.xlim(-10,30)\n",
    "t = plt.ylim(-10,30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of the five points two are clearly outside of both of the exterior rings. The three remaining points are each contained in one of the bounding boxes for an exterior ring. However, one of these points is also contained in the hole ring, and thus is not contained in the exterior ring associated with that hole.\n",
    "\n",
    "We can create a Window object from the parts and holes to demonstrate how to evaluate these containment checks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pointpats import Window\n",
    "window = Window(parts, holes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[(0.0, 0.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0), (0.0, 0.0)],\n",
       " [(11.0, 11.0), (11.0, 20.0), (20.0, 20.0), (20.0, 11.0), (11.0, 11.0)]]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "window.parts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[(3.0, 3.0), (3.0, 6.0), (6.0, 6.0), (6.0, 3.0), (3.0, 3.0)]]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "window.holes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.0, 0.0, 20.0, 20.0]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "window.bbox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "172.0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "window.area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFkFJREFUeJzt3X+QH3V9x/HXy3CnHicBL4GccEyIJVHC1MPeMF6sztX4AxkYpKJVZoTM2DmmlY46MCPitMVOK+qITv9wrLEwYCsqFa1AHRVvOKhDilzMAYnpAR4pFziSu2AuxGu5JL77x3cvXpLv/Ui+u9/d++7zMfOd7373u7ufN8s3r+/efvf7/joiBAAol1fkXQAAoP4IfwAoIcIfAEqI8AeAEiL8AaCECH8AKKGaw9/2q2z/wvZjtrfZ/mwy/xzbj9h+yvZ3bTfXXi4AIA1pHPm/LOkdEfEmSZ2SLrL9FklfkPSViDhX0m8kfTSFsQAAKag5/KNif/KwKbmFpHdI+l4y/w5J76t1LABAOk5KYyO2l0jaLOkPJH1V0q8l7Y2Ig8kiOyWdOcu6vZJ6Jenkk0/+oze84Q1plAQApbF58+bxiFh+POukEv4RcUhSp+1TJf1A0hurLTbLuhslbZSkrq6uGBgYSKMkACgN2/9zvOukerVPROyV1C/pLZJOtT395nKWpOfTHAsAcOLSuNpneXLEL9uvlvROSdslPSDpimSxqyX9sNaxAADpSOO0T7ukO5Lz/q+QdFdE3Gf7V5K+Y/vvJW2RdGsKYwEAUlBz+EfE45IuqDJ/WNKFtW4fAJA+vuELACVE+ANACRH+AFBChD8AlBDhDwAlRPgDQAkR/gBQQoQ/AJQQ4Q8AJUT4A0AJEf4AUEKEPwCUEOEPACVE+ANACRH+AFBChD8AlBDhDwAlRPgDQAml8QPuHbYfsL3d9jbbH0/m32T7OduDye3i2ssFAKQhjR9wPyjpuoj4pe3XSNps+/7kua9ExJdSGAMAkKI0fsB9VNJoMv2S7e2Szqx1uwCA7KR6zt/2SkkXSHokmXWt7cdt32b7tDTHAgCcuNTC33arpLslfSIi9kn6mqTXS+pU5S+DW2ZZr9f2gO2BsbGxtMoBAMwhlfC33aRK8H8rIr4vSRGxKyIORcTvJH1D0oXV1o2IjRHRFRFdy5cvT6McAMA80rjax5JulbQ9Ir48Y377jMUul7S11rEAAOlI42qft0r6iKQnbA8m826U9GHbnZJC0g5J16QwFgAgBWlc7fNzSa7y1I9q3TYAIBt8wxcASojwB4ASIvwBoIQIfwAoIcIfAEqI8AeAEiL8AaCECH8AKCHCHwBKiPAHgBIi/AGghAh/ACghwh8ASojwB4ASIvwBoIQIfwAoIcIfAEqI8AeAEiL8AaCEag5/2x22H7C93fY22x9P5r/W9v22n0ruT6u9XABAGtI48j8o6bqIeKOkt0j6mO3zJN0gqS8izpXUlzwGjrBpZJNu/s+btWlkU96lAKVyUq0biIhRSaPJ9Eu2t0s6U9JlknqSxe6Q1C/pU7WOh4V5fuPzevEnL+r8u8+XJA1/elgTmybmXKepremI5Q/sOaA1G9dIkoZ6hzT55OSc67esbjli+aa2Jq26eZUkaev7t+rAngNHLL/v5X167IXH1KIWPapHdfDig3rbV992ePlTuk/R2defLUna0rNl3v/mtkvajlh+xYYVat/QrqnxKW27Ytu86x+9fMd1HVp26TJNDk1q6JqhOdcdHZX+valDTy5bprbJSV365JD6zlmlkaVL1TExofXPDM87/tHL37t6jfa0tGj1+LjW7RyZd/2jl79r7VpNNjWr84VRdb7wwrzrH7387Z0XSJLWjTyr1Xv2zLv+zOXP2rdPd62tvJbWDw+rY9+Rr71XX3KG/uxfXzfvNpGdVM/5214p6QJJj0g6I3ljmH6DOH2WdXptD9geGBsbS7OcUju476B+u+23eZcxp73/t1ehUEQoFHp24tm8Szphu3dJzzyTdxWLw/KJ/frf+3blXUbpOSLS2ZDdKulBSf8QEd+3vTciTp3x/G8iYs7z/l1dXTEwMJBKPSi+TSObtP6b6zV1aErNS5rVd1Wfuju68y7rhPT0VO77+/OsYnG4/dTKX3Eb9l6QcyWNw/bmiOg6nnVqPu2TDNwk6W5J34qI7yezd9luj4hR2+2SdqcxFhpHd0e3+q7qU/+OfvWs7Fm0wQ8sRjWHv21LulXS9oj48oyn7pF0taTPJ/c/rHUsLNz0OfIL+ot9dNXd0U3oAzlI48j/rZI+IukJ24PJvBtVCf27bH9U0rOSPpDCWAAWucmmprxLgNK52ufnkjzL0+tr3T6AxjJ9FdBf5lxH2fENXwAooVQ+8AWAhVo/PP2dh1W51lF2hD+Aumo5eGD+hZA5wh9AXd27uvIt8OtyrqPsOOcPACXEkT+Aurr0yek+SWtyraPsCP8G1XZJW94lAFW1Tc7dIBD1Qfg3qOnulgBQDef8AaCECP8GtaVny4J64AMoJ077NKgVG1bkXQKAAiP8G1T7hva8SwBQYJz2aVBT41OaGp/KuwwABcWRf4Oa/s3aovfzB5APwh9AXe1pacm7BIjwB1Bn9PYpBs75A0AJceQPoK7o7VMMqRz5277N9m7bW2fMu8n2c7YHk9vFaYwFYHGbPKlJkyfxO755S+u0z+2SLqoy/ysR0ZncfpTSWAAWsb5Vq9S3il/xylsq4R8RD0l6MY1tAQCyl/UHvtfafjw5LXRatQVs99oesD0wNjaWcTkA8vbBbVv1wW1b518Qmcoy/L8m6fWSOiWNSrql2kIRsTEiuiKia/ny5RmWUy4rNqygvw8KqeXAAbUc4Hd885bZ1T4RsWt62vY3JN2X1Vg4Fr19AMwlsyN/2zPT53JJ/J1XR/T2ATCXVI78bX9bUo+kZbZ3SvpbST22OyWFpB2SrkljLCwMvX0AzCWV8I+ID1eZfWsa28aJ6biuI+8SABQY3/BtUMsuXZZ3CQAKjN4+DWpyaFKTQ5N5lwGgoDjyb1BD11T6p3DOH0A1hD+Auho5ZWneJUCEP4A6m+7r89c511F2nPMHgBLiyB9AXf2+r8/5udZRdoQ/gLraecopeZcAEf5AXWwa2aT+Hf3qWdmj7o7uvMvJ1cMdZ+ddAkT4A5nbNLJJ67+5XlOHptS8pFl9V/WV/g0A+SP8gYz17+jX1KEpHYpDmjo0pf4d/aUO/w2DW5IpvoOSJ8K/QdHbpzh6VvaoeUnz4SP/npU9eZcEEP6Nit4+xdHd0a2+q/o4549CIfwb1HRfn5Y1LTlXAqnyBkDoo0gI/wZFbx8AcyH8G9Sqz63Ku4RSefDByn1PT65lLApX7JdaW/OuAoR/g1q6juZZKKbWVun0M/KuAoR/g5p4eEISbwL1EpF3BYvHlp68K4CUUmM327fZ3m1764x5r7V9v+2nkvvT0hgLCzN847CGbxzOuwwABZVWV8/bJV101LwbJPVFxLmS+pLHAEqu7ZI2tV3SlncZpZfWD7g/ZHvlUbMvk9STTN8hqV/Sp9IYD/MbHZV275I+2ZN3JY3jyiul3t68q1j8zr6e3j5FkGU//zMiYlSSkvvTqy1ku9f2gO2BsbGxDMspl927pP37866icQwOSnfemXcVQHpy/8A3IjZK2ihJXV1dfGyWotZWqb8/7yoaA5dwpmdLT6W3D99ByVeW4b/LdntEjNpul7Q7w7EALBIrNqzIuwQo2/C/R9LVkj6f3P8ww7EALBLtG9rzLgFK71LPb0vaJGmN7Z22P6pK6L/L9lOS3pU8BlByU+NTmhqfyruM0kvrap8Pz/LU+jS2D6BxbLtimyTO+ect9w98kY2+cyq9fTbkWwaAgiL8G9TIUto6AJgd4d+gOiYmkineBAAci/BvUOufme7rw3lVAMci/BvUvavXSOKcP4DqCP8GtaeFn28EMDvCv0GtHh9PpvghdwDHIvwb1LqdI8kU4Q/gWIQ/gLqit08xEP4A6orePsWQZT9/ADgGvX2KgSN/AHVFb59iIPwB1FXHdR15lwAR/gDqbNmlXIFWBJzzB1BXk0OTmhyazLuM0uPIHydu48b6/Kr5lVdKvb3Zj4O6GLpmSBLn/PPGkX+Dunf1msP9fTJz553S4GC2YwwO1ucNBigZjvwbVN16+3R2Sv392W2/pye7bQMllnn4294h6SVJhyQdjIiurMcEvX0AzK1eR/5/EhHj8y+GtNDbB8BcOO3ToO5au1YS/fwBVFePD3xD0k9tb7Z9zCUbtnttD9geGBsbq0M55TDZ1KzJpua8ywBQUPU48n9rRDxv+3RJ99v+74h4aPrJiNgoaaMkdXV1RR3qKYXOF0aTKZpoAThW5uEfEc8n97tt/0DShZIemnst1KrzhReSKcIfwLEyDX/bJ0t6RUS8lEy/W9LfZTkmgGKjt08xZH3kf4akH9ieHuvOiPhxxmOiIDaNbFL/jn71rOxRd0d33uWgIOjtUwyZhn9EDEt6U5ZjoJg2jWzS+m+u19ShKTUvaVbfVX28AUCSDvf1aVlTpy8ioiraOyAT/Tv6NXVoSofikKYOTal/R3/eJaEghq4ZOtzfB/nhOn9komdlj5qXNB8+8u9Z2ZN3SSiIVZ9blXcJEOGPjHR3dKvvqj7O+eMYS9ctzbsEiPBHhro7ugl9HGPi4QlJvAnkjfDHiXvwwcp9lp03BwcrnUPRMIZvHJZEP/+8Ef4NqmF6+3R2Vn7MBUCqCP8GVZe+PkE3DmCxIvwbFL19AMyF8G9Q9PYBMBfCv0Hd3ln5MG1DDdugPQPQuAh/VEV7BqCx0d6hQa0beVbrRp494fVpzwA0No78G9TqPXuSqbNPaH3aMwCNjfBHVbRnQFbo7VMMhD9mRXsGZGGxtHVo9AseCH8AdbUYevuU4YIHPvAFUFfDNw4f7u9TVGW44IEjfwB1tebra/IuYV5luOAh8/C3fZGkf5S0RNI/R8Tnsx4TQHEthp9vLMMFD5mGv+0lkr4q6V2Sdkp61PY9EfGrLMcFUFzj945LKv4PuTf6BQ9ZH/lfKOnp5IfcZfs7ki6TRPhnbO+E1KkJ3X7qljmXe/isDj25bJnaJid16ZND6jtnlUaWLlXHxITWPzP/edmjl7939RrtaWnR6vFxrds5Mu/6Ry9/19q1mmxqVucLozP6E83u6OWn21qsG3l2xncdZjdz+bP27dNda8+XJK0fHlbHvonDy12xX2ptlbb0/H7dprYmnX93ZfnhTw/rwJ4DWrOxckpjqHdIk09Ozjl2y+qWI5ZvamvSqpsrl0Fuff9WHdhzYM71l3YvPWL5U7pP0dnXV77XsaVn7v/vktR2SdsRy6/YsELtG9o1NT6lbVdsm3f9o5fvuK5Dyy5dpsmhyTl/o3f/4H61drYWPvwbXdYf+J4paWYC7EzmHWa71/aA7YGxsbGMyymPT+oCDaq4V1MsNq2t0uln5F1FY2jtbNUZV7Iz8+bIsCe77Q9Iek9E/Hny+COSLoyIv6q2fFdXVwwMDGRWDwA0ItubI6LreNbJ+sh/p6SOGY/PkvR8xmMCAOaRdfg/Kulc2+fYbpb0IUn3ZDwmAGAemX7gGxEHbV8r6SeqXOp5W0TM/0kSACBTmV/nHxE/kvSjrMcBACwc7R0AoIQIfwAoIcIfAEqI8AeAEiL8AaCECH8AKCHCHwBKiPAHgBIi/AGghAh/ACghwh8ASojwB4ASIvwBoIQIfwAoIcIfAEqI8AeAEiL8AaCECH8AKKHMwt/2Tbafsz2Y3C7OaiwAwPHJ+jd8vxIRX8p4DADAceK0DwCUUNbhf63tx23fZvu0jMcCACxQTeFv+2e2t1a5XSbpa5JeL6lT0qikW2bZRq/tAdsDY2NjtZQDAFggR0T2g9grJd0XEefPtVxXV1cMDAxkXg8ANBLbmyOi63jWyfJqn/YZDy+XtDWrsQAAxyfLq32+aLtTUkjaIemaDMcCAByHzMI/Ij6S1bYBALXhUk8AKCHCHwBKiPAHgBIi/AGghAh/ACghwh8ASojwB4ASIvwBoIQIfwAoIcIfAEqI8AeAEiL8AaCECH8AKCHCHwBKiPAHgBIi/AGghAh/ACghwh8ASojwB4ASqin8bX/A9jbbv7PdddRzn7b9tO0h2++prUwAQJpq/QH3rZL+VNLXZ860fZ6kD0laK+l1kn5me3VEHKpxPABACmo68o+I7RExVOWpyyR9JyJejohnJD0t6cJaxgIApKfWI//ZnCnpv2Y83pnMO4btXkm9ycOXbW/NqKY0LZM0nncRC0Cd6VoMdS6GGiXqTNua411h3vC3/TNJK6o89ZmI+OFsq1WZF9UWjIiNkjYmYw1ERFe15YqEOtNFnelZDDVK1Jk22wPHu8684R8R7zyBWnZK6pjx+CxJz5/AdgAAGcjqUs97JH3I9ittnyPpXEm/yGgsAMBxqvVSz8tt75TULek/bP9EkiJim6S7JP1K0o8lfWyBV/psrKWeOqLOdFFnehZDjRJ1pu2463RE1VPxAIAGxjd8AaCECH8AKKFChP9ibBNh+ybbz9keTG4X513TNNsXJfvrads35F3PbGzvsP1Esv+O+1K1rNi+zfbumd85sf1a2/fbfiq5Py3PGpOaqtVZuNel7Q7bD9jenvw7/3gyv1D7dI46C7NPbb/K9i9sP5bU+Nlk/jm2H0n25XdtN8+7sYjI/Sbpjap8SaFfUteM+edJekzSKyWdI+nXkpbkXW9S202Srs+7jip1LUn20ypJzcn+Oy/vumapdYekZXnXUaWut0t6s6StM+Z9UdINyfQNkr5Q0DoL97qU1C7pzcn0ayQ9mfzbLtQ+naPOwuxTVb5D1ZpMN0l6RNJbVLnA5kPJ/H+S9BfzbasQR/5Bm4g0XSjp6YgYjogpSd9RZT9igSLiIUkvHjX7Mkl3JNN3SHpfXYuqYpY6CyciRiPil8n0S5K2q/KN/0Lt0znqLIyo2J88bEpuIekdkr6XzF/QvixE+M/hTEkjMx7P2iYiJ9fafjz58zv30wCJou+zmULST21vTtp8FNkZETEqVUJC0uk51zOXIr4uJUm2V0q6QJUj1sLu06PqlAq0T20vsT0oabek+1X5S39vRBxMFlnQv/m6hb/tn9neWuU211HpgttEZGGemr8m6fWSOiWNSrqlXnXNI9d9dpzeGhFvlvReSR+z/fa8C2oARX1dynarpLslfSIi9uVdz2yq1FmofRoRhyKiU5XOCReqctr8mMXm205Wjd2OEYuwTcRCa7b9DUn3ZVzOQi2a1hoR8Xxyv9v2D1R5IT+Ub1Wz2mW7PSJGbberctRVOBGxa3q6SK9L202qBOq3IuL7yezC7dNqdRZ1n0bEXtv9qpzzP9X2ScnR/4L+zRf9tE9h20QkL9Zpl6vy2wZF8Kikc5NP/5tV+V2Fe3Ku6Ri2T7b9mulpSe9WcfZhNfdIujqZvlrSbE0Nc1XE16VtS7pV0vaI+PKMpwq1T2ers0j71PZy26cm06+W9E5VPpt4QNIVyWIL25d5f3qdfDp9uSpHrC9L2iXpJzOe+4wq57SGJL0371pn1PUvkp6Q9LgqL+L2vGuaUdvFqlyp8GtVuq/mXlOVGlepciXSY5K2FalOSd9W5c/7A8nr8qOS2iT1SXoquX9tQess3OtS0h+rchricUmDye3iou3TOeoszD6V9IeStiS1bJX0N8n8VaocGD8t6d8kvXK+bdHeAQBKqOinfQAAGSD8AaCECH8AKCHCHwBKiPAHgBIi/AGghAh/ACih/wf9OCIsBL9jbQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pnts = [(12,12), (4,4), (2,2), (25,1), (5,20)]\n",
    "for pnt in pnts:\n",
    "    plt.plot(pnt[0], pnt[1], 'g.') #plot the five points in green\n",
    "\n",
    "for part in parts:\n",
    "    part = np.asarray(to_ccf(part))\n",
    "    plt.plot(part[:,0], part[:,1], 'b') #plot \"parts\" in blue\n",
    "for hole in holes:\n",
    "    hole = np.asarray(to_ccf(hole))\n",
    "    plt.plot(hole[:,0], hole[:,1], 'r') #plot \"hole\" in red \n",
    "    \n",
    "from pointpats.window import poly_from_bbox\n",
    "poly = np.asarray(poly_from_bbox(window.bbox).vertices)\n",
    "plt.plot(poly[:,0], poly[:,1], 'm-.') #plot the minimum bounding box in magenta\n",
    "\n",
    "plt.xlim(-10,30)\n",
    "t = plt.ylim(-10,30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have extended the figure to include the bounding box for the multi-part window (in cyan). Now we can call the `filter_contained` method of the window on the point sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([12, 12]), array([2, 2])]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pin = window.filter_contained(pnts)\n",
    "pin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This was a lot of code just to illustrate that the methods of a window can be used to identify topological relationships between points and the window's constituent parts. Let's turn to a less contrived example to see this in action.\n",
    "\n",
    "Here we will make use of PySAL's [shapely extension](https://pysal.readthedocs.org/en/latest/users/tutorials/shapely.html) to create a multi-part window from the county shapefile for Virgina."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from libpysal.cg import shapely_ext\n",
    "import numpy as np\n",
    "from pointpats.window import poly_from_bbox, as_window, Window\n",
    "import libpysal as ps\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "va = ps.io.open(ps.examples.get_path(\"vautm17n.shp\")) #open \"vautm17n\" polygon shapefile\n",
    "polys = [shp for shp in va]\n",
    "vapnts = ps.io.open(ps.examples.get_path(\"vautm17n_points.shp\")) #open \"vautm17n_points\" point shapefile\n",
    "points = [shp for shp in vapnts]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "136\n"
     ]
    }
   ],
   "source": [
    "print(len(polys))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The county shapefile `vautm17n.shp` has 136 shapes of the `polygon` type. Some of these are composed of multiple-rings and holes to reflect the [interesting history](https://en.wikipedia.org/wiki/List_of_counties_in_Virginia) of political boundaries in that State.\n",
    "Fortunately, with our window class we can handle these. We will come back to this shortly.\n",
    "\n",
    "First we are going to build up a realistic window for our point pattern based on a *cascaded union* made possible via [Shapely](https://pypi.python.org/pypi/Shapely) through the [PySAL shapely extension](https://pysal.readthedocs.org/en/latest/users/tutorials/shapely.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "cu = shapely_ext.cascaded_union(polys)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This creates a PySAL Polygon:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "libpysal.cg.shapes.Polygon"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(cu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can construct a Window from this polygon instance using the helper function `as_window`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = as_window(cu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[]]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.holes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(w.parts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The window has three parts consisting of the union of mainland counties and two \"island\" parts associated with Accomack and Northampton counties and has no holes.\n",
    "\n",
    "Since this a window, we can access its properties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[260694.99205079858, 4044845.4484747574, 1005496.0048517315, 4370839.043748417]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.bbox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(689097.7340935213, 4155195.0497352206)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.centroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.contains_point(w.centroid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the centroid for our new window is contained by the window.  Such a result is not guaranteed as the geometry of the window could be complex such that the centroid falls outside of the window.\n",
    "\n",
    "Let's continue on with a more interesting query. Since we know the window centroid is contained in the Window, we can find which individual county contains the centroid. \n",
    "\n",
    "Our strategy is a simple one to illustrate the useful nature of the Window. We will create a sequence of Windows, one for each county and use them to carry out a containment test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create a window for each of the individual counties in the state\n",
    "windows = [as_window(county) for county in polys]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(67, <pointpats.window.Window at 0x11e37c9e8>)]"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#check each county for containment of the window's centroid\n",
    "cent_poly = [ (i, county) for i,county in enumerate(windows) if county.contains_point(w.centroid)]\n",
    "cent_poly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "i, cent_poly = cent_poly[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[674997.5183093206, 4119217.2472937624, 713300.2226730094, 4159075.43995212]"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cent_poly.bbox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we did here was create a window for each of the individual counties in the state. With these in hand we checked each one for containment of the window's centroid. The result is we see the window (count) with index 67 is the only one that contains the centroid point.\n",
    "\n",
    "The point of this exercise is not to use an inefficient brute force exhaustive search to find this county. There are more efficient spatial indices in PySAL that we could use for such a query. Rather, we wanted to explicitly check each window to ensure that only one contained the centroid.\n",
    "\n",
    "As we will see in elsewhere in this series of notebooks, this type of decomposition can support highly flexible types of spatial analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Windows and point pattern intensity revisited\n",
    "\n",
    "Returning to the central use of Windows, we saw in the [introductory notebook](pointpattern.ipynb) that the area of the Window is used to form the estimate of intensity for the point pattern:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Point Pattern\n",
      "200 points\n",
      "Bounding rectangle [(273959.664381352,4049220.903414295), (972595.9895779632,4359604.85977962)]\n",
      "Area of window: 216845506675.0557\n",
      "Intensity estimate for window: 9.223156295311261e-10\n",
      "               x             y\n",
      "0  865322.486181  4.150317e+06\n",
      "1  774479.213103  4.258993e+06\n",
      "2  308048.692232  4.054700e+06\n",
      "3  670711.529980  4.258864e+06\n",
      "4  666254.475614  4.256514e+06\n"
     ]
    }
   ],
   "source": [
    "f = ps.examples.get_path('vautm17n_points.shp') #open \"vautm17n_points\" point shapefile\n",
    "fo = ps.io.open(f)\n",
    "pnts = np.asarray([pnt for pnt in fo])\n",
    "fo.close()\n",
    "pp_va = PointPattern(pnts)\n",
    "pp_va.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the default is to form the minimum bounding rectangle and use that as the window for the point pattern and, in turn, to implment the intesity estimation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can override the default by passing a window object in to the constructor for the point pattner. Here we use our window that was formed from the county cascading union above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Point Pattern\n",
      "200 points\n",
      "Bounding rectangle [(273959.664381352,4049220.903414295), (972595.9895779632,4359604.85977962)]\n",
      "Area of window: 103195696155.68987\n",
      "Intensity estimate for window: 1.9380653210407425e-09\n",
      "               x             y\n",
      "0  865322.486181  4.150317e+06\n",
      "1  774479.213103  4.258993e+06\n",
      "2  308048.692232  4.054700e+06\n",
      "3  670711.529980  4.258864e+06\n",
      "4  666254.475614  4.256514e+06\n"
     ]
    }
   ],
   "source": [
    "pp_va_union = PointPattern(pnts, window=w)\n",
    "pp_va_union.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the window is redefined. Thus, window related attributes **Area of window** and **Intensity estimate for window** are changed. However, the **Bounding rectangle** remains unchanged since it is not relavant to the definition of window.\n",
    "\n",
    "Close examination of the summary report from reveals that while the bounding rectangles for the two point pattern instances are identical (as they should be), the area of the windows are substantially different:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.1013037825521717"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_va.window.area / pp_va_union.window.area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "as are the intensity estimates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.47589501732368955"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_va.lambda_window / pp_va_union.lambda_window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}