<jupyter_start><jupyter_text>Once more, I'm trying to show a play-by-play analysis of the Israeli GTFS, with [hasadna/open-bus](https://github.com/hasadna/open-bus) interests in mind.This post is more of a data cleanup diary log. While working on some stuff I found issues with the dataset I was working with, I thought it would be good to document this, and show the process of finding these issues, and how some simple analysis (AKA debugging) helped finding the source of these symptoms and a cure.The agenda:1. [Get up-to-date GTFS files from the MOT FTP server](Get-the-data)2. [Merge and tidy](Tidy-it-up)3. [Cleanup Diary](Cleanup-Diary)<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import datetime

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Get the dataWe get the data straight from MOT's ftp. Just so I don't repeat the last post, this went out to a separate module.<jupyter_code>from gtfs_utils import *
LOCAL_ZIP_PATH = 'data/sample/gtfs_171218.zip' 

conn = ftp_connect()
ftp_dir = get_ftp_dir(conn)
UPTODATE = 90 #days
our_uptodateness = get_uptodateness(ftp_dir, local_zip_path = LOCAL_ZIP_PATH)

if our_uptodateness > UPTODATE:
    get_ftp_file(conn)
    get_ftp_file(conn, file_name = 'Tariff.zip', local_zip_path = 'data/sample/tariff.zip' )

conn.quit()<jupyter_output><empty_output><jupyter_text>We load MOT's tarrif.txt file, which will give us zone names and other local info.Found out it's the file is even more fd-up than I realized, so I add the southern zones manualy<jupyter_code>tariff_df = extract_tariff_df(local_zip_path = 'data/sample/tariff.zip')
south = [
    {
        'zone_name': 'מצפה רמון',
        'zone_id': '903'
    },
    {
        'zone_name': 'ערבה',
        'zone_id': '902'
    },
    {
        'zone_name': 'אילת',
        'zone_id': '901'
    },]
south = pd.DataFrame(south)
tariff_df = tariff_df.append(south)<jupyter_output><empty_output><jupyter_text>Tidy it upAgain I'm using [partridge](https://github.com/remix/partridge/tree/master/partridge) for filtering on dates, and then some tidying up and transformations.<jupyter_code>import partridge as ptg

service_ids_by_date = ptg.read_service_ids_by_date(LOCAL_ZIP_PATH)
service_ids = service_ids_by_date[datetime.date(2017, 12, 21)]

feed = ptg.feed(LOCAL_ZIP_PATH, view={
    'trips.txt': {
        'service_id': service_ids,
    },
})
def to_timedelta(df):
    '''
    Turn time columns into timedelta dtype
    '''
    cols = ['arrival_time', 'departure_time']
    numeric = df[cols].apply(pd.to_timedelta, unit='s')
    df = df.copy()
    df[cols] = numeric
    return df
s = feed.stops
r = feed.routes
t = (feed.trips
     .assign(route_id=lambda x: pd.Categorical(x['route_id'])))
f = (feed.stop_times[['trip_id', 'departure_time', 'arrival_time', 'stop_id', 'stop_sequence']]
     .assign(date = datetime.date(2017, 12, 21))
     .merge(s[['stop_id', 'stop_name', 'stop_lat', 'stop_lon', 'zone_id']], on='stop_id')
     # Much faster joins and slices with Categorical dtypes
     .merge(tariff_df[['zone_id', 'zone_name']], on='zone_id')
     .assign(zone_id=lambda x: pd.Categorical(x['zone_id']))
     .assign(zone_name=lambda x: pd.Categorical(x['zone_name']))
     .merge(t[['trip_id', 'route_id', 'direction_id']], on='trip_id')
     .merge(r[['route_id', 'route_short_name', 'route_long_name']], on='route_id')
     .assign(route_id=lambda x: pd.Categorical(x['route_id']))
     .pipe(to_timedelta)
    )
