03-parcels-attach-blocks-network-nodes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Attach block ID and nearest network node to each parcel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.spatial import cKDTree\n",
    "from sklearn.neighbors import BallTree\n",
    "\n",
    "crs = {'init':'epsg:4326'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the input data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1956207\n",
      "Wall time: 7min 14s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# load parcels shapefile, drop any that lack geometry, and project to 4326\n",
    "parcels = gpd.read_file('data/Parcels/')\n",
    "parcels = parcels.dropna(subset=['geometry'])\n",
    "parcels = parcels.to_crs(crs)\n",
    "assert parcels.index.is_unique\n",
    "print(len(parcels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1956207\n"
     ]
    }
   ],
   "source": [
    "# load parcels-to-blocks lookup table\n",
    "parcels_blocks = pd.read_csv('data/parcels_joined_blocks.csv', dtype={'block_geoid':str})\n",
    "parcels_blocks = parcels_blocks.set_index('parcel_id')\n",
    "assert parcels_blocks.index.is_unique\n",
    "print(len(parcels_blocks))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "223328\n"
     ]
    }
   ],
   "source": [
    "# load the network nodes\n",
    "nodes = pd.read_csv('data/network/bay_area_full_strongly_nodes.csv')\n",
    "nodes = nodes.set_index('osmid')\n",
    "assert nodes.index.is_unique\n",
    "print(len(nodes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add census block ID to each parcel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert integer columns to int\n",
    "to_int = ['PARCEL_ID', 'DEVELOPMEN', 'COUNTY_ID', 'ZONE_ID', 'PROPORTION', 'TAX_EXEMPT', 'ID']\n",
    "for col in to_int:\n",
    "    parcels[col] = parcels[col].astype(int)\n",
    "\n",
    "# drop unused columns and set index\n",
    "parcels = parcels.drop(columns=['CENTROID', 'X', 'Y'])\n",
    "parcels = parcels.set_index('PARCEL_ID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge block IDs into parcels df\n",
    "parcels = pd.merge(left=parcels, right=parcels_blocks, how='left', left_index=True, right_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find network node nearest to each parcel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate parcel centroids\n",
    "parcels['centroid'] = parcels.centroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract x and y from centroids\n",
    "parcels = parcels.set_geometry('centroid')\n",
    "parcels['x'] = parcels.geometry.x\n",
    "parcels['y'] = parcels.geometry.y\n",
    "parcels = parcels.set_geometry('geometry')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### scipy euclidean method\n",
    "\n",
    "faster but less accurate with unprojected coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build the tree for nearest-neighbor search\n",
    "tree = cKDTree(data=nodes[['x', 'y']], compact_nodes=True, balanced_tree=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 2.88 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# query the tree for node nearest to each parcel centroid\n",
    "dist, idx = tree.query(parcels[['x', 'y']], k=1)\n",
    "parcels['node'] = nodes.iloc[idx].index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### scikit-learn haversine method\n",
    "\n",
    "slower but more accurate with unprojected coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# haversine requires data in form of [lat, lng] and inputs/outputs in units of radians\n",
    "nodes_rad = np.deg2rad(nodes[['y', 'x']])\n",
    "parcels_rad = np.deg2rad(parcels[['y', 'x']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build the tree for nearest-neighbor search\n",
    "tree = BallTree(nodes_rad, metric='haversine')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 10min\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# query the tree for node nearest to each parcel centroid\n",
    "idx = tree.query(parcels_rad, k=1, return_distance=False)\n",
    "parcels['node'] = nodes.iloc[idx[:,0]].index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save parcels data to disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# as simple csv\n",
    "df_save = parcels.reset_index()[['PARCEL_ID', 'block_geoid', 'node']]\n",
    "df_save = df_save.rename(columns={'PARCEL_ID':'parcel_id',\n",
    "                                  'block_geoid':'block_id',\n",
    "                                  'node':'node_id'})\n",
    "df_save.to_csv('data/parcels_blocks_nodes.csv', index=False, encoding='utf-8')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 27min 32s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# as full shapefile\n",
    "parcels = parcels.reset_index()\n",
    "parcels = parcels.drop(columns=['centroid', 'x', 'y'])\n",
    "parcels['node'] = parcels['node'].astype(str) # to handle ints bigger than 32-bits\n",
    "parcels.to_file('data/parcels_blocks_nodes')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}