<jupyter_start><jupyter_text>Exploring Israel's MOT GTFS files...Over at [hasadna/open-bus](https://github.com/hasadna/open-bus) we want to build an archive of country-wide bus arrival times. Our way to get there is through a [SIRI interface](https://www.gov.il/he/Departments/General/real_time_information_siri) provided by the Ministry of Transportation. We need to continually make SOAP requests, in such a way so that their aggregated bus location responses will allow us to calculate arrival times for all bus rides across the country. In order to do so, I want to first analyze MOT's [static GTFS data](https://www.gov.il/he/Departments/General/gtfs_general_transit_feed_specifications) to figure out what kind of numbers we're talking about, so we can start talking about design and scale for building an efficient and robust archiving module for SIRI realtime info. The agenda for this notebook is roughly so: (you can skip to item 3 if you're just interested in the stats):1. [Get up-to-date GTFS files from the MOT FTP server](Get-the-data)2. [Extract the needed info from them into tidy pandas DataFrames](Tidy-it-up)3. [Some basic exploratory analysis examples](Basic-exploratory-analysis)4. [One the next episodes](On-the-next-episodes) Yep, not all is covered in this notebook*side note - Hasadna's OpenBus project has the data in a DB. I'm not using it here, since want this notebook to be standalone for further exploration*  Get the dataStart by making an FTP connection and retrieve the gtfs zip file's modification date. We're restricted from using most commands, but we can LIST...<jupyter_code>MOT_FTP = 'gtfs.mot.gov.il'
FILE_NAME = 'israel-public-transportation.zip'
UPTODATE = 150 #days
local_zip_path = 'data/sample/gtfs.zip' 

conn = FTP(MOT_FTP)
conn.login()
ftp_dir = []
conn.retrlines('LIST', lambda x: ftp_dir.append(x)) 
ftp_dir<jupyter_output><empty_output><jupyter_text>Now that we have a list of file records, lets extract the dates and see if our file is up-to-date<jupyter_code>RE_FTP_LINE = re.compile(
    r'(?P<date>\d+-\d+-\d+\s+\d+:\d+[APM]{2})\s+(?P<size><DIR>|[0-9]+)\s+(?P<file_name>.*)')
f = [re.findall(RE_FTP_LINE, line) for line in ftp_dir]
f_dates = {t[0][2]: datetime.datetime.strptime(t[0][0], "%m-%d-%y  %H:%M%p") for t in f}

ftp_date = f_dates[FILE_NAME]

have_zip = os.path.exists(local_zip_path)
if have_zip:
    our_date = datetime.datetime.fromtimestamp(os.path.getmtime(local_zip_path))

    our_uptodateness = (ftp_date - our_date).days
    print(f'Our zip is {our_uptodateness} days old')
    if our_uptodateness > UPTODATE:
        print('Oh snap we need to download a new zip...')
    else:
        print('Oh goody! we\'re up-to-date!')<jupyter_output>Our zip is 37 days old
Oh goody! we're up-to-date!<jupyter_text>And now we shall actually download the zip... if we really want to...<jupyter_code>if (have_zip and our_uptodateness > UPTODATE) or not have_zip:
    temp_local_path = 'data/sample/gtfs_temp.zip'
    fh = open(temp_local_path, 'wb')
    conn.retrbinary('RETR %s' % (FILE_NAME), fh.write)
    fh.close()
    #TODO THIS DOESN'T WORK IF YOU DIDN"T HAVE FILE
    os.remove(local_zip_path)
    os.rename(temp_local_path, local_zip_path)

conn.quit()
    
h5_file_path = 'data/gtfs_store.h5'
if our_uptodateness > UPTODATE:
    os.remove(h5_file_path)
store = pd.HDFStore(h5_file_path)
if our_uptodateness > UPTODATE:
    with zipfile.ZipFile(local_zip_path) as zf:
        for fn in zf.namelist():
            short = fn.split('.')[0]
            if short not in store:
                store[short] = pd.read_csv(zf.open(fn))

store.stops.head()<jupyter_output><empty_output><jupyter_text>Lets also load MOT's tarrif.txt file, which will give us zone names and other local info<jupyter_code>TARIFF_FILE_NAME = 'Tariff.zip'
TARIFF_TXT_NAME = 'Tariff.txt'
TARIFF_TO_REFORM_ZONE = 'StationToReformZone.txt'
local_tariff_path = 'data/sample/tariff.zip' 

conn = FTP(MOT_FTP)
conn.login()
fh = open(local_tariff_path, 'wb')
conn.retrbinary('RETR %s' % (TARIFF_FILE_NAME), fh.write)
fh.close()
conn.quit()<jupyter_output><empty_output><jupyter_text>Tidy it upMany of Seaborn's built-in features call for the data to be in a Pandas DataFrame and to take the form of what Hadley Whickam calls [“tidy” data](http://vita.had.co.nz/papers/tidy-data.pdf). In brief, that means your dataframe should be structured such that each column is a variable and each row is an observation.Moreover, as <NAME> [claims here](https://tomaugspurger.github.io/modern-5-tidy.html), starting an analysis with a tidy dataset should truly help *"facilitate analysis: going as quickly as possible from question to answer."* I'm taking his word on this now, and hope to show this in work here too.Later, We'll also play a little bit with pivoting between long and wide formats.We'll start with one tidy dataset for stop_times, I'll use [partridge](https://github.com/remix/partridge/tree/master/partridge) for filtering on dates.<jupyter_code>import partridge as ptg

service_ids_by_date = ptg.read_service_ids_by_date(local_zip_path)
next_thursday = None #TODO: add this instead of hardcoding the date
service_ids = service_ids_by_date[datetime.date(2018, 3, 1)]

feed = ptg.feed(local_zip_path, view={
    'trips.txt': {
        'service_id': service_ids,
    },
})<jupyter_output><empty_output><jupyter_text>The partridge feed now references all the GTFS files in the zip, and will only load them to memory only once used<jupyter_code>[x for x in dir(feed) if not x.startswith('__')]<jupyter_output><empty_output><jupyter_text>Getting a DataFrame is as easy as...<jupyter_code>s = feed.stops
s.head()<jupyter_output><empty_output><jupyter_text>Before we build our first stop times tidy DataFrame, we need zone name data. We need to combine `Tariff.txt` with `StationToReformZone.txt` to a mapping as described in MOT's developers [pdf](https://www.gov.il/BlobFolder/generalpage/gtfs_general_transit_feed_specifications/he/%D7%A1%D7%98%20%D7%A7%D7%91%D7%A6%D7%99%D7%9D%20-%20GTFS%20-%20%D7%94%D7%A1%D7%91%D7%A8%20%D7%9C%D7%9E%D7%A4%D7%AA%D7%97%D7%99%D7%9D_0.pdf).**Note**: I'm disregarding `FromDate` and `ToDate` and just dropping duplicates. So we will give wrong zone names for zone changes within the GTFS 60 day period.<jupyter_code># not a true csv, so we need to jiggle it a bit
cols = ['ShareCode','ShareCodeDesc','ZoneCodes','Daily','Weekly','Monthly','FromDate','ToDate', 'EXTRA']
reform_cols = ['StationId', 'ReformZoneCode','FromDate','ToDate', 'EXTRA']
with zipfile.ZipFile(local_tariff_path) as zf:
    tariff_df = (pd.read_csv(zf.open(TARIFF_TXT_NAME), header=None, skiprows=[0], names = cols)
                .drop(columns = ['EXTRA']))
    reform_df = (pd.read_csv(zf.open(TARIFF_TO_REFORM_ZONE), header=None, skiprows=[0], names = reform_cols)
                 .drop(columns = ['EXTRA']))

    
# remove ShareCodes which contain multiple zones  e.g. גוש דן מורחב
tariff_df = (tariff_df[~ tariff_df.ZoneCodes.str.contains(';')]
             .rename(columns = {'ShareCodeDesc': 'zone_name',
                               'ZoneCodes': 'zone_id'}))
tariff_df.head()
rs = reform_df[['StationId', 'ReformZoneCode']].drop_duplicates().applymap(str).set_index('StationId').iloc[:,0]
rs.head()
ts = (tariff_df[['zone_id', 'zone_name']].drop_duplicates().set_index('zone_id').iloc[:,0])
ts.head()
zones = rs.map(ts).reset_index().rename(columns={'StationId': 'stop_code', 'ReformZoneCode':'zone_name'})<jupyter_output><empty_output><jupyter_text>And we're ready to build our final stop_times DataFrame. Not a lot of work on it, as it is already in a good tidy format, just select the columns we want, do some transformations and merges.<jupyter_code>def to_timedelta(df):
    '''
    Turn time columns into timedelta dtype
    '''
    cols = ['arrival_time']
    numeric = df[cols].apply(pd.to_timedelta, unit='s')
    df = df.copy()
    df[cols] = numeric
    return df

f = (feed.stop_times[['trip_id', 'arrival_time', 'stop_id']]
     .assign(date = datetime.date(2017, 12, 21))
     .merge(s[['stop_id', 'stop_name', 'stop_lat', 'stop_lon', 'stop_code']], on='stop_id')
     # Much faster joins and slices with Categorical dtypes
     .merge(zones, how='left')
     .assign(zone_name=lambda x: pd.Categorical(x['zone_name']))
     .pipe(to_timedelta)
    )
f.head()<jupyter_output><empty_output><jupyter_text>Basic exploratory analysis  Big numbers Overall daily stops<jupyter_code>f.shape[0]<jupyter_output><empty_output><jupyter_text>That's a big number, from here on we'll try to better understand its distribution  Overall stops per minuteFirst time that I play with this kind of timeseries with pandas, so I'll try a few different ways (including the deprecated seaborn tsplot). Not gonna try to explain how these methods compare... I'm keeping the diversity mostly for educational purposes.<jupyter_code>f_1T = f.set_index(['arrival_time']).resample('1T')
f_1T.size().plot()
sns.despine()
f_1T.size().reset_index().head()
from statsmodels.nonparametric.smoothers_lowess import lowess
counts = f_1T.size().reset_index().rename(columns={0: 'size'})
smooth_counts = lowess(counts.iloc[:,1], counts.arrival_time, frac=0.01) 
counts['smooth'] = smooth_counts[:, 1]
counts = counts.set_index('arrival_time')
counts['ewm'] = counts['size'].ewm(span=10).mean()
counts.head()
counts[['size', 'smooth']].plot()
sns.despine()
counts[['size', 'ewm']].plot()
sns.despine()<jupyter_output><empty_output><jupyter_text>Now build a DF to use with Seaborn tsplot (even though its deprecated)<jupyter_code># one minute counts
c10 = f_1T.size().reset_index().rename(columns={0: 'size'})
# extract minute smalll digit to a new column
c10 = c10.assign(minute=c10.arrival_time.dt.seconds//60%10)
# group the one minute counts into a two level MultiIndex by the 1 minute counts
# and a new resample by 10 minutes. the minute
c10 = c10.set_index('arrival_time').groupby([pd.Grouper(freq='10T'), 'minute']).mean()
# lets look at it in a meaningful way, we can see a peak at the 0  and 5
# (sparkline could be a good viz for this)
c10.reset_index().pivot('arrival_time', 'minute', 'size').iloc[6*7:6*9]
sns.tsplot(data=c10.reset_index(), time='arrival_time', unit='minute', value='size',
           ci='sd')
sns.despine()<jupyter_output>C:\Users\dbare\Anaconda3\lib\site-packages\seaborn\timeseries.py:183: UserWarning: The tsplot function is deprecated and will be removed or replaced (in a substantially altered version) in a future release.
  warnings.warn(msg, UserWarning)<jupyter_text>The lower opacity color represents the standard deviation.That's kind of simple to do on our own, so lets do this and some more.<jupyter_code>f10 = f_1T.size().resample('10T').apply(['max', 'min', 'mean', 'std'])
f10['ci_l'] = f10['mean'] - f10['std']
f10['ci_h'] = f10['mean'] + f10['std']
f10.head()
f10.plot()
sns.despine()<jupyter_output><empty_output><jupyter_text>Again, one thing to note is the noticeable zigzags of the lines. This is due to round 10 minutes or hours have much more buses scheduled.  Photo finishLet's finish with some zone slicing. Starting by getting the top zones overallKind of a teaser for future posts, where I'll dive deeper with more slices, groupings and some dynamic dashboarding (woohoo).Start with getting the top ten zones:<jupyter_code>top_zones = f.groupby('zone_name').size().sort_values(ascending=False).head(10)
top_zones.reset_index()
zone10 = (f[f.zone_name.isin(top_zones.index)].set_index('arrival_time')
          .groupby([pd.Grouper(freq='1T'), 'zone_name'])
          .size())
zones10 = zone10.reset_index().rename(columns={0: 'size'})
#reverse names so they will show nicely in plot
zones10 = zones10.assign(zone_name=lambda x: (x['zone_name'].str[::-1]))
zones10 = zones10.pivot('arrival_time', 'zone_name', 'size')
fig, ax = plt.subplots()
zones10.plot(figsize=(15,10), ax=ax)
handles, labels = ax.get_legend_handles_labels()
# sort both labels and handles by their rank (and do so in a very weird way :))
top_zones_rank = top_zones.reset_index()
labels, handles = zip(*sorted(zip(labels, handles), 
                              key=lambda t: top_zones_rank[top_zones_rank.zone_name==t[0][::-1]].index.tolist()[0]))
ax.legend(handles, labels, title="zone_name")

sns.despine()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12947591145833334, 'fraction_numerical': 0.016276041666666668, 'mean_word_length': 5.141429285357321, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 64, 'https://': 6, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '288048', 'license': 'MIT', 'n_tokens_mistral': 4474, 'n_tokens_neox': 4259, 'n_words': 1177, 'path': 'openbus_1_basics.ipynb', 'repo_name': 'cjer/open-bus-explore', 'size': '12296'}