f.head()<jupyter_output><empty_output><jupyter_text>Cleanup DiaryThe question that ignited this post was about final destinations. This is an important parameter we have to be familiar with - last stop of every trip in `trips.txt`, based on `stop_times.txt` - since our current best polling strategy sends requests for these stops only.  One (wrong) wayI first tried to use a dataframe I constructed when I was trying to find simple loop lines (A --> A). I (naively) thought to myself, lets get the first and last stop for all trips.I did it like so -<jupyter_code>import pickle
import os
# using a custom groupby aggregation function, this is VERY inefficient,
# but far more readable than doing some other magic, so I'm sticking to it for now

def apply_and_pickle(f, pkl_path, apply_func):
    if not os.path.exists(pkl_path):
        ret_df = f.groupby('trip_id').apply(apply_func)
        # since it's so inefficient I pickle it
        with open(pkl_path, 'wb') as pkl_file:
            pickle.dump(orig_dest_only, pkl_file)
    else:
        with open(pkl_path, 'rb') as pkl_file:
            ret_df = pickle.load(pkl_file)
        return ret_df

def first_last(df):
    return df.sort_values(by='stop_sequence').iloc[[0, -1]]

orig_dest_only = apply_and_pickle(f, 'softloop.pkl', first_last)
fig, ax = plt.subplots(figsize=(10,7))
ax.set_yticks(np.arange(32)*50)
orig_dest_only[orig_dest_only.stop_sequence!=1].groupby('stop_id').size().sort_values(ascending=False).hist(bins=147, ax=ax)<jupyter_output><empty_output><jupyter_text>I was all jolly and happy, until I did this:<jupyter_code>orig_dest_only[orig_dest_only.stop_sequence!=1].groupby('stop_id').size().sum()
t.shape<jupyter_output><empty_output><jupyter_text>No good. I've actually found more last stops than trips...This could be due to a gazillion reasons, so lets go ahead and go over them until we nail it. stop_sequence Mystery or "The Bus with No Origin"I used stop_sequence to find first and last stops, maybe there's an issue in the dataLet's just do a simple count -<jupyter_code>f[f.stop_sequence==1].shape
f[f.stop_sequence==2].shape<jupyter_output><empty_output><jupyter_text>There's less first stop than trips, but more 2nd stops than trips. Ok this is crazy, and certainly needs to be looked at in the future, but lets just leave it and use some other field for this task.TODO: dig into `stop_sequence` and talk to MOT about this  Better yet, or "Timing is Everything"Lets use the `arrival_time` field instead, which should be cleaner. I'll skip the inefficient method I used before with the custom groupby aggregationm, and go for some neat trick using the mighty `transform` method.<jupyter_code>idx = f.groupby(['trip_id'], sort=False)['arrival_time'].transform(max) == f['arrival_time']
f[idx].sort_values(by='trip_id').head()<jupyter_output><empty_output><jupyter_text>That took less than a second to run, now lets see if we also got better results<jupyter_code>f[idx].shape<jupyter_output><empty_output><jupyter_text>Ok we're on the good side of the error this time (I think), with some false positives instead of negatives.But now it means we might have duplicate trips in our dataset :(Let's zoom-in on one of the duplicate trips<jupyter_code>f[idx].trip_id.value_counts().head()
st = feed.stop_times
st[((st.trip_id=='24167955_171217') & (st.stop_sequence==31)) ]
f[((f.trip_id=='24167955_171217') & (f.stop_sequence==31)) ]<jupyter_output><empty_output><jupyter_text>Ok, that looks better, this problem is only in my `f` feed merged data frame, not in the original `stop_times.txt`.I debug `f`'s creation to see (by building it one step at a time), I spare you these steps, and give you the problematic step:<jupyter_code>ftest = (feed.stop_times[['trip_id', 'departure_time', 'arrival_time', 'stop_id', 'stop_sequence']]
     .assign(date = datetime.date(2017, 12, 21))
     .merge(s[['stop_id', 'stop_name', 'stop_lat', 'stop_lon', 'zone_id']], on='stop_id')
      # Much faster joins and slices with Categorical dtypes
     .merge(tariff_df[['zone_id', 'zone_name']], on='zone_id')

    )
