step2_py.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Consolidating Route Patterns in GTFS\n",
    "\n",
    "Intended to be run after running the R script that filters the trips, schedules, etc. by a chosen date.\n",
    "\n",
    "The goal of this consolidation step is to do the following:\n",
    "    - All trips with same headsign/route/direction are assigned the same route_pattern_id. The route_pattern_id that represents the group is the route_pattern_id with the most trips per day.\n",
    "    - If a route pattern has less than 3 trips per tod period across all day, assign the route_pattern_id of the dir/route group with the most trips.\n",
    "\n",
    "As a consequence of this consolidation, several connected parts must also be updated: shape_id (on trips table), stop sequence (stops per trip, stop times per trip - e.g. stop_times table), and service_id, trip_headsign, and block_id (on trips table) if the trip was part of a route pattern/RHD that had < 3 trips in each time period.\n",
    "\n",
    "Notes:\n",
    "- The MBTA service day is 3AM - 26:59. For TOD assignment, everything that isn't between the hours of 6:30 - 19:00 is counted as NT. This addresses this issue as it includes 19:00-26:39 AND 3:00 - 6:30.\n",
    "\n",
    "\n",
    "Needs/Steps:\n",
    "- Number of trips per time period per route_pattern_id\n",
    "    - Midpoint time of each trip\n",
    "    - Each trip classified by TOD (based on midpoint)\n",
    "    - Sum of trips per TOD by route_pattern_id\n",
    "- route_pattern_id with most daily trips per Route, Direction, Headsign\n",
    "    - Sum all tod trips per route_pattern_id\n",
    "    - grab just the max per Route, Direction, Headsign (but keep route_pattern_id)\n",
    "- consolidate route patterns by Route, Direction, Headsign\n",
    "    - if route_pattern_id has less than 3 trips in each of the 4 TODs, replace with max trips route_pattern_id\n",
    "- once consolidated, update trips attributes to match new route pattern\n",
    "- update stop sequence and stop times for updated trips\n",
    "\n",
    "*** This script takes about 30 minutes to run with the 2018 MBTA Fall Recap file. 3 fxns take about 10 minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import datetime\n",
    "import copy\n",
    "\n",
    "import matplotlib\n",
    "matplotlib.use('agg')  # allows notebook to be tested in Travis\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy\n",
    "import matplotlib.pyplot as plt\n",
    "import pandana as pdna\n",
    "import time\n",
    "\n",
    "import urbanaccess as ua\n",
    "from urbanaccess.config import settings\n",
    "from urbanaccess.gtfsfeeds import feeds\n",
    "from urbanaccess import gtfsfeeds\n",
    "from urbanaccess.gtfs.gtfsfeeds_dataframe import gtfsfeeds_dfs\n",
    "from urbanaccess.network import ua_network, load_network\n",
    "\n",
    "from tqdm import tqdm\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# required bbox including all of Massachusetts and RI as well as parts of NH, CT, NY\n",
    "bbox = (-73.7207, 41.1198, -69.7876, 43.1161)\n",
    "# path to the downloaded and cleaned gtfs - mbta recap file for fall 2018\n",
    "#   this could also be a folder of gtfs folders (pre merge of multiple gtfs)\n",
    "\n",
    "path_to_gtfs = r\"J:\\Shared drives\\TMD_TSA\\Model\\networks\\Transit\\gtfs\\bat_2019\\1_gtfs_r\" #r\"J:\\Shared drives\\TMD_TSA\\Model\\networks\\Transit\\gtfs\\bnrd\\1_gtfs_r\"\n",
    "out_path = r\"J:\\Shared drives\\TMD_TSA\\Model\\networks\\Transit\\gtfs\\bat_2019\\2_gtfs_py\" #r\"J:\\Shared drives\\TMD_TSA\\Model\\networks\\Transit\\gtfs\\bnrd\\2_gtfs_py\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Checking GTFS text file header whitespace... Reading files using encoding: utf-8 set in configuration.\n",
      "GTFS text file header whitespace check completed. Took 0.05 seconds\n",
      "--------------------------------\n",
      "Processing GTFS feed: 1_gtfs_r\n",
      "The unique agency id: brockton_area_transit_authority was generated using the name of the agency in the agency.txt file.\n",
      "Unique agency id operation complete. Took 0.00 seconds\n",
      "Unique GTFS feed id operation complete. Took 0.00 seconds\n",
      "No GTFS feed stops were found to be outside the bounding box coordinates\n",
      "1_gtfs_r GTFS feed stops: coordinates are in northwest hemisphere. Latitude = North (90); Longitude = West (-90).\n",
      "Appended route type to stops\n",
      "Appended route type to stop_times\n",
      "--------------------------------\n",
      "Added descriptive definitions to stops, routes, stop_times, and trips tables\n",
      "Successfully converted ['departure_time'] to seconds past midnight and appended new columns to stop_times. Took 0.03 seconds\n",
      "1 GTFS feed file(s) successfully read as dataframes:\n",
      "     1_gtfs_r\n",
      "     Took 0.21 seconds\n"
     ]
    }
   ],
   "source": [
    "loaded_feeds = ua.gtfs.load.gtfsfeed_to_df(gtfsfeed_path= path_to_gtfs,\n",
    "                                           validation=True,\n",
    "                                           verbose=True,\n",
    "                                           bbox=bbox,\n",
    "                                           remove_stops_outsidebbox=False,\n",
    "                                           append_definitions=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions Part 1 : Admin\n",
    "\n",
    "1. Separate out Bus Routes\n",
    "2. Get start/stop times per trip\n",
    "3. Assign trip to TOD based on stop/start times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# separate out the bus routes and non-bus routes if multiple mode types\n",
    "def separate_bus(route_table,trip_table,x):\n",
    "    ''' filter out trips by mode to prevent from being consolidated\n",
    "\n",
    "    Parameters\n",
    "    -----------\n",
    "    route_table : df\n",
    "        gtfs routes.txt file in dataframe format\n",
    "    trip_table : df\n",
    "        gtfs trips.txt file in dataframe format\n",
    "    x : int\n",
    "        mode number to keep\n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "    non_bus_trips : df\n",
    "        df of trips to not consolidate\n",
    "    bus_trips : df\n",
    "        df of trips to consolidate\n",
    "    '''\n",
    "    non_bus_routes = list(route_table.query('route_type != @x').route_id.unique())\n",
    "    non_bus_trips = trip_table.query('route_id in @non_bus_routes')\n",
    "    bus_trips = trip_table.query('route_id not in @non_bus_routes')\n",
    "\n",
    "    if len(non_bus_trips) == 0:\n",
    "        print(\"All routes in this GTFS are bus routes.\")\n",
    "    if len(bus_trips) == 0:\n",
    "        print(\"There are no bus routes in this GTFS.\")\n",
    "    return(non_bus_trips, bus_trips)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_start_stop_times(stop_times):    \n",
    "    '''for every trip, grab the start time and stop time of the trip\n",
    "    \n",
    "    Parameters\n",
    "    -----------\n",
    "    stop_times : df\n",
    "        gtfs stop_times.txt in df format\n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "        flintstone : df\n",
    "            df with start and stop times per trip\n",
    "\n",
    "    '''\n",
    "    chocula =0 \n",
    "    for trip_id in stop_times['trip_id'].unique():\n",
    "        max_row = stop_times.query('trip_id==@trip_id').query('stop_sequence == stop_sequence.max()')[['trip_id','arrival_time']]\n",
    "        min_row = stop_times.query('trip_id==@trip_id').query('stop_sequence == stop_sequence.min()')[['trip_id','arrival_time']]\n",
    "        r2 = min_row.merge(max_row, how='left', on='trip_id', suffixes = ('_start','_end'))\n",
    "        if chocula == 0:\n",
    "            flintstone = pd.DataFrame(r2)\n",
    "        else:\n",
    "            flintstone=pd.concat([flintstone,r2])\n",
    "        chocula +=1\n",
    "    return(flintstone)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def assign_tod(start_stop):\n",
    "    '''calculate midpoint of trip, use midpoint to assign TOD\n",
    "    \n",
    "    Parameters\n",
    "    -----------\n",
    "    start_stop : df\n",
    "        df with start and stop times per trip\n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "    start_stop :\n",
    "        df with start time, stop time, midpoint time, and TOD per trip\n",
    "\n",
    "    '''\n",
    "    \n",
    "    start_stop['at_end_dec'] = (\n",
    "        (\n",
    "            (start_stop['arrival_time_end'].str.split(\":\").str[0]).astype('int32')\n",
    "            +\n",
    "            ((start_stop['arrival_time_end'].str.split(\":\").str[1]\n",
    "            ).astype('int32')/60)))\n",
    "    start_stop['at_start_dec'] = (\n",
    "        (\n",
    "            (start_stop['arrival_time_start'].str.split(\":\").str[0]).astype('int32')\n",
    "            +\n",
    "            ((start_stop['arrival_time_start'].str.split(\":\").str[1]\n",
    "            ).astype('int32')/60)))\n",
    "    \n",
    "    start_stop['midpoint'] = start_stop['at_start_dec'] + ((start_stop['at_end_dec']-start_stop['at_start_dec'])/2)\n",
    "    start_stop['tod'] = np.where(start_stop['midpoint'].between(6.50,9.50),'AM', np.where(\n",
    "        start_stop['midpoint'].between(9.50,15.00), 'MD', np.where(\n",
    "            start_stop['midpoint'].between(15.00,19.00),'PM', 'NT' \n",
    "        )\n",
    "            ) \n",
    "        )\n",
    "    \n",
    "    return start_stop\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def assign_trips_tod(trips, start_stop_tod):\n",
    "    \"\"\"\n",
    "    add TOD column from assign_tod() output to trips table\n",
    "        accomodates gtfs without route_pattern_ids or trip_headsigns for buses (e.g. coming out of remix)\n",
    "\n",
    "    Parameters\n",
    "    -----------\n",
    "    trips : df\n",
    "        gtfs trips.txt in df form\n",
    "    start_stop_tod : df\n",
    "        trips with tod by midpoint\n",
    "            (output of assign_tod)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    tod_trips : df\n",
    "        trips with a tod column\n",
    "    count_trips : df\n",
    "        number of trips per route_id, direction_id, and tod combination\n",
    "        - purpose is to be a useful summary file and to join back to imported gtfs\n",
    "\n",
    "    \"\"\"\n",
    "    if len(trips.groupby(by=['route_id','direction_id']).agg({'shape_id' : 'nunique'}).query('shape_id > 1')) > 0:\n",
    "        if (len(trips.query('~ trip_headsign.isna()')) > 0) & ('route_pattern_id' in trips.columns.to_list()):\n",
    "            print(\"Please run the functions that follows - you likely need to consolidate your route patterns\")\n",
    "        else: \n",
    "            print(\"Investigate your data more, you have multiple shapes per RHD but no route_pattern_id\")\n",
    "    else: \n",
    "        print(\"Skip to Clean Up & Export!\")\n",
    "\n",
    "    # add TOD to trip table\n",
    "    tod_trips = trips.merge(start_stop_tod[['trip_id','tod']], how='left', on='trip_id')\n",
    "    count_trips = tod_trips.groupby(by=['route_id','direction_id','tod']).count()\n",
    "    return(tod_trips, count_trips)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions Part 2 : Consolidate by RHD and RD\n",
    "1. Determine route_pattern_id with Max Number of Daily Trips\n",
    "    - Getting duplicate ids where same number of trips for max, so choosing one arbitrarily.\n",
    "    - The section can be expanded to try to determine which has the most trips during peak period, but for now this works.\n",
    "    \n",
    "--\n",
    "\n",
    "2. update the route pattern id field in the trips table based on maximum daily trips per route id/trip headsign/direction id combo\n",
    "    - e.g. consolidate by RHD by creating 1:1 rpid to RHD\n",
    "    \n",
    "--\n",
    "\n",
    "3. Update route/headsign/direction combos where each TOD period has less than three trips per route_pattern_id\n",
    "    - take the route_pattern_id with max number of daily trips per route/direction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FOR R/H/D, select RPID with most trips\n",
    "def max_daily_trips(tod_trips, start_stop_tod):\n",
    "    \"\"\" select the route pattern id with the most daily trips per route id, trip headsign, and direction id combination\n",
    "\n",
    "    Parameters\n",
    "    -----------\n",
    "    trips : df\n",
    "        gtfs trips.txt in df form\n",
    "    start_stop_tod : df\n",
    "        trips with tod by midpoint\n",
    "            (output of assign_tod)\n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "    max_rpid : df\n",
    "        df with route id/trip headsign/direction id and rpid with max trips\n",
    "\n",
    "    \"\"\"\n",
    "    # if rpid exists:\n",
    "    # else - rpid = rhd\n",
    "        # get the number of trips per R/H/D & rpid\n",
    "    day_rpid = tod_trips.groupby(by=['route_id','trip_headsign','direction_id','route_pattern_id']).agg({'trip_id':'nunique'})\n",
    "    day_rpid = day_rpid.rename(columns = {'trip_id':'daily_trips'}).reset_index()\n",
    "\n",
    "        # for each R/H/D, select just the route_pattern_id with the most daily trips\n",
    "    max_rpid = day_rpid.groupby(by=['route_id','trip_headsign','direction_id']).apply(lambda g: g[g['daily_trips'] == g['daily_trips'].max()])[['route_pattern_id','daily_trips']].reset_index()\n",
    "    max_rpid = max_rpid[['route_id','trip_headsign','direction_id','route_pattern_id']].rename(columns = {'route_pattern_id':'route_pattern_id_new'})\n",
    "\n",
    "    # because there are several cases where the max is held by two different route_pattern_ids, chose one arbitrarily\n",
    "    max_rpid = max_rpid.drop_duplicates(subset=['route_id','trip_headsign','direction_id'])\n",
    "    return(max_rpid, day_rpid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# update rpid to rpid with max trips\n",
    "def update_rpid_max(trips, max_trips_rpid, base_trip_columns):\n",
    "    \"\"\" update the route pattern id field in the trips table based on maximum daily trips per route id/trip headsign/direction id combo\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    trips : df\n",
    "        gtfs trips.txt in df form\n",
    "    max_trips_rpid : df\n",
    "        output df from max_daily_trips()\n",
    "        route_id/trip_headsign/direction_id -> rpid\n",
    "    base_trip_columns : idx\n",
    "        columns of trips \n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    trips_update_rpid : df\n",
    "        gtfs trips.txt in df form with updated route pattern id\n",
    "    \"\"\"\n",
    "    # join trip table with the RPID with most trips (join is on R/H/D, RPID is an attribute of the table)\n",
    "    trips_update_rpid = trips.merge(max_trips_rpid, how='left', on=['route_id','trip_headsign','direction_id'])\n",
    "\n",
    "    # update the route_pattern_ids based on R/H/D max\n",
    "        # this may be the same as original RPID as this should cover all R/H/D combos (therefore all trips)\n",
    "    trips_update_rpid['route_pattern_id'] = trips_update_rpid['route_pattern_id_new']\n",
    "    trips_update_rpid = trips_update_rpid[base_trip_columns]\n",
    "    \n",
    "    return(trips_update_rpid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_RHD_under_3(trips_update_rpid, start_stop_tod):\n",
    "    \"\"\" select route patterns with less than 3 trips in every TOD period\n",
    "\n",
    "    Parameters\n",
    "    -----------\n",
    "    trips_update_rpid : df\n",
    "        gtfs trips.txt with updated route_pattern_id (max daily trips) in df form\n",
    "    start_stop_tod : df\n",
    "        df with each trip, its assigned TOD period, and start/stop times\n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "    needs_update : df\n",
    "        route id/trip headsign/direction id with under 3 trips in each TOD period\n",
    "    update_table : df\n",
    "         route id/trip headsign/direction id with under 3 trips in each TOD period with rpid with max daily trips per route_id/t/direction_id\n",
    "    \"\"\"\n",
    "    # add tod back to trips table\n",
    "    tod_trips = trips_update_rpid.merge(start_stop_tod[['trip_id','tod']], how='left', on='trip_id')\n",
    "\n",
    "    # for each TOD and R/H/D count trips per TOD, get number of unique RPIDs, and unique ServiceIDs\n",
    "    tod_rpid = tod_trips.groupby(by=['route_id','trip_headsign','direction_id','tod']).agg({'tod':'count','route_pattern_id':'nunique','service_id':'nunique'})\n",
    "    tod_rpid = tod_rpid.rename(columns = {'tod':'trips_per_tod'}).reset_index()\n",
    "\n",
    "    # select just the trips that need to be updated\n",
    "    tod_rpid_4update = tod_rpid[['route_id','trip_headsign','direction_id','tod','trips_per_tod']]\n",
    "        # for R/H/D, get trips_per_tod separated into columns by TOD\n",
    "    tod_rpid_pivot = tod_rpid_4update.pivot_table(index = ['route_id','trip_headsign','direction_id'], columns = ['tod'])\n",
    "\n",
    "        # get just the trips where all 4 periods have less than 3 trips\n",
    "    needs_update = tod_rpid_pivot['trips_per_tod'].reset_index().fillna(0).query('AM < 3 & MD < 3 & PM < 3 & NT < 3')\n",
    "\n",
    "    # get the number of trips per R/D & rpid\n",
    "    day_rpid = tod_trips.groupby(by=['route_id','direction_id','route_pattern_id']).agg({'trip_id':'nunique'})\n",
    "    day_rpid = day_rpid.rename(columns = {'trip_id':'daily_trips'}).reset_index()\n",
    "\n",
    "        # for each R/D, select just the route_pattern_id with the most daily trips\n",
    "    max_rpid = day_rpid.groupby(by=['route_id','direction_id']).apply(lambda g: g[g['daily_trips'] == g['daily_trips'].max()])[['route_pattern_id','daily_trips']].reset_index()\n",
    "    max_rpid = max_rpid[['route_id','direction_id','route_pattern_id']].rename(columns = {'route_pattern_id':'route_pattern_id_new'})\n",
    "\n",
    "    # because there are several cases where the max is held by two different route_pattern_ids, chose one arbitrarily\n",
    "    max_rpid = max_rpid.drop_duplicates(subset=['route_id','direction_id'])\n",
    "\n",
    "    # update the route_pattern_ids for just R/H/D combos that have < 3 trips in each of the 4 TOD periods\n",
    "    update_table = needs_update.merge(max_rpid, how='left', on=['route_id','direction_id'])\n",
    "\n",
    "    update_table = update_table.drop_duplicates(subset=['route_id','direction_id','trip_headsign'])\n",
    "\n",
    "    return(needs_update, update_table)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions Part 3: Update Trip Attributes Based on Updated RPID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_attributes_new_rpid(trips_updated, update_rpid, orig_trips):\n",
    "    \"\"\" update the trip attributes based on the new route pattern id\n",
    "\n",
    "    Parameters\n",
    "    -----------\n",
    "    trips_updated : df\n",
    "        output from update_rpid_max()\n",
    "        gtfs trips.txt table with updated route pattern ids \n",
    "    \n",
    "    update_rpid : df\n",
    "        output from update_RHD_under_3()\n",
    "        route id/trip headsign/direction id with < 3 trips per TOD period, assigned new route pattern id based on max daily trips per RHD\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    trips_update5 : df\n",
    "        gtfs trip table with route pattern ids consolidated by:\n",
    "            1. route id/trip headsign/direction id (RHD)\n",
    "            2. if RHD has < 3 trips in each TOD period:\n",
    "                - route_pattern_id with max daily trips within route_id/direction_id\n",
    "            and updated attributes based on the new route pattern ids\n",
    "    \"\"\"\n",
    "    trips_updated = trips_updated.merge(\n",
    "        update_rpid[['route_id','direction_id','trip_headsign','route_pattern_id_new']], \n",
    "        how='left',on=['route_id','direction_id','trip_headsign'])\n",
    "\n",
    "    trips_updated['route_pattern_id'] = np.where(\n",
    "        ((trips_updated['route_pattern_id_new'].isna())), \n",
    "            trips_updated['route_pattern_id'], \n",
    "            trips_updated['route_pattern_id_new']\n",
    "            )\n",
    "    \n",
    "    # update headsign for all trips with updated route_pattern_id\n",
    "    rpid_th = trips_updated.groupby(by=['route_pattern_id','trip_headsign','block_id','shape_id','service_id']).agg({'trip_id':'count'}).reset_index()\n",
    "    \n",
    "    rpid_th_max = rpid_th.groupby(by=['route_pattern_id']).apply(lambda g: g[g['trip_id'] == g['trip_id'].max()]).reset_index(drop=True)\n",
    "    \n",
    "    rpid_th_max = rpid_th_max.drop_duplicates(subset='route_pattern_id')\n",
    "    \n",
    "    trips_update4 = trips_updated.merge(\n",
    "        rpid_th_max[['route_pattern_id','trip_headsign','block_id','shape_id','service_id']], \n",
    "        how='left',on=['route_pattern_id'], suffixes=(None, '_to_rplce'))\n",
    "\n",
    "    trips_update4['trip_headsign'] = trips_update4['trip_headsign_to_rplce']\n",
    "    trips_update4['block_id'] = trips_update4['block_id_to_rplce']\n",
    "    trips_update4['service_id'] = trips_update4['service_id_to_rplce']\n",
    "\n",
    "    trips_update5 = trips_update4[orig_trips.columns].merge(\n",
    "        orig_trips[['trip_id','route_pattern_id','trip_headsign']], \n",
    "        how='left', \n",
    "        on='trip_id', \n",
    "        suffixes=(None,'_old'))\n",
    "\n",
    "    return(trips_update5)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions Part 4: Stop Sequences & Stop Times\n",
    "Create Generic Stop Times & Stop Sequence per Route Pattern\n",
    "\n",
    "- Can't just change route_pattern_id as TransCAD does not use this field to combine trips into routes. There is no effect on the import.\n",
    "- Working theory is that to consolidate routes one must use the stop_times.txt table as it defines the stop sequence for every trip. Theoretically, this is being used to consolidate trips into routes based on whether they have the same stop sequence.\n",
    "\n",
    "Plan:\n",
    "- Explore if stop times differ depending on TOD or if only realtime GTFS takes into account traffic.\n",
    "    - For every route_pattern_id, get the average trip length (in minutes).\n",
    "- For every route_pattern_id, get the average number of minutes between each pair of stops in the stop sequence.\n",
    "- Replace the stop times and sequence for trips that had their route_pattern_id updated with the generic stop sequence and times per route_pattern_id created in the previous step. \n",
    "    - Keep the start time and work off of that.\n",
    "    - Arrival time will equal departure time given that the difference is usually less than a minute. Will assume difference in time can be included in the minutes to next arrival time for aggregate modeling purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generic_stop_times_sequence(stop_times, trips):\n",
    "    ''' select the most popular stop sequence per route pattern\n",
    "    ****** THIS EDITS STOP_TIMES (e.g. the input file even though not returned)\n",
    "\n",
    "    Parameters\n",
    "    -----------\n",
    "    stop_times_full : df\n",
    "        gtfs stop_times.txt in df form\n",
    "    trips : df\n",
    "        gtfs trips.txt in df form\n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "    rpid_stop_dict : dict\n",
    "        dictionary with keys = route pattern ids, values = df of trips with most popular stop sequence with that rpid\n",
    "    '''\n",
    "    # flag the first stop in every trip\n",
    "    stop_times['first_stop'] =  0\n",
    "    stop_times.loc[stop_times.groupby('trip_id').stop_sequence.idxmin(),'first_stop'] = 1\n",
    "    stop_times.loc[stop_times.groupby('trip_id').stop_sequence.idxmax(),'first_stop'] = 1000\n",
    "    \n",
    "\n",
    "    rpid_stop_dict = {}\n",
    "    # for each route_pattern_id calculate the most common stop pattern\n",
    "    for rpid in trips['route_pattern_id'].unique(): # 489\n",
    "        patterns = []\n",
    "        rpid_trips = trips.query('route_pattern_id == @rpid')['trip_id'].to_list()\n",
    "        rpid_stops = stop_times.query('@rpid_trips in trip_id').sort_values('stop_sequence')\n",
    "        \n",
    "        # for every trip, get stop pattern and save into a list of lists\n",
    "        for tid in rpid_stops['trip_id']:\n",
    "            patterns.append([tid,':'.join(rpid_stops.query('trip_id == @tid')['stop_id'].to_list())])\n",
    "        \n",
    "        # make data frame of every trip and its pattern of stops\n",
    "        df = pd.DataFrame(patterns, columns = ['trip_id','stop_pattern'])\n",
    "        # get most common stop_pattern\n",
    "        df_grby = df.groupby('stop_pattern').count().reset_index()\n",
    "        max_stop_pattern = df_grby.query('trip_id == trip_id.max()')['stop_pattern'].to_list()\n",
    "        # get all the trip_ids that represent the route_pattern_id's most popular stop pattern\n",
    "        rep_trips = df.query('stop_pattern == @max_stop_pattern[0]')['trip_id'].to_list()\n",
    "\n",
    "        rpid_stop_dict[rpid] = rep_trips\n",
    "\n",
    "    return(rpid_stop_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_to_majority_stop_pattern(trips, rpid_stop_dict, stop_times):\n",
    "    '''for trips without the most popular stop sequence pattern within its route pattern, update to the most popular stop sequence pattern \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    trips : df\n",
    "        gtfs trips.txt in df form\n",
    "    rpid_stop_dict : dict\n",
    "        output of generic_stop_times_sequence()\n",
    "        dictionary with keys = route pattern ids, values = df of trips with most popular stop sequence with that rpid\n",
    "    stop_times : df\n",
    "        gtfs stop_times.txt in df form\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    trips_clean : df\n",
    "        gtfs trips.txt in df form with updated shape_id based on new stop sequence pattern (majority within rpid)\n",
    "    trip_stop_replace : dict\n",
    "        dictionary, key = trip_id, values = trip_id of trip with closest start time to key within the stop sequence majority for the shared rpid\n",
    "\n",
    "    '''\n",
    "    # for every row where the route_pattern_id has been updated & doesn't have the majority stop pattern\n",
    "    trip_stop_replace = {}\n",
    "    trip_shape_replace = {}\n",
    "    for idx, row in trips.iterrows():\n",
    "        # identify trip_id\n",
    "        tid = row['trip_id']\n",
    "        # get the route_pattern_id for the trip\n",
    "        rpid = row['route_pattern_id']\n",
    "\n",
    "        if (tid not in rpid_stop_dict[rpid]): # if trip_id is not in the majority stop pattern\n",
    "            # get all of the trips associated with that route_pattern_id (and stop sequence pattern)\n",
    "            all_trips = rpid_stop_dict[rpid]\n",
    "\n",
    "            # get the start time (for first stop) for trip_id\n",
    "            start_time = stop_times.query('(trip_id == @tid) & (first_stop == 1)')\n",
    "            # get the start time (for first stop) for all trips that share same route_pattern_id\n",
    "            all_start_times = stop_times.query('(trip_id in @all_trips) & (first_stop == 1)')\n",
    "\n",
    "            # create list getting the difference in start times between the selected trip and all the other trips that share the same route_pattern_id\n",
    "            test_list = [[x,(abs(start_time['departure_time_sec']- all_start_times.query('trip_id == @x'))['departure_time_sec'].iloc[0])] for x in all_start_times['trip_id']]\n",
    "            close = {}\n",
    "            close = {sub[0]:sub[1] for sub in test_list}\n",
    "                \n",
    "            # get trip with the minimum difference in start time    \n",
    "            min_t = min(close, key=close.get)\n",
    "            if tid != min_t:\n",
    "                trip_stop_replace[tid] = min_t\n",
    "                trip_shape_replace[tid] = trips.query('trip_id == @min_t')['shape_id'].values[0]\n",
    "        else:\n",
    "            trip_shape_replace[tid] = trips.query('trip_id == @tid')['shape_id'].values[0]\n",
    "    \n",
    "    trip_shape_replace_tab = pd.DataFrame.from_dict(\n",
    "        trip_shape_replace, orient='index').reset_index().rename(\n",
    "            columns = {'index':'trip_id', 0:'shape_id_update'})\n",
    "    \n",
    "    trips['shape_id_new'] = trips['trip_id'].apply(lambda x: trip_shape_replace[x])\n",
    "\n",
    "    trips['shape_id'] = trips['shape_id_new']\n",
    "\n",
    "    # keep only relevant columns\n",
    "    trips_clean = trips.drop(columns=['route_pattern_id_old','trip_headsign_old','shape_id_new'])\n",
    "\n",
    "    return(trips_clean, trip_stop_replace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_stop_times_new_pattern(stop_times, trip_stop_replace):\n",
    "    ''' for trips with stop sequence pattern, update the arrival and departure times associated with each stop\n",
    "\n",
    "    Parameters\n",
    "    -----------\n",
    "    stop_times : df\n",
    "        gtfs stop_times.txt in df form\n",
    "    trip_stop_replace : dict\n",
    "        dictionary, key = trip_id, values = trip_id of trip with closest start time to key within the stop sequence majority for the shared rpid\n",
    "    stop_times_columns : idx\n",
    "        original columns of gtfs stop_times.txt\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    stop_times : df\n",
    "        gtfs stop_times.txt in df form with updated stops times per trip \n",
    "\n",
    "    '''\n",
    "    stop_times = stop_times.sort_values('stop_sequence')\n",
    "    # for every trip (in dictionary with value = trip sharing route_pattern_id with closest start time)\n",
    "    for tid in trip_stop_replace.values():\n",
    "    # update value of time_between_stops for every stop in the selected trip\n",
    "        # time between stops calculated by subtracting the prior departure_time_sec from current (why stop_sequence order is important)\n",
    "        stop_times.loc[stop_times.loc[:,'trip_id'] == str(tid), 'time_between_stops'] = stop_times.loc[stop_times.loc[:,'trip_id'] == str(tid), 'departure_time_sec'].diff()\n",
    "\n",
    "    for trip in trip_stop_replace.keys():\n",
    "        start_time = stop_times.query('(trip_id == @trip) & (first_stop == 1)')['departure_time_sec']\n",
    "        end_time = stop_times.query('trip_id == @trip').loc[stop_times.query('trip_id == @trip')['stop_sequence'].idxmax()]['departure_time_sec']\n",
    "        # drop old stop times\n",
    "        stop_times = stop_times.drop(\n",
    "            stop_times.loc[stop_times['trip_id']==trip].index)\n",
    "        # grab new stop times\n",
    "        new_trip = trip_stop_replace[trip]\n",
    "        nst = stop_times.query('trip_id == @new_trip')\n",
    "        nst['trip_id'] = trip\n",
    "\n",
    "        # replace the start time, then calculate the stop times by the departure_time_sec difference\n",
    "        nst.loc[nst.loc[:,'first_stop']==1,'departure_time_sec'] = int(start_time.iloc[0])\n",
    "        nst.loc[nst.loc[:,'first_stop']==1,'time_between_stops'] = int(start_time.iloc[0])\n",
    "        nst['departure_time_sec'] = nst['time_between_stops'].cumsum()\n",
    "\n",
    "        # recalc arrival/dep times\n",
    "        nst['arrival_time'] = pd.to_datetime(nst['departure_time_sec'],unit='s').astype('str').str[11:19]\n",
    "        # if intermediate times are null, replace last one with end time\n",
    "        nst['arrival_time']=np.where((nst['first_stop'] == 1000) & ((nst['arrival_time'] == \"\") | (nst['arrival_time'].isnull())), end_time, nst['arrival_time'])\n",
    "        #b = nst[(nst['arrival_time'] == \"\") | (nst['arrival_time'].isnull())]['stop_sequence'].max()\n",
    "        #nst['arrival_time']=np.where((nst['stop_sequence'] == b) & ((nst['arrival_time'] == \"\") | (nst['arrival_time'].isnull())), end_time, nst['arrival_time'])\n",
    "        \n",
    "        nst['departure_time'] = nst['arrival_time']\n",
    "\n",
    "        #keep only relevant columns\n",
    "        nst = nst[stop_times.columns]\n",
    "\n",
    "        stop_times = pd.concat([stop_times,nst])\n",
    "    # keep only relevant columns and sort\n",
    "    stop_times = stop_times.sort_values(by=['trip_id','stop_sequence'])\n",
    "\n",
    "    return(stop_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dumb_routes(stop_times):\n",
    "    stop_times.loc[(stop_times['first_stop']==1000) & (stop_times['departure_time_sec'].isnull()),'arrival_time'] = pd.to_timedelta(stop_times[(stop_times['first_stop']==1000) & (stop_times['departure_time_sec'].isnull())]['arrival_time'],unit='s').astype('str').str[8:].str.zfill(8)\n",
    "    stop_times['departure_time'] = stop_times['arrival_time']\n",
    "    return(stop_times)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Functions!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Run Functions to...\n",
    "- separate bus and non-bus trip\n",
    "- get start and stop times for trips\n",
    "- assign each trip a tod based on midpoint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All routes in this GTFS are bus routes.\n"
     ]
    }
   ],
   "source": [
    "non_bus_trips, bus_trips = separate_bus(gtfsfeeds_dfs.routes, gtfsfeeds_dfs.trips,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_stop = get_start_stop_times(gtfsfeeds_dfs.stop_times) # simpson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_stop_tod = assign_tod(start_stop) # smurf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Please run the functions that follows - you likely need to consolidate your route patterns\n"
     ]
    }
   ],
   "source": [
    "trips, count_tod = assign_trips_tod(bus_trips, start_stop_tod)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Run Functions to...\n",
    "- get max daily trips per RHD/rpid\n",
    "- update trips rpid with rpid with max daily trips (per RHD)\n",
    "- if RHD has < 3 trips in each TOD period, update trips rpid with rpid with max daily trips (per Route_ID/Direction_ID) \n",
    "- update trip attributes based on updated rpid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_trips_rpid, day_trips = max_daily_trips(trips, start_stop_tod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "trips_update1 = update_rpid_max(bus_trips, max_trips_rpid, gtfsfeeds_dfs.trips.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "trips_need_update, updated_rpid = update_RHD_under_3(trips_update1, start_stop_tod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "trips_fully_updated = update_attributes_new_rpid(trips_update1, updated_rpid, bus_trips)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. Run Functions to...\n",
    "- select the most popular stop sequence per RPID\n",
    "- replace stop sequence for trips with minority stop sequence (with most popular sequence per RPID)\n",
    "    - also update shape in trips table based on new stop sequence\n",
    "- update the stop times associated with the new stop sequence patterns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "rpid_stop_dict = generic_stop_times_sequence(gtfsfeeds_dfs.stop_times, trips_fully_updated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "trips, trip_stop_replace = update_to_majority_stop_pattern(trips_fully_updated, rpid_stop_dict, gtfsfeeds_dfs.stop_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "stop_times = update_stop_times_new_pattern(gtfsfeeds_dfs.stop_times, trip_stop_replace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "stoptimes = dumb_routes(stop_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "gtfsfeeds_dfs.stop_times = stoptimes.drop(columns=['first_stop'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Clean Up & Export"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge back the filtered out trips\n",
    "gtfsfeeds_dfs.trips = pd.concat([trips, non_bus_trips])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "gtfsfeeds_dfs.stop_times.to_csv(out_path+r\"/stop_times.txt\", index=False)\n",
    "gtfsfeeds_dfs.trips.to_csv(out_path+r\"/trips.txt\", index=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "urbanAccess",
   "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.8.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "2a2b370b38977169157650f5355ad729af7f449719cdba662b087dc855e43e33"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}