1-processing/sample_stop_computation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Processing OpenTransit Bus Location Data - Stop Computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The purpose of this notebook is to present an example of working with the provided bus location data, and our current implementation of stop computation. We'll take a subset of this data (pulled only from buses on routes 1 and 14), along with geographical information about the stops on each route, to find every instance in our data set of a bus arriving at a stop along its route."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import requests\n",
    "\n",
    "from datetime import datetime, timedelta, timezone\n",
    "from geopy.distance import distance\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from typing import List, Union"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load the Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_json(filename = 'sample_routes_data_pst_15s.json'):\n",
    "    with open(filename, 'r') as f:\n",
    "        return json.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "bus_locations = load_json()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get All Stop and Bus IDs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to filter the location data so we're left with every instance of a bus arriving at a bus stop. We'll extract the location data for every stop on a given route (the current implementation was designed with finding these instances on a given route for a short time period). We do the same for bus locations on a given route."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def produce_stops(data: list, route: str) -> pd.DataFrame:\n",
    "    stops = pd.io.json.json_normalize(data,\n",
    "                                      record_path=['stops']) \\\n",
    "            .rename(columns={'lat': 'LAT',\n",
    "                             'lon': 'LON',\n",
    "                             'sid': 'SID'}) \\\n",
    "            .reindex(['SID', 'LAT', 'LON'], axis='columns')\n",
    "    \n",
    "    # obtain stop directions\n",
    "    stops['DID'] = stops['SID'].map({stop: direction['id']\n",
    "                                     for direction in requests\n",
    "                                                      .get(f\"http://restbus.info/api/agencies/sf-muni/routes/{route}\")\n",
    "                                                      .json()['directions']\n",
    "                                     for stop in direction['stops']})\n",
    "    \n",
    "    # remove stops that don't have an associated direction\n",
    "    stops = stops.dropna(axis='index', subset=['DID'])\n",
    "    request = requests.get(f\"http://restbus.info/api/agencies/sf-muni/routes/{route}\")\n",
    "\n",
    "    # obtain stop ordinals\n",
    "    stops['ORD'] = stops['SID'].map({stop_meta['id']: ordinal\n",
    "                                     for ordinal, stop_meta\n",
    "                                     in enumerate(request.json()['stops'])})\n",
    "    \n",
    "    return stops\n",
    "\n",
    "def produce_buses(data: list) -> pd.DataFrame:\n",
    "     return pd.io.json.json_normalize(data,\n",
    "                                      record_path=['routeStates', 'vehicles'],\n",
    "                                      meta=[['routeStates', 'vtime']]) \\\n",
    "            .rename(columns={'lat': 'LAT',\n",
    "                             'lon': 'LON',\n",
    "                             'vid': 'VID',\n",
    "                             'did': 'DID',\n",
    "                             'routeStates.vtime': 'TIME'}) \\\n",
    "            .reindex(['TIME', 'VID', 'LAT', 'LON', 'DID'], axis='columns')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Find Stops"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to extract the arrivals from the location data, we find local minima for the distance, over time, of a given bus from a given stop. If a bus gets close enough to a given stop, we consider that bus as arriving at that particular stop.\n",
    "\n",
    "There are two implementations presented for computing the distance between a bus and a stop - one uses the `geopy` library to compute the geodesic distance, and the other uses the haversine formula to compute the distance. While the former is more accurate, it can be prohibitively slow to use on large portions of the data set. The latter is faster, but less accurate. The code below uses `geopy`, but can be rewritten to use the haversine formula instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# haversine formula for calcuating distance between two coordinates in lat lon\n",
    "# from bird eye view; seems to be +- 8 meters difference from geopy distance\n",
    "def haver_distance(latstop,lonstop,latbus,lonbus):\n",
    "\n",
    "    latstop,lonstop,latbus,lonbus = map(np.deg2rad,[latstop,lonstop,latbus,lonbus])\n",
    "    eradius = 6371000\n",
    "    \n",
    "    latdiff = (latbus-latstop)\n",
    "    londiff = (lonbus-lonstop)\n",
    "    \n",
    "    a = np.sin(latdiff/2)**2 + np.cos(latstop)*np.cos(latbus)*np.sin(londiff/2)**2\n",
    "    c = 2*np.arctan2(np.sqrt(a),np.sqrt(1-a))\n",
    "    \n",
    "    distance = eradius*c\n",
    "    return distance\n",
    "\n",
    "def find_eclipses(buses, stop):\n",
    "    \"\"\"\n",
    "    Find movement of buses relative to the stop, in distance as a function of time.\n",
    "    \"\"\"\n",
    "    def split_eclipses(eclipses, threshold=30*60*1000) -> List[pd.DataFrame]:\n",
    "        \"\"\"\n",
    "        Split buses' movements when they return to a stop after completing the route.\n",
    "        \"\"\"\n",
    "        disjoint_eclipses = []\n",
    "        for bus_id in eclipses['VID'].unique(): # list of unique VID's\n",
    "            # obtain distance data for this one bus\n",
    "            bus = eclipses[eclipses['VID'] == bus_id].sort_values('TIME')\n",
    "            #pprint.pprint(bus)\n",
    "            #pprint.pprint(bus['TIME'].shift())\n",
    "            #pprint.pprint(bus['TIME'].shift() + threshold)\n",
    "            #print('===============')\n",
    "            # split data into groups when there is at least a `threshold`-ms gap between data points\n",
    "            group_ids = (bus['TIME'] > (bus['TIME'].shift() + threshold)).cumsum()\n",
    "\n",
    "            # store groups\n",
    "            for _, group in bus.groupby(group_ids):\n",
    "                disjoint_eclipses.append(group)\n",
    "        return disjoint_eclipses\n",
    "\n",
    "    eclipses = buses.copy()\n",
    "    #eclipses['DIST'] = eclipses.apply(lambda row: distance(stop[['LAT','LON']],row[['LAT','LON']]).meters,axis=1)\n",
    "    \n",
    "    stopcord = stop[['LAT', 'LON']]\n",
    "    buscord = eclipses[['LAT', 'LON']]\n",
    "\n",
    "    # calculate distances fast with haversine function \n",
    "    eclipses['DIST'] = haver_distance(stopcord['LAT'],stopcord['LON'],buscord['LAT'],buscord['LON'])\n",
    "    # only keep positions within 750 meters within the given stop; (filtering out)\n",
    "    eclipses = eclipses[eclipses['DIST'] < 750]\n",
    "    \n",
    "    # update the coordinates list \n",
    "    stopcord = stop[['LAT', 'LON']].values\n",
    "    buscord = eclipses[['LAT', 'LON']].values\n",
    "    \n",
    "    # calculate distances again using geopy for the distance<750m values, because geopy is probably more accurate\n",
    "    dfromstop = []\n",
    "    for row in buscord:\n",
    "        busdistance = distance(stopcord,row).meters\n",
    "        dfromstop.append(busdistance)\n",
    "    eclipses['DIST'] = dfromstop\n",
    "    \n",
    "    # for haversine function:\n",
    "    #stopcord = stop[['LAT', 'LON']]\n",
    "    #buscord = eclipses[['LAT', 'LON']]\n",
    "    #eclipses['DIST'] = haver_distance(stopcord['LAT'],stopcord['LON'],buscord['LAT'],buscord['LON'])\n",
    "    \n",
    "    eclipses['TIME'] = eclipses['TIME'].astype(np.int64)\n",
    "    eclipses = eclipses[['TIME', 'VID', 'DIST']]\n",
    "    \n",
    "    eclipses = split_eclipses(eclipses)\n",
    "    \n",
    "    return eclipses\n",
    "\n",
    "def find_nadirs(eclipses):\n",
    "    \"\"\"\n",
    "    Find points where buses are considered to have encountered the stop.\n",
    "    \n",
    "    Nadir is an astronomical term that describes the lowest point reached by an orbiting body.\n",
    "    \"\"\"\n",
    "    def calc_nadir(eclipse: pd.DataFrame) -> Union[pd.Series, None]:\n",
    "        nadir = eclipse.iloc[eclipse['DIST'].values.argmin()]\n",
    "        if nadir['DIST'] < 100:  # if min dist < 100, then reasonable candidate for nadir\n",
    "            return nadir\n",
    "        else:  # otherwise, hardcore datasci is needed\n",
    "            rev_eclipse = eclipse.iloc[::-1]\n",
    "            rev_nadir = rev_eclipse.iloc[rev_eclipse['DIST'].values.argmin()]\n",
    "            if nadir['TIME'] == rev_nadir['TIME']:  # if eclipse has a global min\n",
    "                return nadir  # then it's the best candidate for nadir\n",
    "            else:  # if eclipse's min occurs at two times\n",
    "                mid_nadir = nadir.copy()\n",
    "                mid_nadir['DIST'] = (nadir['DIST'] + rev_nadir['DIST'])/2\n",
    "                return mid_nadir  # take the midpoint of earliest and latest mins\n",
    "    \n",
    "    nadirs = []\n",
    "    for eclipse in eclipses:\n",
    "        nadirs.append(calc_nadir(eclipse)[['VID', 'TIME']])\n",
    "        \n",
    "    return pd.DataFrame(nadirs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get all stops from a given data set\n",
    "# getting all stops at once from the entire set of location data might take prohibitively long (~3-4 hours)\n",
    "def get_all_stops(data):\n",
    "    bus_stops = pd.DataFrame(columns = [\"VID\", \"TIME\", \"SID\", \"DID\", \"ROUTE\"])\n",
    "    \n",
    "    for route in {ele['rid'] for ele in data}:\n",
    "        print(f\"{datetime.now().strftime('%a %b %d %I:%M:%S %p')}: Starting with {route}.\")\n",
    "        try:\n",
    "            stop_ids = [stop['id']\n",
    "                for stop\n",
    "                in requests.get(f\"http://restbus.info/api/agencies/sf-muni/routes/{route}\").json()['stops']]\n",
    "                 \n",
    "            route_data = [ele for ele in data if ele['rid'] == route]\n",
    "\n",
    "            for stop_id in stop_ids:\n",
    "                try:\n",
    "                    stops = produce_stops(route_data, route)\n",
    "                except:\n",
    "                    print(f\"{datetime.now().strftime('%a %b %d %I:%M:%S %p')}: could not produce stops df for {stop_id} on route {route}. Skipping.\")\n",
    "                    break\n",
    "\n",
    "                buses = produce_buses(route_data)\n",
    "\n",
    "                stop = stops[stops['SID'] == stop_id].drop_duplicates().squeeze()\n",
    "\n",
    "                try:\n",
    "                    buses = buses[buses['DID'] == stop['DID']]\n",
    "                except ValueError as err: # accounts for stops with no associated direction\n",
    "                    print(f\"{datetime.now().strftime('%a %b %d %I:%M:%S %p')}: skipping buses for {stop_id} and {route} due to ValueError: {err}\")\n",
    "                    continue\n",
    "\n",
    "                eclipses = find_eclipses(buses, stop)\n",
    "                nadirs = find_nadirs(eclipses)\n",
    "\n",
    "                try:\n",
    "                    nadirs[\"TIME\"] = nadirs[\"TIME\"].apply(lambda x: datetime.fromtimestamp(x//1000, timezone(timedelta(hours = -8))))\n",
    "                    nadirs[\"SID\"] = stop_id\n",
    "                    nadirs[\"DID\"] = stop[\"DID\"]\n",
    "                    nadirs[\"ROUTE\"] = route\n",
    "                    bus_stops = bus_stops.append(nadirs, sort = True)\n",
    "                except:\n",
    "                    print(f\"{datetime.now().strftime('%a %b %d %I:%M:%S %p')}: could not produce stops df for {stop_id} on route {route}. Skipping.\", end = \"\\r\")              \n",
    "        except KeyError:\n",
    "            print(f\"{datetime.now().strftime('%a %b %d %I:%M:%S %p')}: KeyError at {route}!\")\n",
    "            continue\n",
    "  \n",
    "    return bus_stops"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sun Mar 03 11:58:42 AM: Starting with 14.\n",
      "Sun Mar 03 12:00:47 PM: skipping buses for 37742 and 14 due to ValueError: Can only compare identically-labeled Series objects\n",
      "Sun Mar 03 12:01:00 PM: skipping buses for 7742 and 14 due to ValueError: Can only compare identically-labeled Series objects\n",
      "Sun Mar 03 12:02:47 PM: Starting with 1.\n"
     ]
    }
   ],
   "source": [
    "sample_stops = get_all_stops(bus_locations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`get_all_stops` produces a DataFrame containing every instance of a bus arriving at a stop from the given data. Columns of note:\n",
    "\n",
    "- `DID`: the direction ID of the bus; the O/I in the string is the most relevant part of this value, as it indicates whether the bus is outbound or inbound.\n",
    "\n",
    "- `SID`: the ID of the stop at which the bus arrived.\n",
    "\n",
    "- `VID`: the vehicle ID of the bus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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>DID</th>\n",
       "      <th>ROUTE</th>\n",
       "      <th>SID</th>\n",
       "      <th>TIME</th>\n",
       "      <th>VID</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>79014</th>\n",
       "      <td>1____O_F00</td>\n",
       "      <td>1</td>\n",
       "      <td>6311</td>\n",
       "      <td>2018-10-15 22:54:43-08:00</td>\n",
       "      <td>5633</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39052</th>\n",
       "      <td>1____I_F00</td>\n",
       "      <td>1</td>\n",
       "      <td>6310</td>\n",
       "      <td>2018-10-15 12:12:45-08:00</td>\n",
       "      <td>5596</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31858</th>\n",
       "      <td>1____I_F00</td>\n",
       "      <td>1</td>\n",
       "      <td>3827</td>\n",
       "      <td>2018-10-15 10:53:29-08:00</td>\n",
       "      <td>5545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>58465</th>\n",
       "      <td>14___O_F00</td>\n",
       "      <td>14</td>\n",
       "      <td>5566</td>\n",
       "      <td>2018-10-15 18:19:53-08:00</td>\n",
       "      <td>7228</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45812</th>\n",
       "      <td>1____I_F00</td>\n",
       "      <td>1</td>\n",
       "      <td>4025</td>\n",
       "      <td>2018-10-15 13:26:32-08:00</td>\n",
       "      <td>5545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59382</th>\n",
       "      <td>14___O_F00</td>\n",
       "      <td>14</td>\n",
       "      <td>5536</td>\n",
       "      <td>2018-10-15 18:34:08-08:00</td>\n",
       "      <td>7202</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>38412</th>\n",
       "      <td>1____O_F00</td>\n",
       "      <td>1</td>\n",
       "      <td>4015</td>\n",
       "      <td>2018-10-15 12:05:45-08:00</td>\n",
       "      <td>5556</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11200</th>\n",
       "      <td>1____I_F00</td>\n",
       "      <td>1</td>\n",
       "      <td>4026</td>\n",
       "      <td>2018-10-15 07:03:24-08:00</td>\n",
       "      <td>5553</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55589</th>\n",
       "      <td>1____I_F00</td>\n",
       "      <td>1</td>\n",
       "      <td>3846</td>\n",
       "      <td>2018-10-15 15:16:04-08:00</td>\n",
       "      <td>5617</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41420</th>\n",
       "      <td>1____I_F00</td>\n",
       "      <td>1</td>\n",
       "      <td>4016</td>\n",
       "      <td>2018-10-15 12:38:31-08:00</td>\n",
       "      <td>5553</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              DID ROUTE   SID                       TIME   VID\n",
       "79014  1____O_F00     1  6311  2018-10-15 22:54:43-08:00  5633\n",
       "39052  1____I_F00     1  6310  2018-10-15 12:12:45-08:00  5596\n",
       "31858  1____I_F00     1  3827  2018-10-15 10:53:29-08:00  5545\n",
       "58465  14___O_F00    14  5566  2018-10-15 18:19:53-08:00  7228\n",
       "45812  1____I_F00     1  4025  2018-10-15 13:26:32-08:00  5545\n",
       "59382  14___O_F00    14  5536  2018-10-15 18:34:08-08:00  7202\n",
       "38412  1____O_F00     1  4015  2018-10-15 12:05:45-08:00  5556\n",
       "11200  1____I_F00     1  4026  2018-10-15 07:03:24-08:00  5553\n",
       "55589  1____I_F00     1  3846  2018-10-15 15:16:04-08:00  5617\n",
       "41420  1____I_F00     1  4016  2018-10-15 12:38:31-08:00  5553"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_stops.sample(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reset index to get rid of duplicate indices, then drop redundant columns\n",
    "sample_stops = sample_stops.reset_index().drop(labels = 'index', axis = 'columns')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 24517 entries, 0 to 24516\n",
      "Data columns (total 5 columns):\n",
      "DID      24517 non-null object\n",
      "ROUTE    24517 non-null object\n",
      "SID      24517 non-null object\n",
      "TIME     24517 non-null object\n",
      "VID      24517 non-null object\n",
      "dtypes: object(5)\n",
      "memory usage: 8.5 MB\n"
     ]
    }
   ],
   "source": [
    "sample_stops.info(memory_usage = 'deep')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('sample_routes_stops_pst_15s.json', 'w') as f:\n",
    "    sample_stops.reset_index().to_json(f)"
   ]
  },
  {
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}