cKDTree spatial example

cKDTree.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nearest neighbor searches on geodataframes with cKDTree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.spatial import cKDTree\n",
    "import geopandas as gpd\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Later in the notebook we'll load actual geospatial data in two geodataframes. One `schools` gdf will contain locations and column data on all 1822 NYC schools and another locations and data on 3019 street intersections with Leading Pedestrian Interval Signals (lpis). \n",
    "\n",
    "We want to find the nearest school to each of the 3019 lpis intersections.\n",
    "\n",
    "Right now we'll start with a very simulated small dataset of just N=25 just to visualize what's going on. We store these as the numpy arrays that cKDTree wants to work with, but we'll use `nschools` and `nlpis` names to play along with the example. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 25\n",
    "nschools = np.random.uniform(0., 5000., (N, 2))  \n",
    "nlpis = np.random.uniform(0., 5000., (N, 2))  #"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that we want to find the nearest neighbor to a specific LPIS intersection (let's choose one more or less in middle of the 'map')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lpispoint = np.array([2000,2000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Build the cKDTree**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree = cKDTree(nschools)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The nearest neighbors to lpispoint are found with a query to the tree:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(239.67420989365044, 24)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dist, idx = tree.query(lpispoint,k=1)\n",
    "dist, idx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This returns the distance to the nearest neighbor as well as the index position in the tree/array.\n",
    "\n",
    "We can now retrieve the nearest school location to lpispoint using the returned index position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2038.04675196, 2236.63510212])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nschools[idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a visualization of what's going on. Our lpispoint is the green circle and the nearest neighbor is the star."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD8CAYAAACVZ8iyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF9JJREFUeJzt3X2MXeV94PHvj2GA4S1jwKFmgLWTeL2hSxKzU0Lk1SpLGkxYVCxEVs624HaRLG3SVZJGTm11tSgvKiZuQpptkoYuURyaBihxjJWXdSwculJVXsYxwSHEtaE0eIywG3tgCQMM5rd/3GfsOzC25+XOnPvy/UhX95zffe6d53nse3/nPOc550RmIknqbCdUXQFJUvVMBpIkk4EkyWQgScJkIEnCZCBJwmQgScJkIEnCZCBJAk6sugLHcs455+T8+fOrroYktZRt27b9S2bOncx7mjoZzJ8/n4GBgaqrIUktJSL+ebLvcZhIkmQykCRNMBlExFMRsSMiHomIgRI7KyK2RMSu8jynxCMivhQRuyPi0Yi4pO5zVpTyuyJixcw0SZI0WZPZM/iPmfmuzOwv66uB+zJzIXBfWQf4ALCwPFYCX4Va8gBuAt4NXArcNJpAJEnVms4w0TXA+rK8HlhWF/9m1jwA9EbEPGApsCUzD2TmQWALcOU0/r4kqUEmOpsogR9FRAJfy8zbgHMz8xmAzHwmIt5cyvYBT9e9d0+JHS2uadq4fZB1m3eyd2iY83p7WLV0EcsW27WSJm6iyWBJZu4tP/hbIuIXxygb48TyGPGxb45YSW14iQsvvHCC1etcG7cPsmbDDoZHDgEwODTMmg07AEwIkiZsQsNEmbm3PO8DvkttzP/ZMvxDed5Xiu8BLqh7+/nA3mPEX/+3bsvM/szsnzt3UudMdKR1m3ceTgSjhkcOsW7zzopqJKkVHTcZRMRpEXHG6DJwBfAzYBMwOiNoBXBvWd4E3FBmFV0GPFeGkzYDV0TEnHLg+IoS0zTsHRqeVFySxjORYaJzge9GxGj5v8nM/xMRDwN3R8SNwC+BD5byPwCuAnYDLwJ/AJCZByLiM8DDpdynM/NAw1rSoc7r7WFwnB/+83p7KqiNpFZ13GSQmU8C7xwn/ivgfePEE/jIUT7r68DXJ19NHc2qpYvGHDMA6OnuYtXSRRXWSlKraeprE+n4Rg8SO5tI0nSYDNrAssV9/vhLmhavTSRJMhlIkkwGkiRMBpIkTAaSJEwGkiRMBpIkTAaSJEwGkiQ8A1kdzJsCSUeYDNSRvCmQNJbDROpI3hRIGstkoI7kTYGksUwG6khHu/mPNwVSpzIZqCOtWrqInu6uMTFvCqRO5gFkdSRvCiSNZTJQx/KmQNIRDhNJktwzkCbLk9XUjkwG0iR4spralcmggdxibH/HOlnNf2s1QlW/IyaDBnGLsTN4sppmUpW/Ix5AbhAvb9AZPFlNM6nK3xGTQYO4xdgZPFlNM6nK3xGTQYO4xdgZli3u4+ZrL6avt4cA+np7uPnaix0KVENU+TviMYMGWbV00ZixPnCLsV15sppmSpW/IyaDBvHyBpPjzCvpjar8HYnMnPE/MlX9/f05MDBQdTXUYK+fMQG1rR+HW6TGiIhtmdk/mfe4Z6BZ51x9Vcm90vGZDDTrnHmlqng+0NFNeDZRRHRFxPaI+F5ZXxARD0bEroi4KyJOKvGTy/ru8vr8us9YU+I7I2Jpoxuj1uDMK1XF84GObjJTSz8KPF63fgtwa2YuBA4CN5b4jcDBzHwbcGspR0RcBCwHfhO4EvhKRIydsK2O4Fx9VaXZ9ko3bh9kydqtLFj9fZas3crG7YOV1AMmmAwi4nzgPwH/u6wHcDlwTymyHlhWlq8p65TX31fKXwPcmZkvZ+Y/AbuBSxvRCLUW5+qrKs20Vzo6ZDU4NExyZMiqqoQw0WMGXwQ+CZxR1s8GhjLz1bK+Bxj9JvcBTwNk5qsR8Vwp3wc8UPeZ9e85LCJWAisBLrzwwgk3RK3FufqqQjOdD9RsEymOu2cQEVcD+zJzW314nKJ5nNeO9Z4jgczbMrM/M/vnzp17vOpJ0oQ1015psw1ZTWTPYAnwOxFxFXAKcCa1PYXeiDix7B2cD+wt5fcAFwB7IuJE4E3Agbr4qPr3SNKsaJa90vN6exgc54e/qokUx90zyMw1mXl+Zs6ndgB4a2b+LvBj4LpSbAVwb1neVNYpr2/N2pltm4DlZbbRAmAh8FDDWiJJLaTZJlJM5zyDPwbujIjPAtuB20v8duCOiNhNbY9gOUBmPhYRdwM/B14FPpKZh974sZLU/prtEjZejkJtwbNKpSO8HIU6kmeVStPn/QzU8jyrVJo+k4FaXrNN0ZNakclALa+ZziqVWpXJQC2v2aboSa3IA8hqec02RU9qRSYDtYVmOatUalUOE0mSTAaSJJOBJAmTgSQJk4EkCZOBJAmTgSQJk4EkCZOBJAmTgSQJk4EkCZOBJAmTgSQJr1oqNb2N2we9PLdmnMlAamIbtw+yZsOOw/d4HhwaZs2GHQAmBDWUw0RSE1u3eefhRDBqeOQQ6zbvrKhGalcmA6mJ7R0anlRcmiqTgdTEzuvtmVRcmiqTgdTEVi1dRE9315hYT3cXq5YuqqhGalceQJaa2OhBYmcTaaaZDKQmt2xxnz/+mnEOE0mSTAaSJJOBJIkJJIOIOCUiHoqIn0bEYxHxqRJfEBEPRsSuiLgrIk4q8ZPL+u7y+vy6z1pT4jsjYulMNUqSNDkT2TN4Gbg8M98JvAu4MiIuA24Bbs3MhcBB4MZS/kbgYGa+Dbi1lCMiLgKWA78JXAl8JSLGzpmTJFXiuMkga14oq93lkcDlwD0lvh5YVpavKeuU198XEVHid2bmy5n5T8Bu4NKGtEKSNC0TOmYQEV0R8QiwD9gCPAEMZearpcgeYHTuWx/wNEB5/Tng7Pr4OO+RJFVoQskgMw9l5ruA86ltzb99vGLlOY7y2tHiY0TEyogYiIiB/fv3T6R6kqRpmtRsoswcAu4HLgN6I2L0pLXzgb1leQ9wAUB5/U3Agfr4OO+p/xu3ZWZ/ZvbPnTt3MtWTJE3RRGYTzY2I3rLcA/w28DjwY+C6UmwFcG9Z3lTWKa9vzcws8eVlttECYCHwUKMaIkmauolcjmIesL7M/DkBuDszvxcRPwfujIjPAtuB20v524E7ImI3tT2C5QCZ+VhE3A38HHgV+EhmHkKSVLmobbQ3p/7+/hwYGKi6GpLUUiJiW2b2T+Y9noEsSTIZSJJMBpIkTAaSJEwGkiRMBpIkTAaSJEwGkiRMBpIkTAaSJEwGkiRMBpIkTAaSJEwGkiRMBpIkTAaSJEwGkiRMBpIkTAZqJa+9VnUNpLZlMlBr+NWv4IIL4MCBqmsitSWTgVrDpk2wd2/tWVLDmQx02MbtgyxZu5UFq7/PkrVb2bh9sOoqHfGNb4x9ltRQJ1ZdATWHjdsHWbNhB8MjhwAYHBpmzYYdACxb3Fdl1eD55+GBB2rLDzxQWz/zzGrrJLUZk0GL2rh9kHWbd7J3aJjzentYtXTRtH60123eeTgRjBoeOcS6zTtnNxkcPAi//OWY0I5vbuAt0cVpwK85gSc/9XkuvuHase+78EKYM2f26im1GZNBC5qJrfi9Q8OTis+YtWvhc5+DU06Bk05i5NBrzB95jdNeqdXjtJeHmf8Xf8bI175Ad9cJ8Mor8NJL8MlPwi23zG5dpTbiMYMWdKyt+Kk6r7dnUvEZc/PN8Kd/ChHw/PN0//oFznjlxTFFznjlRbp//UJtuCiiVv7mm2e3nlKbMRm0oJnYil+1dBE93V1jYj3dXaxaumjKnzklJ5wAa9bA/ffDb/wGL3V1j1vspa5umDcP/u7vauVP8L+yNB1+g1rQTGzFL1vcx83XXkxfbw8B9PX2cPO1F1d38PjSS2HnTg6cefa4Lx8482z4xS/gt35rlismtSePGbSgVUsXjTlmAI3Zil+2uK/6mUP1TjqJc188OO5L5744BCefPMsVktqXewYtqOm24mfKli10lR/84e6TeaXrRIa7a+tdJ58EW7ZUWTuprbhn0KKabiu+wZ448ARDt3yMS55/nhdPhG/2v8Y//tHvc9Pml+n51j21g8d33AFXX111VaW24J6Bms4Pd/2Qf/cXF/NvHnqS50+Ga5fDh5eO8OXH19O34Ds8/L/W1E46+/73YWSk6upKbeG4ySAiLoiIH0fE4xHxWER8tMTPiogtEbGrPM8p8YiIL0XE7oh4NCIuqfusFaX8rohYMXPNUqt64sATXPe31/HyK8N8453wr/87/OhttddGXhvhxZEXee+za3nqwR/BDTfAq69WW2GpTUxkz+BV4BOZ+XbgMuAjEXERsBq4LzMXAveVdYAPAAvLYyXwVaglD+Am4N3ApcBNowlEGvX5f/g8I4dGeKkb/vBq2Hf6G8uMHBrhz568A77yFeiZ5fMgpDZ13GSQmc9k5k/K8v8DHgf6gGuA9aXYemBZWb4G+GbWPAD0RsQ8YCmwJTMPZOZBYAtwZUNbo5b314/+NSOvHXvoZ+S1Ee549I5ZqpHUGSZ1zCAi5gOLgQeBczPzGaglDODNpVgf8HTd2/aU2NHir/8bKyNiICIG9u/fP5nqqQ288MoLDS0naWImnAwi4nTgO8DHMvP5YxUdJ5bHiI8NZN6Wmf2Z2T937tyJVk9t4vSTxhkXmkY5SRMzoWQQEd3UEsG3MnNDCT9bhn8oz/tKfA9wQd3bzwf2HiMuHfZ77/g9uk8Y/xIUo7pP6Ob6d1w/SzWSOsNEZhMFcDvweGZ+oe6lTcDojKAVwL118RvKrKLLgOfKMNJm4IqImFMOHF9RYtJhn3jPJ+g+yvWIRnV3dfPxyz4+SzWSOsNE9gyWANcDl0fEI+VxFbAWeH9E7ALeX9YBfgA8CewG/gr4MEBmHgA+AzxcHp8uMemwt571Vu754D2c2n3qG/YQuk/o5tTuU7nng/fw1rPeWlENpfYUmW8Ytm8a/f39OTAwUHU1VIEnDjzBrQ/cyh2P3sELr7zA6SedzvXvuJ6PX/ZxE4F0HBGxLTP7J/Uek4EktZepJIO2vDZRo28JqWp16r9np7Zb1Wi7ZNDUN3bXpHXqv2entlvVabsL1c3ELSFVnU799+zUdqs6bbdn0DQ3dldDdOq/Z1Xtdmiqc7XdnkHT3NhdDdGp/55VtHt0aGpwaJjkyNDUxu2DM/Y31TzaLhk0zY3d1RCd+u9ZRbsdmupsbTdMNLpL665ue+jUf88q2t2pQ3Kq8TwDSQAsWbuVwXF++Pt6e/j71ZdXUCNN1VTOM2i7YSJJU9OpQ3KqabthIklT06lDcqoxGUg6bNniPn/8O5TDRJIk9wwkqdFa8eQ9k4EkNVCrXlfKYSJJaqBWPXnPZCBJDdSqJ++ZDCSpgVr1elomA0lqoFY9ec8DyJLUQK168p7JQJIarBVP3nOYSJJkMpAkmQwkSXjMQNIEtOLlFTQ5JgNJx9Sql1fQ5DhMJOmYWvXyCpock4GkY2rVyytockwGko6pVS+voMkxGUg6pla9vIImxwPIko6pVS+voMk5bjKIiK8DVwP7MvPflthZwF3AfOAp4D9n5sGICODPgauAF4Hfz8yflPesAP5H+djPZub6xjZF0kxpxcsraHImMkz0DeDK18VWA/dl5kLgvrIO8AFgYXmsBL4Kh5PHTcC7gUuBmyJiznQrL0Ft6uOStVtZsPr7LFm7lY3bB6uuktRyjrtnkJn/NyLmvy58DfDesrweuB/44xL/ZmYm8EBE9EbEvFJ2S2YeAIiILdQSzLen3QJ1NOfA6/U8QW5qpnoA+dzMfAagPL+5xPuAp+vK7Smxo8WlaXEOvOqNbhwMDg2THNk4cG/x+Bo9myjGieUx4m/8gIiVETEQEQP79+9vaOXUfpwDr3puHEzdVJPBs2X4h/K8r8T3ABfUlTsf2HuM+Btk5m2Z2Z+Z/XPnzp1i9dQpnAOvem4cTN1Uk8EmYEVZXgHcWxe/IWouA54rw0ibgSsiYk45cHxFibU1D2zOPOfAq950Ng46/ft63GQQEd8G/gFYFBF7IuJGYC3w/ojYBby/rAP8AHgS2A38FfBhgHLg+DPAw+Xx6dGDye3KscvZsWxxHzdfezF9vT0E0Nfbw83XXuwBww411Y0Dv68QtYk/zam/vz8HBgaqrsaULFm7lcFxdk37env4+9WXV1AjqTNMZTZRu31fI2JbZvZP5j2egTxDHLuUqjGVE+T8vnptohnjgU2pdfh9NRnMGA9sSq3D76vDRDPGi3tJrcPvqweQJantTOUAssNEkiSTgSTJZCBJwmQgScJkIEnCZCBJwmQgScJkIEnCZCBJwmQgScJkIEnCZCBJwmQgScJkIEnCZCBJwmQgScJkIEnCZCBJwmQgScJkIEnCZCBJwmQgScJkIEnCZCBJwmQgSQJOrLoCkgSwcfsg6zbvZO/QMOf19rBq6SKWLe6rulozptnaO+t7BhFxZUTsjIjdEbF6tv++pOazcfsgazbsYHBomAQGh4ZZs2EHG7cPVl21GdGM7Z3VZBARXcCXgQ8AFwEfioiLZrMOkprPus07GR45NCY2PHKIdZt3VlSjmdWM7Z3tPYNLgd2Z+WRmvgLcCVwzy3WQ1GT2Dg1PKt7qmrG9s50M+oCn69b3lJikDnZeb8+k4q2uGds728kgxonlmAIRKyNiICIG9u/fP0vVklSlVUsX0dPdNSbW093FqqWLKqrRzGrG9s72bKI9wAV16+cDe+sLZOZtwG0A/f39YxKFpPY0OoummWbXzKRmbG9kzt7vbUScCPwj8D5gEHgY+C+Z+dh45fv7+3NgYGDW6idJ7SAitmVm/2TeM6t7Bpn5akT8IbAZ6AK+frREIEmaPbN+0llm/gD4wWz/XUnS0Xk5CkmSyUCSZDKQJDHLs4kmKyL2A/9cdT0qcg7wL1VXoknYFzX2wxH2Rc3R+uFfZebcyXxQUyeDThYRA5OdGtau7Isa++EI+6Kmkf3gMJEkyWQgSTIZNLPbqq5AE7EvauyHI+yLmob1g8cMJEnuGUiSTAazKiK+HhH7IuJndbGzImJLROwqz3NKPCLiS+X2oI9GxCV171lRyu+KiBVVtGU6IuKCiPhxRDweEY9FxEdLvKP6IiJOiYiHIuKnpR8+VeILIuLB0qa7IuKkEj+5rO8ur8+v+6w1Jb4zIpZW06Lpi4iuiNgeEd8r6x3XFxHxVETsiIhHImKgxGb+u5GZPmbpAfwH4BLgZ3WxzwGry/Jq4JayfBXwQ2r3gLgMeLDEzwKeLM9zyvKcqts2yX6YB1xSls+gdiXbizqtL0p7Ti/L3cCDpX13A8tL/C+B/1aWPwz8ZVleDtxVli8CfgqcDCwAngC6qm7fFPvkj4C/Ab5X1juuL4CngHNeF5vx70blDe+0BzD/dclgJzCvLM8DdpblrwEfen054EPA1+riY8q14gO4F3h/J/cFcCrwE+Dd1E4iOrHE3wNsLsubgfeU5RNLuQDWAGvqPutwuVZ6ULu/yX3A5cD3Sts6ri+Okgxm/LvhMFH1zs3MZwDK85tL/Gi3CG2rW4eW3fvF1LaKO64vyrDII8A+YAu1LdmhzHy1FKlv0+H2ltefA86mDfqh+CLwSeC1sn42ndkXCfwoIrZFxMoSm/HvxqxfwloTdrRbhB731qGtIiJOB74DfCwzn48Yr2m1ouPE2qIvMvMQ8K6I6AW+C7x9vGLluW37ISKuBvZl5raIeO9oeJyibd8XwJLM3BsRbwa2RMQvjlG2Yf3gnkH1no2IeQDleV+JH+0Woce9dWgriIhuaongW5m5oYQ7si8AMnMIuJ/auG9vuSsgjG3T4faW198EHKA9+mEJ8DsR8RRwJ7Whoi/SgX2RmXvL8z5qGwiXMgvfDZNB9TYBo0f6V1AbPx+N31BmC1wGPFd2DzcDV0TEnDKj4IoSaxlR2wW4HXg8M79Q91JH9UVEzC17BERED/DbwOPAj4HrSrHX98No/1wHbM3agPAmYHmZYbMAWAg8NDutaIzMXJOZ52fmfGoHhLdm5u/SYX0REadFxBmjy9T+T/+M2fhuVH2wpJMewLeBZ4ARapn7RmrjnPcBu8rzWaVsAF+mNoa8A+iv+5z/Cuwujz+oul1T6Id/T22X9VHgkfK4qtP6AngHsL30w8+A/1nib6H2A7Yb+Fvg5BI/pazvLq+/pe6z/qT0z07gA1W3bZr98l6OzCbqqL4o7f1peTwG/EmJz/h3wzOQJUkOE0mSTAaSJEwGkiRMBpIkTAaSJEwGkiRMBpIkTAaSJOD/A02I24qnMfD6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(nschools[:,0],nschools[:,1])\n",
    "plt.scatter(lpispoint[0],  lpispoint[1] ,s=100,c='g')\n",
    "plt.scatter(nschools[idx][0], nschools[idx][1], c='r', marker=\"*\",s=200);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The great thing about cKDTree queries is that they can be vectorized. \n",
    "\n",
    "Let's now raise the challenge and searh for nearest neighbors to 10000 locations for which we have 10000 possible neighbors. If we were doing this by brute force we'd have to do 100 million distance calculations which would be very costly.  The spatial index greatly reduces this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10000\n",
    "nschools = np.random.uniform(0., 5000., (N, 2))  \n",
    "nlpis = np.random.uniform(0., 5000., (N, 2))  #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree = cKDTree(nschools)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11.3 ms ± 85.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "dist, idx = tree.query(nlpis,k=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Real world geodataframe data\n",
    "\n",
    "Usually we'd have geospatial data stored in geopandas dataframes, so let's write a wrapper function to work with some real world data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>school_nam</th>\n",
       "      <th>ft_lon</th>\n",
       "      <th>ft_lat</th>\n",
       "      <th>ID</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>P.S. 015 Roberto Clemente</td>\n",
       "      <td>990141.0</td>\n",
       "      <td>202349.0</td>\n",
       "      <td>0</td>\n",
       "      <td>POINT (990141 202349)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>P.S. 019 Asher Levy</td>\n",
       "      <td>988547.0</td>\n",
       "      <td>205239.0</td>\n",
       "      <td>1</td>\n",
       "      <td>POINT (988547 205239)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  school_nam    ft_lon    ft_lat  ID               geometry\n",
       "0  P.S. 015 Roberto Clemente  990141.0  202349.0   0  POINT (990141 202349)\n",
       "1        P.S. 019 Asher Levy  988547.0  205239.0   1  POINT (988547 205239)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "schools = gpd.read_file('schools.shp')\n",
    "schools.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LPID</th>\n",
       "      <th>main_stree</th>\n",
       "      <th>cross_stre</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>Co-op City Blvd</td>\n",
       "      <td>Dreiser loop East</td>\n",
       "      <td>POINT (1031739.000154228 259373.000004255)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>Lenox Avenue</td>\n",
       "      <td>West 119 St</td>\n",
       "      <td>POINT (998570.9998893011 232184.9999184268)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   LPID       main_stree         cross_stre  \\\n",
       "0     0  Co-op City Blvd  Dreiser loop East   \n",
       "1     1     Lenox Avenue        West 119 St   \n",
       "\n",
       "                                      geometry  \n",
       "0   POINT (1031739.000154228 259373.000004255)  \n",
       "1  POINT (998570.9998893011 232184.9999184268)  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lpis = gpd.read_file('lpis.shp')\n",
    "lpis.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A `ckdnearest` function\n",
    "\n",
    "Let's write a wrapping function for all of this to make nearest neighbor queries easy using geodaframes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ckdnearest(gdA, gdB, bcol):\n",
    "    \"\"\"\n",
    "    This function takes geodataframes: `gdA` and `gdB` and \n",
    "    a column name `bcol`. Both dataframes are assumed to have a `geometry` column. \n",
    "    It finds the nearest neighbor from each location in `gdA` to a \n",
    "    nearest neighbor in `gdB`.  \n",
    "\n",
    "    It returns a two-column pandas dataframe with a 'distance' (here rounded to nearest foot)\n",
    "    and the value of the `bcol` in `gdB'  (e.g. 'school_name')\n",
    "    \"\"\"\n",
    "    \n",
    "    nA = np.array(list(zip(gdA.geometry.x, gdA.geometry.y)) )\n",
    "    nB = np.array(list(zip(gdB.geometry.x, gdB.geometry.y)) )\n",
    "    btree = cKDTree(nB)\n",
    "    dist, idx = btree.query(nA,k=1)\n",
    "    df = pd.DataFrame.from_dict({'distance': dist.astype(int),\n",
    "                             'bcol' : gdB.loc[idx, bcol].values })\n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So if we wanted to insert the \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LPID</th>\n",
       "      <th>main_stree</th>\n",
       "      <th>cross_stre</th>\n",
       "      <th>geometry</th>\n",
       "      <th>distance_to_school</th>\n",
       "      <th>nearest_school</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>Co-op City Blvd</td>\n",
       "      <td>Dreiser loop East</td>\n",
       "      <td>POINT (1031739.000154228 259373.000004255)</td>\n",
       "      <td>1699</td>\n",
       "      <td>Cornerstone Academy for Social Action</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>Lenox Avenue</td>\n",
       "      <td>West 119 St</td>\n",
       "      <td>POINT (998570.9998893011 232184.9999184268)</td>\n",
       "      <td>479</td>\n",
       "      <td>Success Academy Charter School - Harlem 1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   LPID       main_stree         cross_stre  \\\n",
       "0     0  Co-op City Blvd  Dreiser loop East   \n",
       "1     1     Lenox Avenue        West 119 St   \n",
       "\n",
       "                                      geometry  distance_to_school  \\\n",
       "0   POINT (1031739.000154228 259373.000004255)                1699   \n",
       "1  POINT (998570.9998893011 232184.9999184268)                 479   \n",
       "\n",
       "                              nearest_school  \n",
       "0      Cornerstone Academy for Social Action  \n",
       "1  Success Academy Charter School - Harlem 1  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lpis[['distance_to_school','nearest_school']] = ckdnearest(lpis, schools,'school_nam')\n",
    "lpis.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "105 ms ± 2.35 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "ckdnearest(lpis, schools,'school_nam')"
   ]
  }
 ],
 "metadata": {
  "gist": {
   "data": {
    "description": "cKDTree spatial example",
    "public": true
   },
   "id": ""
  },
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python [default]",
   "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.6"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}