Intersect points with a regular grid without a spatial index.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Intersect points with a regular grid without a spatial index\n",
    "- Source: https://github.com/carlosg-m/goose-lab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.6/dist-packages/geopandas/_compat.py:110: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.1-CAPI-1.14.2). Conversions between both will be slow.\n",
      "  shapely_geos_version, geos_capi_version_string\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pyproj\n",
    "import matplotlib.pyplot as plt\n",
    "import geopandas as gpd\n",
    "import shapely "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Intersect points with a regular grid using Numpy\n",
    "- Intersecing points with a regular grid on a cartesian plane can be done efficiently using only NumPy and vectorized operations, without relying on spatial indexes or brute force approaches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_grid(xmin, ymin, xmax, ymax, xstep, ystep):\n",
    "    \"\"\"\n",
    "    Create bounding coordinates for a regular 2d grid with numpy.\n",
    "    These coordinates can then be used with shapely, pygeos and geopandas.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Generate two evenly spaced ranges\n",
    "    x = np.linspace(xmin, xmax, int(xstep+1))\n",
    "    y = np.linspace(ymin, ymax, int(ystep+1))\n",
    "    \n",
    "    # Get center points for each range\n",
    "    x = np.mean((x[1:], x[:-1]), axis=0)\n",
    "    y = np.mean((y[1:], y[:-1]), axis=0)\n",
    "    \n",
    "    # Cartesian product of both ranges \n",
    "    x, y = np.meshgrid(x, y)\n",
    "    x, y = x.ravel(), y.ravel()\n",
    "\n",
    "    # Calculate delta x and delta y\n",
    "    xd = (xmax - xmin) / (xstep * 2)\n",
    "    yd = (ymax - ymin) / (ystep * 2)\n",
    "    \n",
    "    # Repeat center points and add respective delta\n",
    "    grid = np.column_stack((x,y))\n",
    "    grid = np.tile(grid, 2) + np.array([[-xd, -yd, xd, yd]])\n",
    "    \n",
    "    # Generate index\n",
    "    grid_index = np.arange(0, xstep * ystep)\n",
    "    \n",
    "    # Output coordinates: xmin, ymin, xmax, ymax\n",
    "    return grid, grid_index\n",
    "\n",
    "def intersect_points_grid(xy, xmin, ymin, xmax, ymax, xstep, ystep, include_xmax=False, include_ymax=False):\n",
    "    \"\"\"\n",
    "    Intersect 2d Points with a regular 2d Grid using NumPy.\n",
    "    Returns the index of the Grid rectangle that each Point belongs to.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Calculate delta x and delta y\n",
    "    xd = (xmax - xmin) / xstep\n",
    "    yd = (ymax - ymin) / ystep\n",
    "    \n",
    "    # Center points at zero\n",
    "    xy = xy - np.array([[xmin, ymin]])\n",
    "    \n",
    "    # Get index of x and y coordinates\n",
    "    xi = np.floor(xy[:,0] / xd) \n",
    "    yi = np.floor(xy[:,1] / yd)\n",
    "    \n",
    "    # Cast x and y index to 64 bit integer\n",
    "    xi = xi.astype(np.int64)\n",
    "    yi = yi.astype(np.int64)\n",
    "    \n",
    "    # Opened interval for x\n",
    "    if include_xmax == True:\n",
    "        xi = xi - (xy[:,0] % xd == 0)\n",
    "    \n",
    "    # Opened interval for y\n",
    "    if include_ymax == True:\n",
    "        yi = yi - (xy[:,1] % yd == 0)\n",
    "        \n",
    "    # Mask out-of-bound indexes on the x and y axis\n",
    "    xi = np.ma.array(xi, mask=((xi < 0) | (xi >= xstep)))\n",
    "    yi = np.ma.array(yi, mask=((yi < 0) | (yi >= ystep)))\n",
    "        \n",
    "    # Calculate index\n",
    "    index = xi + (yi * xstep)\n",
    "    \n",
    "    # Unmask NumPy array\n",
    "    index = index.filled(np.iinfo(np.int64).min)\n",
    "    \n",
    "    return index\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Create example grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid, grid_index = create_grid(xmin=10, ymin=50, xmax=20, ymax=60, xstep=10, ystep=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Create example points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_points = np.array([[10, 50],\n",
    "                        [10, 60],\n",
    "                        [20, 50],\n",
    "                        [20, 60],\n",
    "                        [11, 51],\n",
    "                        [14, 54],\n",
    "                        [11.5, 57.7],\n",
    "                        [17.5, 54.5]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Plot grid and Points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fdd51361780>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAFlCAYAAADGe3ILAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAz7klEQVR4nO3de3RU5cH+/e8dUmAFoWppOAUyKEJiQjJACvIWUeQNKuKBQzk4XU8kaJQ+y1JUlP5oVfwZUagiFotEUFiQAp4gvBYB5aA8KmLQIKgRpAQSQAgiB4lAQu73jxweoAk5MHtnm1yftbIys2dmX/cdwjV77jnEWGsRERHvCanrAYiISMVU0CIiHqWCFhHxKBW0iIhHqaBFRDxKBS0i4lGhboa1bNnS+nw+NyNFRDxt8+bNh6y1v67oMlcL2ufzkZmZ6WakiIinGWN2V3aZljhERDxKBS0i4lEqaBERj1JBi4h4lApaRMSjVNAiIh6lghYR8SgVtIiIR6mgRUQ8qloFbYy51BjzhjEm2xjztTGmtzHmcmPMu8aYHaXfL3NigOnp6fh8PkJCQvD5fKSnpzsRIyJSY073U3WPoGcAK621UUA88DUwEVhjrb0KWFN6PqjS09NJSUlh9+7dWGvZvXs3KSkpKmkRqXNu9JOp6m8SGmN+CWQBV9izrmyM+Qa43lq73xjTBlhvre1yoX0lJCTYmnwWh8/nY/fu3VwSfyPNrr6+fHvTpk3o1euaau+nNg4eOwnGEN68iaM5DTm3Ic21rnIb0lzdzv3kk42cPHkKgBNfrefHLasAiIyMJCcnp9r7McZsttYmVHRZdY6gOwL5wKvGmM+NMXOMMc2AVtba/aXX+Q5oVUl4ijEm0xiTmZ+fX+1BA+zZsweAZldfT+PwjuXby34oTtpzuIA9hwscz2nIuQ1prnWV25Dm6nZuWQ81Du94zgFkWW8FQ3U+zS4U6A7cb639xBgzg/OWM6y11hhT4aG4tTYNSIOSI+iaDK5Dhw7s3l3yQU+nD+7iwKI/AyX3UEv+GfQVlXN0fbzk3nDJvb0dzWnIuQ1prnWV25Dm6naub8oodu/eTatRU87Z3qFDh6BlVOcIOg/Is9Z+Unr+DUoK+0Dp0gal3w8GbVSlUlNTCQsLO2dbWFgYqampwY4SEakRN/qpyoK21n4H5BpjytaX+wNfAcuBpNJtSUBG0EZVKhAIkJaWRtOmJetJkZGRpKWlEQgEgh0lIlIjbvRTdT+w/34g3RjTGPg3MJqScn/NGDMG2A0MD9qozhIIBFj+48cAji9riIjUhNP9VK2CttZmARU9y9g/qKMREZFyeiehiIhHqaBFRDxKBS0i4lEqaBERj1JBi4h4lApaRMSjVNAiIh6lghYR8SgVtIiIR6mgRUQ8SgUtIuJRKmgREY9SQYuIeJQKWkTEo1TQIiIepYIWEfEoFbSIiEepoEVEPEoFLSLiUSpoERGPqpcFfebMGbp168agQYMA2LVrF7169aJTp06MGDGC06dPO545c+ZMOnXqhDGGQ4cOBT2vstxAIECXLl2IjY0lOTmZwsJCV3LHjBlDfHw8cXFxDBs2jB9//NHxzDJ//OMfueSSS4KeV1nuXXfdRceOHfH7/fj9frKyslzJtdYyadIkOnfuTHR0NC+88ILjmddee235PNu2bcsdd9wR9MyKctesWUP37t3x+/306dOHb7/91pXctWvX0r17d2JjY0lKSqKoqMiR3OqqlwU9Y8YMoqOjy88/8sgjjB8/nm+//ZbLLruMuXPnOp7529/+lvfee4/IyMigZ10oNxAIkJ2dzdatW/npp5+YM2eOK7nTp09ny5YtfPHFF3To0IGZM2c6ngmQmZnJDz/8EPSsqnKnTZtGVlYWWVlZ+P1+V3LnzZtHbm4u2dnZfP3114wcOdLxzA0bNpTPs3fv3gwZMiTomRXljh07lvT0dLKysrjzzjt58sknHc8tLi4mKSmJxYsXs23bNiIjI5k/f74judVV7wo6Ly+Pf/3rX9x9991AyVHH2rVrGTZsGABJSUksW7bM0UyAbt264fP5gppTndyBAwdijMEYQ8+ePcnLy3Mlt0WLFkDJz/unn37CGON45pkzZ5gwYQJTp04NalZVuW6oKHfWrFk8+uijhISU/LcNDw93PLPMsWPHWLt2rSNH0BXlGmM4duwYAEePHqVt27aO537//fc0btyYzp07A5CYmMibb74Z9NyaqHcF/ac//YmpU6eW/xJ///33XHrppYSGhgIQERHB3r17Hc10y4VyCwsLWbBgATfddJNruaNHj6Z169ZkZ2dz//33O545c+ZMbrvtNtq0aRPUrKpyASZNmkRcXBzjx4/n1KlTruTu3LmTJUuWkJCQwM0338yOHTsczyyzbNky+vfvX35H7HTunDlzGDhwIBERESxYsICJEyc6ntuyZUuKiorIzMwE4I033iA3NzfouTVRrwr67bffJjw8nB49etTrzOrk/uEPf6Bv375ce+21ruW++uqr7Nu3j+joaJYsWeJo5r59+3j99deDfkdQVS7AlClTyM7O5tNPP+Xw4cM888wzruSeOnWKpk2bkpmZyT333ENycrLjmWUWLVrEqFGjgpZXVe706dNZsWIFeXl5jB49mgceeMDxXGMMixcvZvz48fTs2ZPmzZvTqFGjoObWVGidpgfZhx9+yPLly1mxYgUnT57k2LFjjBs3jiNHjlBUVERoaCh5eXm0a9fO0czf//73LFy4MGgZNc2dPHky+fn5zJ4929VcgEaNGjFy5EimTp3K6NGjHcuMiYmhSZMmdOrUCYCCggI6deoU1CeTqpprkyZNGD16NH/729+Clnmh3IiIiPI14MGDBwft53uhzIULF3Lo0CE2bdrE0qVLg5Z3odxbbrmF7OxsevXqBcCIESOC/kjwQvPdsGEDAKtXr2b79u1Bza0xa61rXz169LC1Mfylj+zwlz6q0W3WrVtnb7nlFmuttcOGDbOLFi2y1lp777332hdffLHK28c+ttLGPray1pllIiMjbX5+frX3cbG5L7/8su3du7ctKCio0T4uJre4uNju2LHDWmttcXGxffDBB+2DDz7oaOb5mjVrVu19XGzuvn37rLUlcx03bpx95JFHXMl95JFH7Ny5c8u3JyQkOJ5prbWzZs2y//Vf/1WjfVxMbmFhof3Vr35lv/nmG2uttXPmzLFDhgxxPNdaaw8cOGCttfbkyZP2hhtusGvWrKny9rXpp7MBmbaSzqxXSxyVeeaZZ3juuefo1KkT33//PWPGjHE884UXXiAiIoK8vDzi4uJce5Lpvvvu48CBA/Tu3Ru/388TTzzheKa1lqSkJLp27UrXrl3Zv38/jz76qOO5dSUQCJTP9dChQ/zlL39xJXfixIm8+eabdO3alT//+c+OvULnfIsXL3ZkeaMyoaGhvPzyywwdOpT4+HgWLFjAtGnTXMmeNm0a0dHRxMXFceutt3LDDTe4klsZU1Lg7khISLBlC/A1MWL2xwAsubd3sIdUqa6PrwJg6+M3upbZ0HIb0lzrKrchzbWuci+2n4wxm621CRVd1iCOoEVEfo5U0CIiHqWCFhHxKBW0iIhHqaBFRDxKBS0i4lEqaBERj1JBi4h4lApaRMSjVNAiIh6lghYR8SgVtIiIR6mgRUQ8SgUtIuJRKmgREY9SQYuIeJQKWkTEo1TQIiIepYIWEfEoFbSIiEepoEVEPCq0OlcyxuQAx4EzQJG1NsEY4wdeApoCRcAfrLWbHBqniEiDU62CLtXPWnvorPNTgcnW2neMMQNLz18fzMGJiDRkF7PEYYEWpad/Cey7+OGIiEgZY62t+krG7AJ+oKSUZ1tr04wx0cAqwFBS9P+PtXb3hfaTkJBgMzMzazzIfn9bz57DBYQ1blTj29bW8ZNFADRvWpMHGcr1emZDy21Ic62r3OMni+jYshnrHrq+Vrc3xmy21iZUdFl1Z9HHWrvXGBMOvGuMyQaGAeOttW8aY4YDc4H/t4LwFCAFoEOHDrWaANaWfImIeEwjA+HNmziy72oVtLV2b+n3g8aYpUBPIAkYV3qV14E5ldw2DUiDkiPo2gwyvEVTwls0Zcm9vWtz81rp+vgqALY+fqNrmQ0ttyHNta5yG9Jc6yp3xOyPHdt3lWvQxphmxpjmZaeBAcA2Stacryu92g3ADqcGKSLSEFXnCLoVsNQYU3b9f1prVxpjfgRmGGNCgZOULmOIiEhwVFnQ1tp/A/EVbP8foIcTgxIREb2TUETEs1TQQZaeno7P5yMkJASfz0d6enpdD0lEfqbcfZFiPZeenk5KSgoFBQUA7N69m5SUkqX5QCBQl0MTkZ8hHUEH0aRJk8rLuUxBQQGTJk2qoxGJyM+ZCjqI9uzZU6PtIiIXooIOosreKVnrd1CKSIOmgg6i1NRUwsLCztkWFhZGampqHY1IRH7OVNBBFAgESEtLIzIyEmMMkZGRpKWl6QlCEakVvYojyAKBgApZRIJCR9AiIh6lghYR8SgVtIiIR6mgRUQ8SgUtIuJRKmgREY9SQYuIeJQKWkTEo1TQIiIepYIWEfEoFbSIiEfVq4I+efIkPXv2JD4+npiYGB577DEAdu3aRa9evejUqRMjRozg9OnTruTOnDmTTp06YYzh0KFDrmQGAgG6dOlCbGwsycnJFBYWupI7ZswY4uPjiYuLY9iwYfz444+u5Jb54x//yCWXXOJK5l133UXHjh3x+/34/X6ysrJcybXWMmnSJDp37kx0dDQvvPCCK7nXXntt+Vzbtm3LHXfc4XjmmjVr6N69O36/nz59+vDtt98GLfNCuWvXrqV79+7ExsaSlJREUVFRUHNrql59WFKTJk1Yu3Ytl1xyCYWFhfTp04ebb76Z5557jvHjxzNy5Ejuu+8+5s6dy9ixYx3P/e1vf8ugQYO4/vrrg5ZVVWYgEGDhwoUA3HnnncyZM8eVuU6fPp0WLVoA8MADDzBz5kwmTpzoeO4111xDZmYmP/zwQ9CyqsoEmDZtGsOGDQt65oVyv/76a3Jzc8nOziYkJISDBw+6krthw4by6wwdOpTbb7/d8cyxY8eSkZFBdHQ0//jHP3jyySeZN2+eo7k33ngjSUlJrFmzhs6dO/Poo48yf/58xowZE7TcmqpXR9DGmPKjqMLCQgoLCzHGsHbt2vL/TElJSSxbtsyV3G7duuHz+YKaVVXmwIEDMcZgjKFnz57k5eW5kltWztZafvrpJ4wxruSeOXOGCRMmMHXq1KDmXSjTaZXlzpo1i0cffZSQkJL/tuHh4a7kljl27Bhr164N6hF0ZZnGGI4dOwbA0aNHadu2bdAyK8tt1KgRjRs3pnPnzgAkJiby5ptvBjW3pupVQQOcOXMGv99PeHg4iYmJXHnllVx66aWEhpY8WIiIiGDv3r2O5/bq1SvoGTXJLCwsZMGCBdx0002u5Y4ePZrWrVuTnZ3N/fff70ruzJkzue2222jTpk3Q8yrLhJK/PxkXF8f48eM5deqUK7k7d+5kyZIlJCQkcPPNN7Njxw5XcsssW7aM/v37l98ZO5k5Z84cBg4cSEREBAsWLAjqo7HKcnv27ElRURGZmZkAvPHGG+Tm5gY9tybqXUE3atSIrKws8vLy2LRpE9nZ2XWSu23btjrN/MMf/kDfvn259tprXct99dVX2bdvH9HR0SxZssTx3A8++IDXX3/dkTuDyjK3bdvGlClTyM7O5tNPP+Xw4cM888wzruSeOnWKpk2bkpmZyT333ENycrIruWUWLVrEqFGjXMmcPn06K1asIC8vj9GjR/PAAw84nvvll1+yePFixo8fT8+ePWnevDmNGjUKem5N1LuCLnPppZfSr18/Pv74Y44cOVK+2J+Xl0e7du0cz125cqVjGVVlTp48mfz8fJ577jlXc6Hkl37kyJGOPjQsy123bh3ffvstnTp1wufzUVBQQKdOnRzNXLlyJW3atMEYQ5MmTRg9ejSbNm1yJPP83IiICIYMGQLA4MGD+eKLL1zJBTh06BCbNm3illtucTzznXfeYcuWLeVH7yNGjOCjjz5yPHflypX07t2bDRs2sGnTJvr27Vu+3FFX6lVB5+fnc+TIEQB++ukn3n33XaKjo+nXrx9vvPEGAPPnzw/qkxyV5UZFRQU1o7qZc+bMYdWqVSxatKh8rdLp3C5dupQ/y26tZfny5UGff0W5PXr04LvvviMnJ4ecnBzCwsKC+mx/ZT/j/fv3AyVzXbZsGbGxsUHLvFDuHXfcwbp16wB4//33g14eF/o9fuONNxg0aBBNmzZ1PDM6OpqjR4+yfft2gPJtTudGRUWVP/F66tQpnnnmGe67776g5tZUvXoVx/79+0lKSuLMmTMUFxczfPhwBg0axNVXX83IkSP5y1/+Qrdu3YL+rGxluS+88AJTp07lu+++Iy4ujoEDBzJnzhxHM0NDQ4mMjKR3794ADBkyhEcffTQomZXl3nLLLVx77bUcO3YMay3x8fHMmjUraJmV5Q4aNCioGdXNvOGGG8jPz8dai9/v56WXXnIlt0+fPgQCAaZPn84ll1wStN+lqnIBFi9e7Mg6cGWZL7/8MkOHDiUkJITLLruMV155xZXcCRMm8Pbbb1NcXMzYsWO54YYbgppbU8Za61pYQkKCLVuAr4kRsz8GYMm9vYM9pEp1fXwVAFsfv9G1zIaW25DmWle5DWmudZV7sf1kjNlsrU2o6LJ6tcQhIlKfqKBFRDxKBS0i4lEqaBERj1JBi4h4lApaRMSjVNAiIh6lghYR8SgVtIiIR6mgRUQ8SgUtIuJRKmgREY9SQYuIeJQKWkTEo1TQIiIepYIWEfEoFbSIiEepoEVEPEoFLSLiUSpoERGPUkGLiHiUClpExKOqVdDGmBxjzFZjTJYxJvOs7fcbY7KNMV8aY6Y6N0wRkYYntAbX7WetPVR2xhjTD7gdiLfWnjLGhAd9dCIiDZix1lZ9JWNygITzCvo1IM1a+151wxISEmxmZmbVVzxPv7+tZ8/hAsIaN6rxbWvr+MkiAJo3rcl9mHK9ntnQchvSXOsq9/jJIjq2bMa6h66v1e2NMZuttQkVXVbdNWgLrDbGbDbGpJRu6wxca4z5xBjzvjHmN5WEpxhjMo0xmfn5+TUfPYC1JV8iIh7TyEB48yaO7Lu6dzN9rLV7S5cx3jXGZJfe9nLgGuA3wGvGmCvseYfk1to0IA1KjqBrM8jwFk0Jb9GUJff2rs3Na6Xr46sA2Pr4ja5lNrTchjTXusptSHOtq9wRsz92bN/VOoK21u4t/X4QWAr0BPKAt2yJTUAx0NKpgYqINDRVFrQxppkxpnnZaWAAsA1YBvQr3d4ZaAwcqmQ3IiJSQ9VZ4mgFLDXGlF3/n9balcaYxsArxphtwGkg6fzlDRERqb0qC9pa+28gvoLtp4HfOzEoERHROwlFRDxLBS0i4lEqaBERj1JBi4h4lApaRMSjVNAiIh6lghYR8SgVtIiIR6mgRUQ8SgUtIuJRKmgREY9SQYuIeJQKWkTEo1TQIiIepYIWEfEoFbSIiEepoEVEPEoFLSLiUSpoERGPqlcFnZubS79+/bj66quJiYlhxowZABw+fJjExESuuuoqEhMT+eGHH1zJff3114mJiSEkJITMzExXMidMmEBUVBRxcXEMHjyYI0eOuJL717/+lbi4OPx+PwMGDGDfvn2u5JZ59tlnMcZw6FDw/rB8ZZmPP/447dq1w+/34/f7WbFiRdAyL5QL8Pe//52oqChiYmJ4+OGHXckdMWJE+Vx9Ph9+v9/xzKysLK655hr8fj8JCQls2rQpaJkXyt2yZQu9e/ema9eu3HrrrRw7diyouTVmrXXtq0ePHrY2hr/0kR3+0kdVXm/fvn128+bN1lprjx07Zq+66ir75Zdf2gkTJtgpU6ZYa62dMmWKffjhh6vcV+xjK23sYyurNb7Kcr/66iubnZ1tr7vuOvvpp59Wa1/Vza0sc9WqVbawsNBaa+3DDz9crbkGI/fo0aPl15kxY4a99957g5Z5oVxrrd2zZ48dMGCA7dChg83Pzw9abmWZjz32mJ02bVq1xh3M3LVr19r+/fvbkydPWmutPXDgQNAyL5R7tgceeMBOnjw5aLmVZSYmJtoVK1ZYa63917/+Za+77rpqzeFicxMSEuz69euttdbOnTvX/uUvf6lyX9Xtp8oAmbaSzqxXR9Bt2rShe/fuADRv3pzo6Gj27t1LRkYGSUlJACQlJbFs2TJXcqOjo+nSpUtQs6rKHDBgAKGhJX+s/ZprriEvL8+V3BYtWpRf58SJExhjXMkFGD9+PFOnTnU100mV5c6aNYuJEyfSpEkTAMLDw13JLWOt5bXXXmPUqFGOZxpjyo9ejx49Stu2bYOWeaHc7du307dvXwASExN58803g5pbU/WqoM+Wk5PD559/Tq9evThw4ABt2rQBoHXr1hw4cMCVXLdUlvnKK69w8803u5Y7adIk2rdvT3p6Ok888YQruRkZGbRr1474+HjH8s7PBJg5cyZxcXEkJycHfcmsstzt27ezYcMGevXqxXXXXcenn37qSm6ZDRs20KpVK6666irHM59//nkmTJhA+/bteeihh5gyZYojmefnxsTEkJGRAZQsUebm5jqWWx31sqB//PFHhg4dyvPPP3/OkR2AMSboR1rVyXVKZZmpqamEhoYSCARcy01NTSU3N5dAIMDMmTMdzw0NDeWpp55y9M7g/MwWLVowduxYdu7cSVZWFm3atOHBBx90JbeoqIjDhw+zceNGpk2bxvDhwyl5hOxsbplFixYF9ej5QpmzZs1i+vTp5ObmMn36dMaMGeNK7iuvvMI//vEPevTowfHjx2ncuLEjudVV7wq6sLCQoUOHEggEGDJkCACtWrVi//79AOzfvz/oDw0ry3VaZZnz5s3j7bffJj093ZE7o6rmGggEHHloeH7uzp072bVrF/Hx8fh8PvLy8ujevTvfffedY5lQ8vvUqFEjQkJCuOeee4L+BFZluREREQwZMgRjDD179iQkJCSoT4pWlgtQVFTEW2+9xYgRI4KaV1nm/Pnzy0//7ne/c+1nHBUVxerVq9m8eTOjRo3iyiuvDHpuTdSrgrbWMmbMGKKjo3nggQfKt992223Mnz8fKPmHv/32213JdVJlmStXrmTq1KksX76csLAw13J37NhRfjojI4OoqCjHc7t27crBgwfJyckhJyeHiIgIPvvsM1q3bu1YJlB+Zw+wdOlSYmNjg5JXVe4dd9zBunXrANi+fTunT5+mZcuWjucCvPfee0RFRRERERG0vAtltm3blvfffx+AtWvXBn1ZpbLcgwcPAlBcXMyTTz7JfffdF9Tcmgqt0/Qg+/DDD1mwYAFdu3YtfynQU089xcSJExk+fDhz584lMjKS1157zZXcU6dOcf/995Ofn88tt9yC3+9n1apVjmb+8Y9/5NSpUyQmJgIlTxS+9NJLQcm8UO7cuXP55ptvCAkJITIyMqiZF8odOHBgUHOqk7lo0SKysrIwxuDz+Zg9e7YrucnJySQnJxMbG0vjxo2ZP39+UB8hXehnvHjxYkeWNyrLfPnllxk3bhxFRUU0bdqUtLQ0V3J37NjBiy++CMCQIUMYPXp0UHNrql4VdJ8+fSpdk1uzZk2d5A4ePNjVTCcLy4u5Z8vJyXElsy7nunDhwjrJnTdvnuuZmzdvdiSzqtxx48Y5lltT9WqJQ0SkPlFBi4h4lApaRMSjVNAiIh6lghYR8SgVtIiIR6mgRUQ8SgUtIuJRKmgR+Q/p6en4fD5CQkLw+Xykp6fX9ZAapHr1TkIRuXjp6emkpKRQUFAAwO7du0lJSQFw7NMRpWI6ghaRc0yaNKm8nMsUFBQwadKkOhpRw6WCFpFz7Nmzp0bbxTkqaBE5R4cOHWq0XZyjghaRc6Smpv7HZ4mHhYWRmppaRyNquFTQInKOQCBAWloakZGRGGOIjIwkLS1NTxDWAb2KQ0T+QyAQUCF7gI6gRUQ8SgUtIuJRKmgREY9SQYuIeJQKWkTEo1TQIiIeVa2CNsbkGGO2GmOyjDGZ5132oDHGGmNaOjNEqY/S09M5fuwYR48c0aeliVSiJkfQ/ay1fmttQtkGY0x7YACgN+lLtZV9WlpxcTHwv5+WppIWOZex1lZ9JWNygARr7aHztr8B/F8go6LLz5eQkGAzMzMvdJUK9fvbevYcLiCscaMa37a2jp8sAqB5U3ffy9MQco8fO0ZxcTEhTS/BWos9dQKAkJAQmrdo4Xx+A/gZ12VmQ8s9frKIji2bse6h62t1e2PM5rMPfM9W3SNoC6w2xmw2xqSU7vR2YK+1dksV4SnGmExjTGZ+fn6NBv6/6bbkS+qFsiPn8w8OyraL/Jw0MhDevIkj+67u3Uwfa+1eY0w48K4xJhv4P5Qsb1yQtTYNSIOSI+jaDDK8RVPCWzRlyb29a3PzWun6+CoAtj5+o2uZDSXX5/Oxe/du2o9bDEDujJEAREZGsjUnx/H8hvAzrsvMhpY7YvbHju27WkfQ1tq9pd8PAkuB64COwJbS5Y8I4DNjTGuHxin1iD4tTaR6qjyCNsY0A0KstcdLTw8AnrDWhp91nRyqsQYtAv/7Z5MmbQ6huLiYyMhIUlNT9eE8IuepzhJHK2CpMabs+v+01q50dFRS7wUCAZ7eUfpw1IVlDZGfoyoL2lr7byC+iuv4gjUgEREpoXcSioh4lApaRMSjVNAiIh6lghYR8SgVtIiIR6mgRUQ8SgUtIuJRKmgREY9SQYuIeJQKWkTEo1TQIiIepYIWEfEoFbSIiEepoEVEPEoFLSLiUSpoERGPUkGLiHiUClpExKNU0CIiHqWCFhHxqHpV0MnJyYSHhxMbG1u+7fDhwyQmJnLVVVeRmJjIDz/84Eru66+/TkxMDCEhIWRmZgY9s7LcCRMmEBUVRVxcHIMHD+bIkSOOZ/71r38lLi4Ov9/PgAED2LdvX1AzK8st8+yzz2KM4dChQ67kPv7447Rr1w6/34/f72fFihWOZwL8/e9/JyoqipiYGB5++OGgZlaWO2LEiPJ5+nw+/H6/K7lZWVlcc801+P1+EhIS2LRpk+OZW7ZsoXfv3nTt2pVbb72VY8eOBTWzNupVQd91112sXLnynG1PP/00/fv3Z8eOHfTv35+nn37aldzY2Fjeeust+vbtG/S8C+UmJiaybds2vvjiCzp37syUKVMcz5wwYQJffPEFWVlZDBo0iCeeeCKomZXlAuTm5rJ69Wo6dOgQ9MwL5Y4fP56srCyysrIYOHCg45nr1q0jIyODLVu28OWXX/LQQw8FNbOy3CVLlpTPc+jQoQwZMsSV3IcffpjHHnuMrKwsnnjiiaDfIVWUeffdd/P000+zdetWBg8ezLRp04KaWRv1qqD79u3L5Zdffs62jIwMkpKSAEhKSmLZsmWu5EZHR9OlS5egZ1WVO2DAAEJDQwG45ppryMvLczyzRYsW5adPnDiBMSaomZXlQklRTp061ZHMC+U6qaLMWbNmMXHiRJo0aQJAeHi4K7llrLW89tprjBo1ypVcY0z5EezRo0dp27at45nbt28vP6BKTEzkzTffDGpmbdSrgq7IgQMHaNOmDQCtW7fmwIEDdTwi97zyyivcfPPNrmRNmjSJ9u3bk56e7sgRdEUyMjJo164d8fHxruSdbebMmcTFxZGcnOzIstn5tm/fzoYNG+jVqxfXXXcdn376qeOZZ9uwYQOtWrXiqquuciXv+eefZ8KECbRv356HHnoo6I8EKxITE0NGRgZQskSZm5vreGZV6n1Bn80Y49iRltekpqYSGhpKIBBwLS83N5dAIMDMmTMdzysoKOCpp55y7c7gbGPHjmXnzp1kZWXRpk0bHnzwQcczi4qKOHz4MBs3bmTatGkMHz4ca63juWUWLVrkyNFzZWbNmsX06dPJzc1l+vTpjBkzxvHMV155hX/84x/06NGD48eP07hxY8czq1LvC7pVq1bs378fgP379zvy0NBr5s2bx9tvv016errrd0iBQMCVh4Y7d+5k165dxMfH4/P5yMvLo3v37nz33XeOZ7dq1YpGjRoREhLCPffcE/QnsCoSERHBkCFDMMbQs2dPQkJCHHlStCJFRUW89dZbjBgxwpU8gPnz55evd//ud79z5WccFRXF6tWr2bx5M6NGjeLKK690PLMq9b6gb7vtNubPnw+U/KPffvvtdTwiZ61cuZKpU6eyfPlywsLCXMncsWNH+emMjAyioqIcz+zatSsHDx4kJyeHnJwcIiIi+Oyzz2jdurXj2WV3+ABLly6t8JUlwXbHHXewbt06oGS54/Tp07Rs2dLxXID33nuPqKgoIiIiXMkDaNu2Le+//z4Aa9eudWVp5eDBgwAUFxfz5JNPct999zmeWZXQuh5AMI0aNYr169dz6NAhIiIimDx5MhMnTmT48OHMnTuXyMhIXnvtNVdyL7/8cu6//37y8/O55ZZb8Pv9rFq1yvHcKVOmcOrUKRITE4GSJwpfeuklRzNXrFjBN998Q0hICJGRkUHNu1CuGw97K8pdv349WVlZGGPw+XzMnj3b8czk5GSSk5OJjY2lcePGzJ8/P+iPjir7GS9evNjR5Y2Kcl9++WXGjRtHUVERTZs2JS0tzfHMH3/8kRdffBGAIUOGMHr06KBm1ka9KuhFixZVuH3NmjV1kjt48GDXc50urbrIrCz3bDk5Oa7l1sXPGGDhwoV1kjtv3rw6yd28ebPrmePGjXMsszbq/RKHiMjPlQpaRMSjVNAiIh6lghYR8SgVtIiIR6mgRUQ8SgUtIuJRKmgREY9SQYuIeJQKWkTEo1TQIiIepYIWEfEoFbSIiEepoEVEPEoFLSLiUSpoERGPUkGLiHiUClpExKNU0CIiHqWCFhHxKBW0iIhHVeuvehtjcoDjwBmgyFqbYIyZBtwKnAZ2AqOttUccGqeISINTkyPoftZav7U2ofT8u0CstTYO2A78OeijExFpwGq9xGGtXW2tLSo9uxGICM6QREQEwFhrq76SMbuAHwALzLbWpp13+f8HLLHWLrzQfhISEmxmZmaNB9nvb+vZc7iAsMaNanzb2jp+suS+p3nTaq0CKfdnktnQchvSXOsq9/jJIjq2bMa6h66v1e2NMZvPWpk4R3Vn0cdau9cYEw68a4zJttZ+ULrzSUARkF5JeAqQAtChQ4caDx4Aa0u+REQ8ppGB8OZNHNl3tQraWru39PtBY8xSoCfwgTHmLmAQ0N9WciheerSdBiVH0LUZZHiLpoS3aMqSe3vX5ua10vXxVQBsffxG1zIbWm5Dmmtd5TakudZV7ojZHzu27yrXoI0xzYwxzctOAwOAbcaYm4CHgdustQWOjVBEpIGqzhF0K2CpMabs+v+01q40xnwLNKFkyQNgo7X2PsdGKiLSwFRZ0NbafwPxFWzv5MiIREQE0DsJRUQ8SwUtIuJRKmgREY9SQYuIeJQKWkTEo1TQIiIepYIWEfEoFbSIiEepoEVEPEoFLSLiUSpoERGPUkGLiHiUClpExKNU0CIiHqWCFhHxKBW0iIhHqaBFRDxKBS0i4lEqaBERj1JBi4h4VL0u6JUrV9KlSxc6derE008/7UpmcnIy4eHhxMbGupIHkJubS79+/bj66quJiYlhxowZruSePHmSnj17Eh8fT0xMDI899pgruWXOnDlDt27dGDRokCt5Pp+Prl274vf7SUhIcCUT4MiRIwwbNoyoqCiio6P5+OOPHc/85ptv8Pv95V8tWrTg+eefdzx3+vTpxMTEEBsby6hRozh58qTjmQAzZswgNjaWmJgYV+ZZXfW2oM+cOcN///d/88477/DVV1+xaNEivvrqq2rdNj09nePHjnH0yBF8Ph/p6enVzr3rrrtYuXJlbYddK6GhoTz77LN89dVXbNy4kRdffLHac70YTZo0Ye3atWzZsoWsrCxWrlzJxo0bHc8tM2PGDKKjo13LA1i3bh1ZWVlkZma6ljlu3DhuuukmsrOz2bJliytz7tKlC1lZWWRlZbF582bCwsIYPHiwo5l79+7lhRdeIDMzk23btnHmzBkWL17saCbAtm3bePnll9m0aRNbtmzh7bff5ttvv3U8tzrqbUFv2rSJTp06ccUVV9C4cWNGjhxJRkZGlbdLT08nJSWF4uJiAHbv3k1KSkq1S7pv375cfvnlFzX2mmrTpg3du3cHoHnz5kRHR7N3717Hc40xXHLJJQAUFhZSWFiIMcbxXIC8vDz+9a9/cffdd7uSV1eOHj3KBx98wJgxYwBo3Lgxl156qatjWLNmDVdeeSWRkZGOZxUVFfHTTz9RVFREQUEBbdu2dTzz66+/plevXoSFhREaGsp1113HW2+95XhuddTbgt67dy/t27cvPx8REVGt0po0aRIFBQXnbCsoKGDSpElBH6MTcnJy+Pzzz+nVq5creWfOnMHv9xMeHk5iYqJruX/605+YOnUqISHu/QobYxgwYAA9evQgLS3Nlcxdu3bx61//mtGjR9OtWzfuvvtuTpw44Up2mcWLFzNq1CjHc9q1a8dDDz1Ehw4daNOmDb/85S8ZMGCA47mxsbFs2LCB77//noKCAlasWEFubq7judVRbwu6tvbs2VOj7V7y448/MnToUJ5//nlatGjhSmajRo3IysoiLy+PTZs2sW3bNscz3377bcLDw+nRo4fjWWf7n//5Hz777DPeeecdXnzxRT744APHM4uKivjss88YO3Ysn3/+Oc2aNXPt+RSA06dPs3z5cn73u985nvXDDz+QkZHBrl272LdvHydOnGDhwoWO50ZHR/PII48wYMAAbrrpJvx+P40aNXI8tzrqbUG3a9funHvBvLw82rVrV+XtOnToUKPtXlFYWMjQoUMJBAIMGTLE9fxLL72Ufv36ubL+/uGHH7J8+XJ8Ph8jR45k7dq1/P73v3c8t+z3Jzw8nMGDB7Np0ybHMyMiIoiIiCh/ZDJs2DA+++wzx3PLvPPOO3Tv3p1WrVo5nvXee+/RsWNHfv3rX/OLX/yCIUOG8NFHHzmeCzBmzBg2b97MBx98wGWXXUbnzp1dya1KvS3o3/zmN+zYsYNdu3Zx+vRpFi9ezG233Vbl7VJTUwkLCztnW1hYGKmpqU4N9aJZaxkzZgzR0dE88MADruXm5+dz5MgRAH766SfeffddoqKiHM+dMmUKeXl55OTksHjxYm644QbHj7ROnDjB8ePHy0+vXr3alVfqtG7dmvbt2/PNN98AJevBV199teO5ZRYtWuTK8gaUHARt3LiRgoICrLWsWbPGtSeBDx48CJQ8Un7rrbe48847XcmtSmhdD8ApoaGhzJw5kxtvvJEzZ86QnJxMTExMlbcLBAIATNocQnFxMZGRkaSmppZvr8qoUaNYv349hw4dIiIigsmTJ5c/weOUDz/8kAULFpS/BAzgqaeeYuDAgY7m7t+/n6SkJM6cOUNxcTHDhw937SVvbjtw4ED5qxiKioq48847uemmm1zJ/vvf/04gEOD06dNcccUVvPrqq67knjhxgnfffZfZs2e7kterVy+GDRtG9+7dCQ0NpVu3bqSkpLiSPXToUL7//nt+8Ytf8OKLL7r+RGxl6m1BAwwcOLBWJRUIBHh6xyoAtubk1Oi2ixYtqnHexerTpw/WWtdz4+Li+Pzzz13PPdv111/P9ddf73jOFVdcwZYtWxzPqYjf73f1ZX1lmjVrxvfff+9q5uTJk5k8ebKrmQAbNmxwPbM66u0Sh4jIz50KWkTEo1TQIiIepYIWEfEoFbSIiEepoEVEPEoFLSLiUSpoERGPUkGLiHiUClpExKNU0CIiHqWCFhHxKBW0iIhHqaBFRDxKBS0i4lEqaBERj1JBi4h4lApaRMSjVNAiIh6lghYR8SgVtIiIR6mgRUQ8qloFbYzJMcZsNcZkGWMyS7ddbox51xizo/T7ZU4MMD09nU8+2cj777+Pz+cjPT3diRgRkRpzup9qcgTdz1rrt9YmlJ6fCKyx1l4FrCk9H1Tp6emkpKRw8uQpAHbv3k1KSopKWkTqnBv9ZKy1VV/JmBwgwVp76Kxt3wDXW2v3G2PaAOuttV0utJ+EhASbmZlZ7cH5fD52795Nq1FTaBzekdMHdwHQtGkTevW6ptr7qY1Pdh2mUYghIdKRBwbKraPMhpbbkObqdu4nn2zk5MlTNO3QlZN7tnJg0Z8BiIyMJCcnp9r7McZsPuvA9xyh1dyHBVYbYyww21qbBrSy1u4vvfw7oFUl4SlACkCHDh2qPWiAPXv2AHDiq/XnbC+7x3JSx1+FEd6iqeM5DTm3Ic21rnIb0lzdzi3robPLGf63t4KhukfQ7ay1e40x4cC7wP3AcmvtpWdd5wdr7QXvtmp7BH2+mt5DiYgEW7D66UJH0NVag7bW7i39fhBYCvQEDpQubVD6/WC1R1RNqamphIWFnbMtLCyM1NTUYEeJiNSIG/1UZUEbY5oZY5qXnQYGANuA5UBS6dWSgIygjapUIBAgLS2NyMhIjDFERkaSlpZGIBAIdpSISI240U9VLnEYY66g5KgZStas/2mtTTXG/Ap4DegA7AaGW2sPX2hfNV3iEBGp7y7qSUJr7b+B+Aq2fw/0v/jhiYhIRfROQhERj1JBi4h4lApaRMSjVNAiIh6lghYR8SgVtIiIR6mgRUQ8SgUtIuJRKmgREY+q1qfZBS3MmHxK3hZeGy2BQ1Veq37QXOuvhjRfzbV6Iq21v67oAlcL+mIYYzIre796faO51l8Nab6a68XTEoeIiEepoEVEPOrnVNBpdT0AF2mu9VdDmq/mepF+NmvQIiINzc/pCFpEpEHxZEEbY14xxhw0xmw7a9vlxph3jTE7Sr+7+/fcHVLJXKcZY7KNMV8YY5YaYy6twyEGTUVzPeuyB40x1hjTsi7GFmyVzdUYc3/pv+2XxpipdTW+YKvk99hvjNlojMkyxmQaY3rW5RiDxRjT3hizzhjzVem/47jS7UHvKE8WNDAPuOm8bROBNdbaq4A1pefrg3n851zfBWKttXHAduDP59/oZ2oe/zlXjDHtKflbl8H7e/V1bx7nzdUY0w+4HYi31sYAf6uDcTllHv/5bzsVmGyt9QOPlp6vD4qAB621VwPXAP9tjLkaBzrKkwVtrf0AOP/vG94OzC89PR+4w80xOaWiuVprV1tri0rPbgQiXB+YAyr5dwWYDjwM1JsnRCqZ61jgaWvtqdLrHHR9YA6pZL4WaFF6+pfAPlcH5RBr7X5r7Welp48DXwPtcKCjPFnQlWhlrd1fevo7oFVdDsZFycA7dT0Ipxhjbgf2Wmu31PVYXNAZuNYY84kx5n1jzG/qekAO+xMwzRiTS8mjhfrySLCcMcYHdAM+wYGO+jkVdDlb8tKTenO0VRljzCRKHk6l1/VYnGCMCQP+DyUPfxuCUOBySh4WTwBeM8aYuh2So8YC46217YHxwNw6Hk9QGWMuAd4E/mStPXb2ZcHqqJ9TQR8wxrQBKP1ebx4eVsQYcxcwCAjY+vtayCuBjsAWY0wOJUs5nxljWtfpqJyTB7xlS2wCiin5DIf6Kgl4q/T060C9eJIQwBjzC0rKOd1aWzbHoHfUz6mgl1PyD07p94w6HIujjDE3UbIme5u1tqCux+MUa+1Wa224tdZnrfVRUmDdrbXf1fHQnLIM6AdgjOkMNKZ+f5jQPuC60tM3ADvqcCxBU/qoZy7wtbX2ubMuCn5HWWs99wUsAvYDhZT8px0D/IqSZ0Z3AO8Bl9f1OB2c67dALpBV+vVSXY/Tqbmed3kO0LKux+ngv2tjYCGwDfgMuKGux+nwfPsAm4EtlKzR9qjrcQZprn0oWb744qz/owOd6Ci9k1BExKN+TkscIiINigpaRMSjVNAiIh6lghYR8SgVtIiIR6mgRUQ8SgUtIuJRKmgREY/6/wHuXq/vUQLpFgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "grid_shapely = list(map(lambda x: shapely.geometry.box(*x), grid))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6,6))\n",
    "\n",
    "gpd.GeoDataFrame(geometry=grid_shapely).exterior.plot(ax=ax)\n",
    "\n",
    "for coos, index in zip(grid, grid_index):\n",
    "    ax.text(coos[[0,2]].mean(), coos[[1,3]].mean(), str(index), horizontalalignment='center', verticalalignment='center')\n",
    "    \n",
    "ax.scatter(test_points[:,0], test_points[:,1], color='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Intersect Points with Grid\n",
    "- The parameters that define the grid must be known.\n",
    "- Points that are out of bounds have as index the smallest possible value that can be represented by a 64 bit integer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([                   0, -9223372036854775808, -9223372036854775808,\n",
       "       -9223372036854775808,                    1,                   24,\n",
       "                         31,                   27])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "intersect_points_grid(xy=test_points, \n",
    "                      xmin=10, ymin=50, xmax=20, ymax=60, xstep=10, ystep=5, \n",
    "                      include_xmax=False, include_ymax=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Next steps\n",
    "- Now it's possible to use NumPy indexing or a database join (Pandas, Dask or Pyspark), to merge 2d Points with their respective rectangle on a predefined regular grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = intersect_points_grid(xy=test_points, \n",
    "                              xmin=10, ymin=50, xmax=20, ymax=60, xstep=10, ystep=5, \n",
    "                              include_xmax=False, include_ymax=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10. , 50. ],\n",
       "       [11. , 51. ],\n",
       "       [14. , 54. ],\n",
       "       [11.5, 57.7],\n",
       "       [17.5, 54.5]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_points[index>=0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10., 50., 11., 52.],\n",
       "       [11., 50., 12., 52.],\n",
       "       [14., 54., 15., 56.],\n",
       "       [11., 56., 12., 58.],\n",
       "       [17., 54., 18., 56.]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid[index[index>=0]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Another example - single square"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid, grid_index = create_grid(xmin=0, ymin=0, xmax=1, ymax=1, xstep=1, ystep=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_points = np.array([[0,0],\n",
    "                        [0,1],\n",
    "                        [1,0],\n",
    "                        [1,1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fdd512a6b70>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAFlCAYAAAAtTMkIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAASqElEQVR4nO3df4zfdX3A8ecLutKcoi72MEvbuyuhGCtuAU/AmCiLbJT+Uf7QmTY16iQ2OjBLNEsgLsxgGnTL3GLsorfEbJPSgv6hl1nsMgVNDC0cQRCOoLX0+kMnJ2MkaEqFvfbH99vuy3HX+7b93H376vf5SC75fj7fd7+f16fXPvvt53vfu8hMJElnt/N6PYAkaX7GWpIKMNaSVICxlqQCjLUkFWCsJamAJb068PLly3NkZKRXh5eks87DDz/868wcnO2+nsV6ZGSEiYmJXh1eks46ETE1131eBpGkAoy1JBVgrCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVICxlqQC5o11RHwtIp6JiMfnuD8i4ksRsS8iHouIK5ofs2X79u2MjIxw3nnnMTIywvbt2xfqUJJ0Sha8T5l50g/g3cAVwONz3L8euBcI4Gpg73yPmZm8/e1vz1Nx55135sDAQAInPgYGBvLOO+88pceRpKY11SdgIudoZmQXP4MxIkaAf8/My2a576vA/Zm5o739FHBNZv7yZI85Ojqap/K9QUZGRpiamuK1f3Qdr1l7zYn9y5ZdwFVXXd3140hS0/bu3cPRoy8C8JvJ+3nh0d0ADA8Pc+DAga4fJyIezszR2e5r4pr1CuBQx/bh9r7ZBtkSERMRMTE9PX1KBzl48CAAr1l7DUsvWn1i//HfIEnqleMdWnrR6lc8mTzerSYs6nfdy8wxYAxaz6xP5dcODQ0xNdX6hlTHnnmaX+24FWj9y3X3Xbc0PKkkdW/kjk1MTU3xpk13vGL/0NBQY8do4pn1EWBVx/bK9r5Gbd26lYGBgVfsGxgYYOvWrU0fSpJOyWL0qYlYjwMfan9VyNXA8/Ndrz4dmzdvZmxsjGXLLgBaz6jHxsbYvHlz04eSpFOyGH2a9zJIROwArgGWR8Rh4G+A3wPIzK8Au2h9Rcg+4LfAnzc23QybN29m/IUHALz0IemsstB9mjfWmblpnvsTuKmxiSRJr+I7GCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVICxlqQCjLUkFWCsJakAYy1JBRhrSSrAWEtSAcZakgow1pJUgLGWpAKMtSQVYKwlqQBjLUkFGGtJKsBYS1IBxlqSCjDWklSAsZakAoy1JBVgrCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVICxlqQCjLUkFWCsJakAYy1JBRhrSSrAWEtSAcZakgow1pJUgLGWpAKMtSQVYKwlqQBjLUkFGGtJKsBYS1IBxlqSCjDWklSAsZakAoy1JBVgrCWpAGMtSQV0FeuIWBcRT0XEvoi4ZZb7hyLivoh4JCIei4j1zY8qSf1r3lhHxPnANuB6YC2wKSLWzlj218A9mXk5sBH4p6YHlaR+1s0z6yuBfZm5PzOPATuBG2asSeB17duvB37R3IiSpCVdrFkBHOrYPgxcNWPNZ4H/iIhPAq8Brm1kOkkS0NwLjJuAf8nMlcB64OsR8arHjogtETERERPT09MNHVqSzn3dxPoIsKpje2V7X6cbgXsAMvMBYBmwfOYDZeZYZo5m5ujg4ODpTSxJfaibWD8ErImI1RGxlNYLiOMz1hwE3gsQEW+hFWufOktSQ+aNdWa+BNwM7AaepPVVH09ExO0RsaG97NPAxyLiUWAH8JHMzIUaWpL6TTcvMJKZu4BdM/bd1nF7EnhXs6NJko7zHYySVICxlqQCjLUkFWCsJakAYy1JBRhrSSrAWEtSAcZakgow1pJUgLGWpAKMtSQVYKwlqQBjLUkFGGtJKsBYS1IBxlqSCjDWklSAsZakAoy1JBVgrCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVICxlqQCjLUkFWCsJakAYy1JBRhrSSrAWEtSAcZakgow1pJUgLGWpAKMtSQVYKwlqQBjLUkFGGtJKsBYS1IBxlqSCjDWklSAsZakAoy1JBVgrCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVICxlqQCjLUkFdBVrCNiXUQ8FRH7IuKWOdZ8ICImI+KJiLir2TElqb8tmW9BRJwPbAP+BDgMPBQR45k52bFmDXAr8K7MfC4iLlqogSWpH3XzzPpKYF9m7s/MY8BO4IYZaz4GbMvM5wAy85lmx5Sk/tZNrFcAhzq2D7f3dboUuDQifhQReyJi3WwPFBFbImIiIiamp6dPb2JJ6kNNvcC4BFgDXANsAv45It4wc1FmjmXmaGaODg4ONnRoSTr3dRPrI8Cqju2V7X2dDgPjmfm7zHwa+CmteEuSGtBNrB8C1kTE6ohYCmwExmes+RatZ9VExHJal0X2NzemJPW3eWOdmS8BNwO7gSeBezLziYi4PSI2tJftBp6NiEngPuCvMvPZhRpakvrNvF+6B5CZu4BdM/bd1nE7gU+1PyRJDfMdjJJUgLGWpAKMtSQVYKwlqQBjLUkFGGtJKsBYS1IBxlqSCjDWklSAsZakAoy1JBVgrCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVICxlqQCjLUkFWCs1de++93v8uY3v5lLLrmEz3/+870eR5qTsVbfevnll7npppu49957mZycZMeOHUxOTvZ6LGlWxlp968EHH+SSSy7h4osvZunSpWzcuJFvf/vbvR5LmpWxVt86cuQIq1atOrG9cuVKjhw50sOJpLkZa0kqwFirb61YsYJDhw6d2D58+DArVqzo4UTS3Iy1+tY73vEOfvazn/H0009z7Ngxdu7cyYYNG3o9ljSrJb0eQOqVJUuW8OUvf5nrrruOl19+mY9+9KO89a1v7fVY0qyMtfra+vXrWb9+fa/HkOblZRBJKsBYS1IBxlqSCjDWklSAsZakAoy1JBVgrCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVICxlqQCjLUkFWCsJakAYy1JBRhrSSrAWEtSAcZakgow1pJUgLGWpAKMtSQV0FWsI2JdRDwVEfsi4paTrHtfRGREjDY3oiRp3lhHxPnANuB6YC2wKSLWzrLuQuAvgb1NDylJ/a6bZ9ZXAvsyc39mHgN2AjfMsu5zwBeAow3OJ0miu1ivAA51bB9u7zshIq4AVmXmd072QBGxJSImImJienr6lIeVpH51xi8wRsR5wBeBT8+3NjPHMnM0M0cHBwfP9NCS1De6ifURYFXH9sr2vuMuBC4D7o+IA8DVwLgvMkpSc7qJ9UPAmohYHRFLgY3A+PE7M/P5zFyemSOZOQLsATZk5sSCTCxJfWjeWGfmS8DNwG7gSeCezHwiIm6PiA0LPaAkCZZ0sygzdwG7Zuy7bY6115z5WJKkTr6DUZIKMNaSVICxlqQCjLUkFWCsJakAYy1JBRhrSSrAWEtSAcZakgow1pJUgLGWpAKMtSQVYKwlqQBjLUkFGGtJKsBYS1IBxlqSCjDWklSAsZakAoy1JBVgrCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVICxlqQCjLUkFWCsJakAYy1JBRhrSSrAWEtSAcZakgow1pJUgLGWpAKMtSQVYKwlqQBjLUkFGGtJKsBYS1IBxlqSCjDWklSAsZakAoy1JBVgrCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVEBXsY6IdRHxVETsi4hbZrn/UxExGRGPRcT3ImK4+VElqX/NG+uIOB/YBlwPrAU2RcTaGcseAUYz8w+BbwJ/2/SgktTPunlmfSWwLzP3Z+YxYCdwQ+eCzLwvM3/b3twDrGx2TEnqb93EegVwqGP7cHvfXG4E7j2ToSRJr7SkyQeLiA8Co8B75rh/C7AFYGhoqMlDS9I5rZtn1keAVR3bK9v7XiEirgU+A2zIzBdne6DMHMvM0cwcHRwcPJ15JakvdRPrh4A1EbE6IpYCG4HxzgURcTnwVVqhfqb5MSWpv80b68x8CbgZ2A08CdyTmU9ExO0RsaG97O+A1wLfiIgfR8T4HA8nSToNXV2zzsxdwK4Z+27ruH1tw3NJkjr4DkZJKsBYS1IBxlqSCjDWklSAsZakAoy1JBVgrCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVICxlqQCjLUkFWCsJakAYy1JBRhrSSrAWEtSAcZakgow1pJUgLGWpAKMtSQVYKwlqQBjLUkFGGtJKsBYS1IBxlqSCjDWklSAsZakAoy1JBVgrCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVICxlqQCjLUkFWCsJakAYy1JBRhrSSrAWEtSAcZakgow1pJUgLGWpAKMtSQVYKwlqQBjLUkFGGtJKsBYS1IBxlqSCugq1hGxLiKeioh9EXHLLPdfEBF3t+/fGxEjjU8KbN++nb179/CDH/yAkZERtm/fvhCHkaRTttB9mjfWEXE+sA24HlgLbIqItTOW3Qg8l5mXAP8AfKHRKWn9RmzZsoWjR18EYGpqii1bthhsST23GH2KzDz5goh3Ap/NzOva27cCZOYdHWt2t9c8EBFLgP8CBvMkDz46OpoTExNdDzoyMsLU1BRv2nQHSy9azbFnngZg2bILuOqqq7t+HElq2t69ezh69EWWDb2Nowd/wq923ArA8PAwBw4c6PpxIuLhzByd7b5uLoOsAA51bB9u75t1TWa+BDwPvHGWQbZExERETExPT3cz+wkHDx4E4DeT958INXDiXzJJ6pXjHeoMNfx/t5qwpLFH6kJmjgFj0HpmfSq/dmhoiKmpKV54dDcvPLr7xP7h4WHuvutVl9EladGM3LGJqampV+0fGhpq7BjdPLM+Aqzq2F7Z3jfrmvZlkNcDzzYx4HFbt25lYGDgFfsGBgbYunVrk4eRpFO2GH3qJtYPAWsiYnVELAU2AuMz1owDH27ffj/w/ZNdrz4dmzdvZmxsjOHhYSKC4eFhxsbG2Lx5c5OHkaRTthh9mvcFRoCIWA/8I3A+8LXM3BoRtwMTmTkeEcuArwOXA/8NbMzM/Sd7zFN9gVGSznUne4Gxq2vWmbkL2DVj320dt48Cf3YmQ0qS5uY7GCWpAGMtSQUYa0kqwFhLUgHGWpIKMNaSVICxlqQCjLUkFWCsJamArt5uviAHjpgGXv1tqrqzHPh1g+OczfrpXKG/ztdzPTedybkOZ+bgbHf0LNZnIiIm5nr//Lmmn84V+ut8Pddz00Kdq5dBJKkAYy1JBVSN9VivB1hE/XSu0F/n67memxbkXEtes5akflP1mbUk9ZWzOtYRsS4inoqIfRHxqp+KGxEXRMTd7fv3RsRID8ZsRBfn+qmImIyIxyLiexEx3Is5mzDfuXase19EZESU/SqCbs41Ij7Q/tw+ERF3LfaMTeriz/FQRNwXEY+0/yyv78WcTYiIr0XEMxHx+Bz3R0R8qf178VhEXHFGB8zMs/KD1o8Q+zlwMbAUeBRYO2PNXwBfad/eCNzd67kX8Fz/GBho3/7EuXyu7XUXAj8E9gCjvZ57AT+va4BHgN9vb1/U67kX+HzHgE+0b68FDvR67jM433cDVwCPz3H/euBeIICrgb1ncryz+Zn1lcC+zNyfmceAncANM9bcAPxr+/Y3gfdGRCzijE2Z91wz877M/G17cw+tnzJfUTefV4DPAV8Aji7mcA3r5lw/BmzLzOcAMvOZRZ6xSd2cbwKva99+PfCLRZyvUZn5Q1o/c3YuNwD/li17gDdExB+c7vHO5livAA51bB9u75t1TWa+BDwPvHFRpmtWN+fa6UZa/2JXNO+5tv+7uCozv7OYgy2Abj6vlwKXRsSPImJPRKxbtOma1835fhb4YEQcpvVzXT+5OKP1xKn+vT6prn5grs4eEfFBYBR4T69nWQgRcR7wReAjPR5lsSyhdSnkGlr/W/phRLwtM/+nl0MtoE3Av2Tm30fEO4GvR8Rlmfm/vR7sbHc2P7M+Aqzq2F7Z3jfrmohYQuu/Vc8uynTN6uZciYhrgc8AGzLzxUWarWnzneuFwGXA/RFxgNa1vvGiLzJ283k9DIxn5u8y82ngp7TiXVE353sjcA9AZj4ALKP1vTTORV39ve7W2Rzrh4A1EbE6IpbSegFxfMaaceDD7dvvB76f7Sv7xcx7rhFxOfBVWqGufF3zpOeamc9n5vLMHMnMEVrX5zdk5kRvxj0j3fwZ/hatZ9VExHJal0X2L+KMTermfA8C7wWIiLfQivX0ok65eMaBD7W/KuRq4PnM/OVpP1qvX1Gd59XW9bSeafwc+Ex73+20/vJC6xP9DWAf8CBwca9nXsBz/U/gV8CP2x/jvZ55oc51xtr7KfrVIF1+XoPWZZ9J4CfAxl7PvMDnuxb4Ea2vFPkx8Ke9nvkMznUH8Evgd7T+h3Qj8HHg4x2f223t34ufnOmfY9/BKEkFnM2XQSRJbcZakgow1pJUgLGWpAKMtSQVYKwlqQBjLUkFGGtJKuD/ABzL8K3jzMJQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "grid_shapely = list(map(lambda x: shapely.geometry.box(*x), grid))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6,6))\n",
    "\n",
    "gpd.GeoDataFrame(geometry=grid_shapely).exterior.plot(ax=ax)\n",
    "\n",
    "for coos, index in zip(grid, grid_index):\n",
    "    ax.text(coos[[0,2]].mean(), coos[[1,3]].mean(), str(index), horizontalalignment='center', verticalalignment='center')\n",
    "    \n",
    "ax.scatter(test_points[:,0], test_points[:,1], color='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Intersect Points with Grid\n",
    "- The boolean parameters include_xmax and include_ymax determine which bounds are inclusive and exclusive. \n",
    "- If both are false the bottom left corner of the grid (origin) will be considered inclusive, and the other 3 corners exclusive.\n",
    "- For the other points that lie inside the rectangles the index will remain the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([                   0, -9223372036854775808, -9223372036854775808,\n",
       "       -9223372036854775808])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "intersect_points_grid(xy=test_points, \n",
    "                      xmin=0, ymin=0, xmax=1, ymax=1, xstep=1, ystep=1, \n",
    "                      include_xmax=False, include_ymax=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-9223372036854775808,                    0, -9223372036854775808,\n",
       "       -9223372036854775808])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "intersect_points_grid(xy=test_points, \n",
    "                      xmin=0, ymin=0, xmax=1, ymax=1, xstep=1, ystep=1, \n",
    "                      include_xmax=False, include_ymax=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-9223372036854775808, -9223372036854775808,                    0,\n",
       "       -9223372036854775808])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "intersect_points_grid(xy=test_points, \n",
    "                      xmin=0, ymin=0, xmax=1, ymax=1, xstep=1, ystep=1, \n",
    "                      include_xmax=True, include_ymax=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-9223372036854775808, -9223372036854775808, -9223372036854775808,\n",
       "                          0])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "intersect_points_grid(xy=test_points, \n",
    "                      xmin=0, ymin=0, xmax=1, ymax=1, xstep=1, ystep=1, \n",
    "                      include_xmax=True, include_ymax=True)"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}