<jupyter_start><jupyter_text>MBTA Vehicle UpdatesI've been archiving the output of the MBTA's realtime vehicle updates API and making it available for download from http://mbta-history.apptic.xyz. This notebook illustrates how to retrieve and use the data programmatically with Python 3 and Pandas.Ultimately, we'll create a visualization of actual arrival time vs. scheduled arrival time that looks something like this:![Red Line delays](https://www.dropbox.com/s/px7u6oriz7dqv1j/Unknown.png?dl=1)<jupyter_code>from datetime import datetime, timedelta
from io import BytesIO, TextIOWrapper
import os
import tarfile
from urllib.error import HTTPError
from urllib.request import urlopen, urlretrieve

import numpy as np
import pandas as pd
import pytz
try:
  from ipywidgets import IntProgress, HBox, HTML
  from IPython.display import display

  def retrieve(*args, reporthook=None, **kwargs):
      "Wrapper around urlretrieve that displays a progress bar."
      bar = IntProgress(min=0, max=1, value=0)
      display(bar)
      def do_report(block, read, total):
          if total > bar.min:
              bar.max = total
              bar.value = (block+1)*read
          if reporthook:
              reporthook(block, read, total)

      return urlretrieve(*args, reporthook=do_report, **kwargs)
except ImportError:
  retrieve = urlretrieve<jupyter_output><empty_output><jupyter_text>We'll define a function `get_updates` that takes a `datetime` argument and returns a pandas DataFrame containing all reported vehicle updates for that day (assuming they're available). You can also visit http://mbta-history.apptic.xyz directly, download the data for a particular day, then use `pandas.read_csv(path_to_csv_file)` to construct a dataframe.<jupyter_code>def archive_url(dt):
    csv_filename = dt.strftime("%Y-%m-%d")
    return dt.strftime(f"http://mbta-history.apptic.xyz/%Y/%m/{csv_filename}.csv.gz")

def get_updates(dt):
    """Retrieves the vehicle updates for the given date from
    the archive site and read them into a pandas DataFrame.
    
    :param dt: a date-like object
    """
    # This reads the gzipped tarfile into memory, extracts the file,
    # and reads it directly into a dataframe.
    
    # Pandas can also read directly from URLs and file paths
    return pd.read_csv(archive_url(dt), dtype="unicode")

def download_updates(dt, verbose=True):
    """A convenience function that downloads and creates a DataFrame
    for data from the given datetime `dt`, with visual feedback for
    download progress and local caching.
    
    Creates a local directory called 'data' for storing downloaded data
    files. Checks there first for data before attempting to download.
    
    """
    os.makedirs("data", exist_ok=True)
    url = archive_url(dt)
    basename = os.path.basename(url)
    outfile = os.path.join("data", basename)
    log = lambda *m: verbose and print(*m)
    if not os.path.exists(outfile):
        log(f"Downloading {basename}...")
        retrieve(url, outfile)
        log("Download complete")
    else:
        log(f"{basename} exists")
    log("Reading file...")
        
    return pd.read_csv(outfile, dtype="unicode")<jupyter_output><empty_output><jupyter_text>And let's retrieve the updates for February 9th:<jupyter_code># now_utc = pytz.utc.localize(datetime.utcnow())
# yesterday = now_utc.astimezone(pytz.timezone("US/Eastern")) - timedelta(days=3)

day = pytz.utc.localize(datetime(2018, 2, 9))