ftest[ ((ftest.trip_id=='10230038_171217') & (ftest.stop_sequence==31)) ]
tariff_df[tariff_df.zone_id=='232']<jupyter_output><empty_output><jupyter_text>Got it! Source of all evil - the super-messy tariff file where we get the zone names from.<jupyter_code>tariff_df.groupby(['zone_id', 'zone_name']).size()<jupyter_output><empty_output><jupyter_text>See if you can find the needle...Now lets rebuild our `f` dataset with deduped zones:<jupyter_code>ftry = (feed.stop_times[['trip_id', 'departure_time', 'arrival_time', 'stop_id', 'stop_sequence']]
     .assign(date = datetime.date(2017, 12, 21))
     .merge(s[['stop_id', 'stop_name', 'stop_lat', 'stop_lon', 'zone_id']], on='stop_id')
      # Much faster joins and slices with Categorical dtypes
     .merge(tariff_df.groupby(['zone_id', 'zone_name']).size().reset_index()[['zone_id', 'zone_name']], on='zone_id')
    )
ftry [ ((ftry.trip_id=='10230038_171217') & (ftry.stop_sequence==31)) ]
new_f = (feed.stop_times[['trip_id', 'departure_time', 'arrival_time', 'stop_id', 'stop_sequence']]
     .assign(date = datetime.date(2017, 12, 21))
     .merge(s[['stop_id', 'stop_name', 'stop_lat', 'stop_lon', 'zone_id']], on='stop_id')
     # Much faster joins and slices with Categorical dtypes
     .merge(tariff_df.groupby(['zone_id', 'zone_name']).size().reset_index()[['zone_id', 'zone_name']], on='zone_id')
     .assign(zone_id=lambda x: pd.Categorical(x['zone_id']))
     .assign(zone_name=lambda x: pd.Categorical(x['zone_name']))
     .merge(t[['trip_id', 'route_id', 'direction_id']], on='trip_id')
     .merge(r[['route_id', 'route_short_name', 'route_long_name']], on='route_id')
     .assign(route_id=lambda x: pd.Categorical(x['route_id']))
     .pipe(to_timedelta)
    )
new_f.head()<jupyter_output><empty_output><jupyter_text>TODO: apply this fix to previous and in-progress notebooks<jupyter_code>idx = new_f.groupby(['trip_id'], sort=False)['arrival_time'].transform(max) == new_f['arrival_time']
new_f[idx].trip_id.value_counts().head()
new_f[idx].shape<jupyter_output><empty_output><jupyter_text>Trip duplicates - **check.**And now we're actually missing 9. hmmmm... Lets see which ones we're missing.<jupyter_code>t = t.set_index('trip_id')
missing_trips = t[~t.index.isin(new_f.trip_id)]
missing_trips<jupyter_output><empty_output><jupyter_text>Maybe they just don't have stop times?<jupyter_code>st[st.trip_id.isin(missing_trips.index)].shape<jupyter_output><empty_output><jupyter_text>Nope. Good!To be safe lets test we are getting `stop_times` rows from the same line of code, on another series of trips<jupyter_code>test_trips = t.head()
test_trips
st[st.trip_id.isin(test_trips.index)].shape<jupyter_output><empty_output><jupyter_text>Missing trip - **check.**Woohoo! looks like we're good to go (For now :O ) GoodbyeTo shortly sum up - cleanup is an ongoing process. I have found an issue in the original data (stop_sequence) and an issue with how I merged tables (double zone in the very very messy tariffs table). See you soon when I test [GTFS best practices](http://gtfs.org/best-practices) on our GTFS files (with regards to loops routes, lasso routes and some more).<jupyter_code>idx = new_f.groupby(['trip_id'], sort=False)['arrival_time'].transform(min) == new_f['arrival_time']
new_f[idx].trip_id.value_counts().head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13091743119266055, 'fraction_numerical': 0.013302752293577982, 'mean_word_length': 5.3823185011709604, 'pattern_counts': {'":': 0, '<': 73, '<?xml version=': 0, '>': 75, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '414567', 'license': 'MIT', 'n_tokens_mistral': 4014, 'n_tokens_neox': 3789, 'n_words': 1010, 'path': 'openbus_4_cleanup.ipynb', 'repo_name': 'cjer/open-bus-explore', 'size': '10900'}