df = download_updates(day)
df.drop_duplicates(inplace=True)
df.head()<jupyter_output>2018-02-09.csv.gz exists
Reading file...<jupyter_text>We can use the `DataFrame`'s `describe()` method to show summary statistics for each column.<jupyter_code>desc = df.describe()
desc<jupyter_output><empty_output><jupyter_text>Note that pandas is treating all the columns as categorical columns. Instead of getting the median, mean, max (&c) that we might expect, we see the count of non-empty values (`count`); the number of unique values in each column; the most frequent value for each column (`top`); and the number of times the top column occurs (`freq`).The reason is that when we imported the CVS file, we specified a dtype of `"unicode"`, which prevents pandas from inferring the type of every column. On this dataset, certain columns can look numeric (`vehicle_id`, `stop_id`) when they ought—according to the GTFS spec—to be treated as strings. (Note: if we knew the column types in advance, we could save a lot of memory by providing converters or types at import time. I prefer `"unicode"` for exploration, since it behaves predictably)  DataFrameBriefly, a `DataFrame` consists of a number of `Series` with a shared index. You can think of a `Series` as a column of data with some extra metadata attached--a name and an index, for example--and lots of helper methods.In our data, the `stop_sequence` column is an example of a `Series`.We can pick out a single `Series` by name:<jupyter_code>df["stop_sequence"].head()<jupyter_output><empty_output><jupyter_text>And we can summarize it, if we'd like:<jupyter_code>df["stop_sequence"].describe()<jupyter_output><empty_output><jupyter_text>Note that when we run `describe()` on a `Series`, we get another `Series` back, with an index that should look familiar. Indexes can take many forms; among other things, they can be numeric, categorical, or hierarchical. In the case of the summary `Series`, it's a named index containing "count", "unique", "top", and "freq". So when we ran `describe()` on the DataFrame earlier and got another `DataFrame` back, pandas effectively ran `describe()` on each column (that is, each `Series`), then bundled all the summary `Series` together, using their shared index for the row names in the summary `DataFrame` we got back.  Interlude: VectorizationOperations on core pandas data structures like `DataFrame` and `Series` are **vectorized**. That means, in general, that when you perform an operation on a structure, it is applied across that structure's elements. What that looks like concretely depends on the operator (or method) and the operands.Let's see a few quick examples using randomly generated dummy data:<jupyter_code># 5 x 5 DataFrame of pseudorandom integers between 1-99
from ipywidgets import HTML
n_df = pd.DataFrame(np.random.randint(1, 100, size=(5, 5)), columns=list("abcde"))
n_df2 = pd.DataFrame(np.random.randint(1, 100, size=(5, 5)), columns=list("abcde"))
display(n_df, n_df2)<jupyter_output><empty_output><jupyter_text>Adding two dataframes: element-wise addition **by index**.<jupyter_code>n_df + n_df2<jupyter_output><empty_output><jupyter_text>Multiplying by a scalar:<jupyter_code>5 * n_df<jupyter_output><empty_output><jupyter_text>Element-wise comparison with a scalar:<jupyter_code>n_df > 50
# Element-wise multiplication on each row
rvec = np.random.randint(1, 100, size=5)
print(rvec)
n_df * rvec<jupyter_output>[57 46 50 31 34]<jupyter_text>We can also apply aggregation operations. Often, these will reduce the dimensions. For instance, summing a vector will result in a scalar. Summing a `DataFrame` will produce a `Series`.By default, a dataframe applies functions across columns. You can change the axis using the `axis` parameter. `0` means columns, and `1` means rows. (You can also specify `"index"` for `0` or `"columns"` for `1`. That's not a typo. It's specifying the axis _along which_ the aggregate is calculated. For example, `axis=0` for `sum` means that we go *across* the columns and calculate the sum *along* the rows of that column. Frankly, this is why I prefer `0` and `1`; you can map those to your preferred way of thinking about it.)<jupyter_code>n_df.sum()
n_df.sum(axis=1)<jupyter_output><empty_output><jupyter_text>The generalization of applying an operation along an axis is the `apply` method.<jupyter_code>n_df.apply(lambda column: column.sum())<jupyter_output><empty_output><jupyter_text>The generalization of applying an element-wise operation is `applymap`.<jupyter_code>n_df.applymap(lambda n: 1000+n)
# Or: n_df.applymap((1000).__add__)<jupyter_output><empty_output><jupyter_text>`Series` objects behave similarly.<jupyter_code>a_column = n_df["a"]
display(a_column)
a_column.sum()
# Adding two columns (Series with a shared index)
n_df["a"] + n_df["b"]
# Or two rows (note that the index corresponds to the column names in the DataFrame)
n_df.loc[0] + n_df.loc[1]<jupyter_output><empty_output><jupyter_text>Cleaning upAs we saw above, we can access an individual column as a `Series`:<jupyter_code>df["stop_sequence"].head()<jupyter_output><empty_output><jupyter_text>And we can get descriptive statistics for the column:<jupyter_code>df["stop_sequence"].describe()<jupyter_output><empty_output><jupyter_text>We can also select a subset of the columns as a DataFrame:<jupyter_code>df[["stop_sequence", "stop_id"]].head()<jupyter_output><empty_output><jupyter_text>Or we can _replace_ a column in a DataFrame. For instance, we can convert the `stop_sequence` column, which is currently a string, to an integer column. ([spec](https://developers.google.com/transit/gtfs/reference/stop_timestxt))<jupyter_code>df["stop_sequence"] = df.stop_sequence.astype(int)<jupyter_output><empty_output><jupyter_text>And we can convert the coordinates to floats:<jupyter_code>df[["lat", "lon"]] = df[["lat", "lon"]].astype(float)<jupyter_output><empty_output><jupyter_text>If we then run `describe` again...<jupyter_code>df.describe()<jupyter_output><empty_output><jupyter_text>... we **only** see a summary of the columns we've just cleaned up.Pandas will, by default, summarize just the numerical columns of the DataFrame. Only if there are no numerical columns will it fall back to the behavior we saw previously.  Timestamp CleanupNow let's clean up the `timestamp` column.We imported the data from CSV, which doesn't have any standard representation for dates. We need to convert the string timestamps to datetimes that can be used in comparisons.The recorded timestamps are in UTC. UTC doesn't observe DST, so a UTC timestamp unambiguously identifies a moment in time. But we're interested in Eastern time, EDT or EST, so we'll need to first turn the strings into UTC datetimes, then convert them into Eastern.We'll put the resulting Series in a new column of the Dataframe.(If you wanted to generalize this process to other systems that are not on Eastern time, you can use the `agency_timezone` field from `agencies.txt` in the [GTFS static feed](https://developers.google.com/transit/gtfs/reference/agencytxt). More on this in a bit.)<jupyter_code>df["local_timestamp"] = pd.to_datetime(df.timestamp).dt.tz_localize("UTC").dt.tz_convert("US/Eastern")
df.head()<jupyter_output><empty_output><jupyter_text>FilteringThe MBTA's feed may publish more than one update for each trip and stop, because the position (lat/lng) and status may vary.We can pick out just the updates that represent vehicle stops in a couple ways.We can use a mask, or we can use the `.query()` method.<jupyter_code>stops = df[df.status == "STOPPED_AT"] # mask
stops.head()<jupyter_output><empty_output><jupyter_text>What's going on here?Well, the equality comparison with `"STOPPED_AT"` is vectorized, so `df.status == "STOPPED_AT"` gives us back a `Series` containing `True` or `False` at each index. When you pass that Series inside square brackets to a dataframe, it acts as a filter.<jupyter_code>stops = df.query("status == 'STOPPED_AT'") # query
stops.head()<jupyter_output><empty_output><jupyter_text>MergingWe have a `trip_id` column, but what we're often interested in is **routes**, not **trips**. The terminology is explained thoroughly in the [GTFS Reference](https://developers.google.com/transit/gtfs/reference/), but briefly, a trip is a collection of stops and times associated with a schedule and a route, and ordered by stop sequence.So, the Red Line train to Ashmont leaving Alewife on non-holiday weekdays at 6:00 AM is a trip.As you may have noted, the trip updates currently omit route_ids. In order to connect our trip_ids to route_ids, we have to join our vehicle update data to the static data provided by the MBTA.To access the static data, we'll need to first define some additional functions. The MBTA publishes the static schedule information in the form of a zipped directory of CSV files. (The files also include some other things we won't be using today, like fares). If you're curious about the full contents of the manifest, you can download the current one [here](http://www.mbta.com/uploadedfiles/MBTA_GTFS.zip) or read the [docs](https://developers.google.com/transit/gtfs/). For the purposes of this tutorial, we'll manipulate the zipfile in memory.The MBTA updates the file several times a year, so you should be careful when analyzing historical data to retrieve the feed corresponding to the date(s) whose data you're viewing. We'll define a function `mbta_feed_url_for` that takes a datetime and returns the URL of the feed that was active on that date, based on the [list of archived feeds](https://www.mbta.com/gtfs_archive/archived_feeds.txt). Below, we'll define `get_zip`, which we'll use to load a GTFS zip file from a URL into memory and `get_zip_item`, which will take a  zip and the name of an entry (`"agency"`, `"stops"`, etc.) and return a dataframe from the contents of the corresponding CSV file.*There are also some functions here for finding the appropriate feed URL for a specific day. For the purposes of this tutorial, you can safely ignore them. The reason they're here is that the MBTA frequently updates the static feed to account for scheduled station closures, detours, new trips, and other changes. Fortunately, they also publish a [CSV file](https://www.mbta.com/gtfs_archive/archived_feeds.txt) documenting all the historical feed URLs and the date range for which each is valid. `mbta_feed_urls_for` and `mbta_feed_url_for` return the appropriate URLs or URL for a date range or single date, respectively.*<jupyter_code>import csv
import zipfile

def mbta_feed_urls():
    """Returns a generator of (feed_start_date, feed_end_date, archive_url)
    tuples from the MBTA's archived feeds site.
    """
    u = urlopen("https://www.mbta.com/gtfs_archive/archived_feeds.txt")
    for l in csv.DictReader(TextIOWrapper(BytesIO(u.read()))):
        yield (datetime.strptime(l["feed_start_date"], "%Y%m%d"),
               datetime.strptime(l["feed_end_date"], "%Y%m%d"),
               l["archive_url"])

def mbta_feed_urls_for(range_start=None, range_end=None):
    eastern = pytz.timezone("US/Eastern")
    range_start = range_start or datetime.now()
    range_end = range_end or range_start
    for start, end, url in mbta_feed_urls():
        if range_start.tzinfo:
            start = start.astimezone(eastern)
            end = end.astimezone(eastern)
        if start <= range_end:
            if end < range_start:
                continue
            yield url
            
def mbta_feed_url_for(when):
    """Get the URL for the MBTA's GTFS feed active at the 
    datetime `when`."""
    return next(mbta_feed_urls_for(when, when), None)

def get_zip(url="http://www.mbta.com/uploadedfiles/MBTA_GTFS.zip"):
    with urlopen(url) as u:
        return zipfile.ZipFile(BytesIO(u.read()))
    
def get_zip_item(feed, name):
    data = TextIOWrapper(BytesIO(feed.read(name + ".txt")), 
                         encoding="utf-8", line_buffering=True)
    # Specify 'unicode' as the datatype so that Pandas
    # doesn't try to infer cell types. It will get it
    # wrong anyway.
    return pd.read_csv(data, dtype="unicode")
static_feed = get_zip(mbta_feed_url_for(day))<jupyter_output><empty_output><jupyter_text>Let's check out `trips.txt` first.<jupyter_code>trips = get_zip_item(static_feed, "trips")
trips.head()<jupyter_output><empty_output><jupyter_text>We want to match up the `trip_id` from each row in our `stops` DataFrame with a row in the `trips.txt` that has the same `trip_id` and add columns back into our data.In Pandas, we do this with `pd.merge`. It's a lot like a **join** in SQL, and as with joins, there are a few different ways to do it, depending on how you want to handle rows for which there is no match.We'll be using the default join method, namely `"inner"`, which will omit rows in both `stops` (our left-hand DataFrame) and `trips` (right-hand) for which there is no match.Briefly, the other types of merge are:`"left"`  preserves all rows in our left-hand DataFrame, and insert nulls for the new, right-hand columns when there is no match. We'll usually find a few unknown trip_ids in `stops`, if there's a detour, for instance, or track work.`"right"`, predictably, preserves all the rows in our right-hand DataFrame. Because we're looking at data for one day, there will be more trips in `trips.txt` that we won't see in `stops`. Some trips run only on weekdays, or holidays, or Fridays, and so on.And `"outer"` just preserves all rows from the left-hand and right-hand frames, inserting nulls in right-hand columns where the left join column has no match, and vice versa.<jupyter_code>stops_with_routes = pd.merge(stops, trips[["trip_id", "route_id", "direction_id"]], on="trip_id")
stops_with_routes.head()<jupyter_output><empty_output><jupyter_text>We can get stop names by merging in data from `stops.txt` in the static feed.Here's a peek at the structure of `stops`:<jupyter_code>stops_info = get_zip_item(static_feed, "stops")
stops_info.head()<jupyter_output><empty_output><jupyter_text>Let's perform the merge.  (We're also going to drop a few columns now, since I don't plan to use them. Obviously, this may not always be a good idea.)<jupyter_code>named_stops = pd.merge(stops_with_routes.drop(["lat", "lon", "timestamp", "status", "vehicle_id"], axis=1), 
                       stops_info[["stop_id", "stop_name"]], on="stop_id")
named_stops.head()<jupyter_output><empty_output><jupyter_text>We're going to compare the actual arrival times to the scheduled arrival times.To accomplish that, we'll first need to merge in the `"stop_times.txt"`.<jupyter_code>stop_times = get_zip_item(static_feed, "stop_times")
stop_times.head()<jupyter_output><empty_output><jupyter_text>We'll merge on `trip_id` and `stop_sequence` and add the `arrival_time` column.<jupyter_code>stop_times["stop_sequence"] = pd.to_numeric(stop_times["stop_sequence"])
stops_with_times = pd.merge(named_stops, stop_times[["trip_id", "stop_sequence", "arrival_time"]], on=["trip_id", "stop_sequence"])
stops_with_times.head()<jupyter_output><empty_output><jupyter_text>Let's merge the schedule information with the observed data. Join the rows together wherever the `trip_id` and `stop_sequence` are the same. If everything always went according to plan, we would have the same number of stops after the join as we did before. That's usually not the case. Temporary service disruptions mean that some of the trip ids present in our original `stops` dataframe won't correspond to trips in the static feed. The MBTA provides another API, the [realtime API](http://realtime.mbta.com/Portal/Home/Documents), where you can retrieve information about those missing trips. Implementing that here would be a distraction, so we'll make do with what we have.  More CleanupMaybe you noticed that the values in the  `arrival_time` column are recorded as wall-clock times in the format `hh:mm:ss`. For instance, `03:27:00` for 3:27 am.That makes sense; the schedule doesn't give us exact timestamps, because trips typically recur. Some run every weekday. Some run every Sunday. And so forth.But it does mean we have to convert them to timestamps if we want to calculate delays. Per the spec, the times are relative to noon - 12 hours on the service day, in order to account for days when DST begins or ends.<jupyter_code>from datetime import datetime, timedelta
import pytz

timezone = pytz.timezone("US/Eastern")

def get_date(trip_start):
    y, M, d = map(int, trip_start.split("-"))
    return timezone.localize(datetime(y, M, d).replace(hour=12)) - timedelta(hours=12)

def convert_clock_time(row, start):
    h, m, s = map(int, row.arrival_time.split(":", 2))
    return timezone.normalize(
        start + timedelta(hours=h, minutes=m, seconds=s)
    )<jupyter_output><empty_output><jupyter_text>Now let's use the function to convert the scheduled arrival time to a timestamp, using the trip start date. It may take a few moments.<jupyter_code>trip_start = stops_with_times.trip_start.unique()[0]
scheduled_times = stops_with_times.apply(convert_clock_time, axis=1, start=get_date(trip_start))
stops_with_schedule = stops_with_times.assign(scheduled_arrival_time=scheduled_times)
stops_with_schedule.head()<jupyter_output><empty_output><jupyter_text>We can calculate the delay now simply by subtracting the scheduled arrival time from the local timestamp:<jupyter_code>delay = (stops_with_schedule.local_timestamp - stops_with_schedule.scheduled_arrival_time).dt.total_seconds()
# The .dt accessor on the Series provides access to vectorized datetime methods.
delay.head()
delay.describe()
stops_with_schedule["delay"] = delay
routes = stops_with_schedule.groupby("route_id")
# We could also use `.describe()` for full descriptive statistics
routes["delay"].median()
red_line = stops_with_schedule.query("route_id == 'Red'")
red_line.head()<jupyter_output><empty_output><jupyter_text>We got a brief glimpse of `GroupBy` above. We can use `apply` on them as well. If we think of `Series.apply` apply a function down the length of a Series and calling it on each member; and we think of `DataFrame.groupby` as traversing the columns or rows of the frame; `GroupBy.apply` applies a function across each group. So `Series.apply` takes a function that receives a scalar input; `DataFrame.apply` takes a function that receives a `Series`; and `GroupBy.apply` takes a function that receives a `DataFrame` (technically, it depends what you're grouping).Let's create an 'hour_start' column for our red line data. This will represent the hour during which each trip began. To do that, we'll have to group by `trip_id`.<jupyter_code>def add_hour_start(group):
    start_time = group.scheduled_arrival_time.min()
    group["hour_start"] = start_time.hour
    group["schedule_offset"] = (group.scheduled_arrival_time - start_time).dt.total_seconds()
    return group

red_line = red_line.groupby("trip_id").apply(add_hour_start)
red_line.head()
orange_line = stops_with_schedule.query("route_id == 'Orange'").groupby("trip_id").apply(add_hour_start)
blue_line = stops_with_schedule.query("route_id == 'Blue'").groupby("trip_id").apply(add_hour_start)<jupyter_output><empty_output><jupyter_text>Let's take a look at how many trips start during each hour. We'll be using [matplotlib](http://matplotlib.org), a high-level charting API, and [Seaborn](http://seaborn.pydata.org), which adds styles and additional chart types to matplotlib.<jupyter_code>%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns

red_line.groupby("hour_start").trip_id.nunique().plot(color="red")
orange_line.groupby("hour_start").trip_id.nunique().plot(color="orange")
blue_line.groupby("hour_start").trip_id.nunique().plot(color="blue")<jupyter_output><empty_output><jupyter_text>Now let's visualize the trip delays.<jupyter_code>sns.stripplot(x="hour_start", y="delay", data=red_line)<jupyter_output><empty_output><jupyter_text>Note that we're collapsing all the stops a trip makes into a line, so we can't distinguish between a train that leaves the first station late from a train that starts on time and gets progressively later the further it gets from its origin.We can address this in a few ways.We could take the mean delay of each trip:<jupyter_code>red_line_trips = red_line.groupby("trip_id")
red_line_trip_delays = pd.DataFrame({
    "delay": red_line_trips.delay.mean(), 
    "hour_start": red_line_trips.hour_start.first()})
sns.stripplot(x="hour_start", y="delay", data=red_line_trip_delays)<jupyter_output><empty_output><jupyter_text>This gives an overall impression of how early or late entire trips are, but we can do better. We can use a heatmap to show the average lateness at each stop for each hour:<jupyter_code>def reshape(stops):
    return stops.groupby("stop_sequence").delay.mean()

f, ax = plt.subplots(figsize=(9, 6))
plt.title(day.strftime("%a %b %-d, %Y"))
sns.heatmap(red_line.groupby("hour_start").apply(reshape), 
            ax=ax)
# orange_line.groupby("hour_start").apply(reshape)<jupyter_output><empty_output><jupyter_text>We could also take a closer look at the trips for a particular hour. In particular, let's look at delays plotted against `stop_sequence`.<jupyter_code>red_line_7 = red_line.query("hour_start == 7")
sns.stripplot(x="stop_sequence", 
              y="delay", 
              data=red_line_7)
plt.figure(1)

for trip_id, stops in red_line_7.groupby("trip_id"):
    plt.plot(stops.schedule_offset/60, stops.delay, label=trip_id)

plt.show()
red_line_7_stops = red_line_7.groupby("stop_sequence").delay.describe()
red_line_7_stops.head()
plt.figure(1)
plt.plot(red_line_7_stops["50%"], color="red")
plt.plot(red_line_7_stops["25%"], color="pink")
plt.plot(red_line_7_stops["75%"], color="pink")
plt.show()<jupyter_output><empty_output><jupyter_text>This shows us the bounds of the middle 50 of the delays on red line trains starting between 7 and 8 am.We can show the<jupyter_code>def reshape(stops):
    return stops.groupby("stop_sequence").delay.mean()

f, ax = plt.subplots(figsize=(9, 6))
sns.heatmap(red_line.groupby("hour_start").apply(reshape), 
            ax=ax)
red_line_7.head()<jupyter_output><empty_output><jupyter_text>How about some readable stop names?<jupyter_code>stop_info = get_zip_item(static_feed, "stops")
stop_info.head()
red_line_stops_named = pd.merge(red_line_7, stop_info[["stop_id", "stop_name"]], on=["stop_id"])
red_line_stops_named.head()
def reshape(stops):
    return stops.groupby("stop_sequence").delay.mean()

f, ax = plt.subplots(figsize=(9, 6))
reshaped_red_line_stops_named = red_line_stops_named.groupby("hour_start").apply(reshape)
# sns.heatmap(reshaped_red_line_stops_named, xticklabels=red_line
#             ax=ax
reshaped_red_line_stops_named.stop_name
routes = get_zip_item(static_feed, "routes")
routes.head()
routes[routes.agency_id == 2]
# 
# recent_updates.head()
# trips
# recent_updates.trip_id = recent_updates.trip_id.str.strip()
# pd.merge(recent_updates, trips, on="trip_id")
recent_updates = get_updates(datetime.now() - timedelta(days=1))
recent_updates = pd.merge(recent_updates, trips, on="trip_id")
recent_updates.route_id.nunique()
grouped_routes = recent_updates.groupby("route_id")
distinct_trips = grouped_routes.trip_id.nunique()
distinct_updates = grouped_routes.count().trip_id
route_stats = pd.concat({"trips": distinct_trips, "updates": distinct_updates}, axis=1)
route_stats["updates per trip"] = route_stats["updates"]/route_stats["trips"]
route_stats.sort_values("updates", ascending=False)
orange_line = stops_with_schedule.query("route_id == 'Orange'").groupby("trip_id").apply(add_hour_start)
blue_line = stops_with_schedule.query("route_id == 'Blue'").groupby("trip_id").apply(add_hour_start)<jupyter_output><empty_output><jupyter_text>Vehicle Clusters<jupyter_code>named_stops.sort_values("local_timestamp", inplace=True)
red_line = named_stops.query("route_id == 'Red'")
orange_line = named_stops.query("route_id == 'Orange'")
blue_line = named_stops.query("route_id == 'Blue'")
# Group by stop_id and direction_id?

rl_grouped = red_line.groupby(["stop_id", "direction_id"])

def add_headway(df):
  df["headway"] = (df.local_timestamp.diff(1)).dt.total_seconds()
  return df

rl_with_headway = rl_grouped.apply(add_headway)
rl_with_headway.tail()
one_stop = rl_with_headway.query("stop_sequence == 100 and direction_id == '1'")
one_stop.head(10)<jupyter_output><empty_output><jupyter_text>Let's see if we can find clusters at this stop.We'll define a cluster as one or more vehicles that have less than a 300 second (five minute) interval between arrivals at a given stop. We'll give each cluster a unique number, starting at 0.We'll take advantage of the fact that `False` = 0 and `True` = 1, and we'll use the `cumsum()` Series method to get the cumulative sum of a boolean column. The intuition is that each time we exceed the cluster threshold for time between vehicles, we increment `cluster`.<jupyter_code>one_stop["cluster"] = (one_stop.headway >= 300).cumsum()
one_stop<jupyter_output>/opt/conda/lib/python3.6/site-packages/ipykernel_launcher.py:1: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  """Entry point for launching an IPython kernel.<jupyter_text>How many clusters do we have?<jupyter_code>print("Total:", len(one_stop))
print("Clusters:", one_stop.cluster.nunique())<jupyter_output>Total: 82
Clusters: 80<jupyter_text>In other words, there are only two clusters of two or more vehicles. Pretty boring.We can quickly experiment with some other threshold headways:<jupyter_code>print("300 second headway:", sum(one_stop.headway < 300)) # Note that we invert the condition
print("360 second headway:", sum(one_stop.headway < 360))
print("420 second headway:", sum(one_stop.headway < 420))<jupyter_output>300 second headway: 2
360 second headway: 6
420 second headway: 10<jupyter_text>We can see how many clusters we get at each stop on the red line for a given threshold;<jupyter_code>one_direction = rl_with_headway.query("direction_id == '1'")
# 3 minute headway:
one_direction.groupby("stop_sequence").apply(lambda grp: pd.Series({"stop_name": grp.stop_name.iloc[0], "clusters": sum(grp.headway < 180)}))
n_df["odd"] = (n_df.a % 2)
n_df
n_df.groupby("odd").size()
# n_df.groupby("odd", as_index=False).size()<jupyter_output><empty_output><jupyter_text>If we want to get a little fancy, we can also look up the maximum timestamp for each `trip_id` and `stop_sequence`. Here, we use a `groupby` operation to create groups where the trip and stop sequence are the same and find the index of the maximum timestamp within each group.The call to `idxmax()`, in conjunction with the groupby, creates a `MultiIndex` that we can then use as a label-based index into the original dataframe.The effect is to select the last update for a stop on a particular trip on the chosen day. Occasionally, vehicles may not send an update with status STOPPED_AT for a particular stop. It's worth considering this as a way to capture more stop times.<jupyter_code>stops = df.loc[df.groupby(["trip_id","stop_sequence"])["local_timestamp"].idxmax()]
stops.head()<jupyter_output><empty_output><jupyter_text>Remove the commuter rail and ferries: What is the relationship of tweet volume/minute and delays?<jupyter_code>!pip install TwitterAPI
from TwitterAPI import TwitterAPI, TwitterRestPager
consumer_key = "iJcZc53yULImIO7P0GteQzkdy"
consumer_secret = "<KEY>"
access_token = "<KEY>"
access_token_secret = "<KEY>"
api = TwitterAPI(consumer_key, consumer_secret, access_token, access_token_secret)
r = api.request("search/tweets", {"q": "@MBTA"})
# [status["text"] for status in r.json()["statuses"]]
r.json()["statuses"][0].keys()
r.json()["statuses"][0]
from datetime import datetime
datetime.strptime("")
r.json()["statuses"][0]["created_at"]
"%a %b %d %H:%M:%S %z %Y"

api.request("statuses/filter", {"track": "@MBTA"})
for item in r:
    print(item)
datetime.strptime("Thu Nov 09 22:42:25 -0000 2017", "%a %b %d %H:%M:%S %z %Y"
                 )
import pytz
help(pytz)<jupyter_output>Help on package pytz:

NAME
    pytz

DESCRIPTION
    datetime.tzinfo timezone definitions generated from the
    Olson timezone database:
    
        ftp://elsie.nci.nih.gov/pub/tz*.tar.gz
    
    See the datetime section of the Python Library Reference for information
    on how to use these modules.

PACKAGE CONTENTS
    exceptions
    lazy
    reference
    tzfile
    tzinfo

CLASSES
    builtins.Exception(builtins.BaseException)
        pytz.exceptions.InvalidTimeError
            pytz.exceptions.AmbiguousTimeError
            pytz.exceptions.NonExistentTimeError
    builtins.KeyError(builtins.LookupError)
        pytz.exceptions.UnknownTimeZoneError
    
    class AmbiguousTimeError(InvalidTimeError)
     |  Exception raised when attempting to create an ambiguous wallclock time.
     |  
     |  At the end of a DST transition period, a particular wallclock time will
     |  occur twice (once before the clocks are set back, once after). Both
     |  possibilities may be correct,[...]
{'fraction_non_alphanumeric': 0.10276003301036159, 'fraction_numerical': 0.006846593514075252, 'mean_word_length': 5.35671264814455, 'pattern_counts': {'":': 9, '<': 196, '<?xml version=': 0, '>': 193, 'https://': 8, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '48', 'dir': 'jupyter-structured-clean-dedup', 'id': '631052', 'license': 'MIT', 'n_tokens_mistral': 10586, 'n_tokens_neox': 10141, 'n_words': 3442, 'path': 'Crash Course.ipynb', 'repo_name': 'bdougsand/data-tutorial', 'size': '32847'}