0_Taxi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CSS Lab: Extra Topics in CSS\n",
    "\n",
    "## Section 1: Background\n",
    "\n",
    "Much of the data studied in other labs is digital trace data from online platforms. However, computational methods can be used to study many types of data. In this lab, you will examine records of taxi trips, compiled by the New York City [Taxi and Limousine Commission](https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page). Many cities compile data on city functions and make them available through APIs and open data portals.\n",
    "\n",
    "Taxi trips vary from data studied so far because they combine space, time, and network data. In NYC, geographic ___zones___ are used to classify the pick-up and drop-off location of each trip. Taxi trips form a network of interconnections between zones, with each pick-up and drop-off representing a point in space and time.\n",
    "\n",
    "In this lab, you will combine the taxi data with US Census data to study how socioeconomic factors relate to interconnections between taxi zones.\n",
    "\n",
    "## Section 2: Imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- These cells take a while to run! Probably 10 mins for each of the following 5 cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!conda install -c conda-forge cartopy -y -q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install descartes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!conda install -c conda-forge geopandas -y -q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install Fiona"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip uninstall shapely -y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install shapely==1.8.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bz2\n",
    "import cartopy.io.shapereader as shpreader\n",
    "from datetime import datetime\n",
    "from descartes import PolygonPatch\n",
    "import fiona\n",
    "import io\n",
    "import math\n",
    "import matplotlib as mpl\n",
    "from matplotlib.collections import PatchCollection\n",
    "import matplotlib.lines as mlines\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import requests\n",
    "import scipy.stats as spstats\n",
    "import shapefile\n",
    "from shapely.geometry import MultiPoint, MultiPolygon, Polygon, Point, shape\n",
    "import sys\n",
    "from tqdm import tqdm\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 3: Data\n",
    "### 3.1 American Community Survey Data\n",
    "- The American Community Survey (ACS) is a part of the US Census that gives us more detailed information about the US population. We will use some ACS data (population, household income, and fraction of non-hispanic white people in population) to characterize the neighborhoods where taxis go. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load raw census data\n",
    "df_census = pd.read_csv('census.tsv', delimiter='\\t')\n",
    "\n",
    "# clean up the data\n",
    "to_rename = {\n",
    "    'Families: Median Family Income in the Past 12 Months (In 2016 Inflation-Adjusted Dollars)': 'income',\n",
    "    'Total Population.1': 'population',\n",
    "    'Not Hispanic or Latino: White Alone': 'non_hispanic_white',\n",
    "}\n",
    "\n",
    "df_census = df_census.rename(columns=to_rename)\n",
    "df_census['index'] = df_census.FIPS\n",
    "df_census = df_census.set_index('index')\n",
    "df_census['non_hispanic_white_frac'] = df_census['non_hispanic_white'] / df_census['population']\n",
    "df_census = df_census[['FIPS', 'income', 'non_hispanic_white_frac']]\n",
    "df_census = df_census.dropna()\n",
    "\n",
    "#reduce memory use by saving data in smallest dtypes\n",
    "icols = ['FIPS','income']\n",
    "for c in icols:\n",
    "    df_census[c] = pd.to_numeric(df_census[c], downcast='integer')\n",
    "df_census['non_hispanic_white_frac'] = pd.to_numeric(df_census.non_hispanic_white_frac, downcast='float')\n",
    "\n",
    "df_census.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Taxi Trip Data\n",
    "- Here we load data on all yellow taxi trips in NYC between January 1 and January 17, 2017. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_taxi = pd.read_csv('taxi.csv.gz')\n",
    "df_taxi.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# clean up data \n",
    "df_taxi = df_taxi[df_taxi.tpep_pickup_datetime.notna()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reduce memory use by saving data in the smallest types\n",
    "icols = ['DOLocationID','PULocationID']\n",
    "for c in icols:\n",
    "    df_taxi[c] = pd.to_numeric(df_taxi[c], downcast='integer')\n",
    "df_taxi['trip_distance'] = pd.to_numeric(df_taxi['trip_distance'], downcast='float')\n",
    "df_taxi['tpep_pickup_datetime'] = pd.to_datetime(df_taxi.tpep_pickup_datetime)#format='%Y-%m-%d %H:%M:%S'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select just the first half of the month's data reduce memory use\n",
    "df_taxi = df_taxi.sort_values(by='tpep_pickup_datetime')\n",
    "df_taxi = df_taxi.head(n=int(df_taxi.shape[0]/2))\n",
    "df_taxi.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Load NYC Taxi Zones\n",
    "- Here we load data on the geographical boundaries of taxi zones, or neighborhoods, in NYC. This will let us draw maps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "def log_ticks(low, high):\n",
    "    lowexp = math.floor(low)\n",
    "    highexp = math.ceil(high)\n",
    "    ticks = []\n",
    "    labels = []\n",
    "    for e in range(lowexp, highexp + 1):\n",
    "        for i in range(1, 10):\n",
    "            v = i * 10**e\n",
    "            logv = np.log10(v)\n",
    "            if logv > lowexp and logv < highexp:\n",
    "                ticks.append(logv)\n",
    "                labels.append(v)\n",
    "            if highexp - lowexp > 2:\n",
    "                break\n",
    "    return ticks, labels\n",
    "\n",
    "def plot_geo_values(zones, key=None, label=None, transform=None, cmap='Blues'):\n",
    "    \n",
    "    try:\n",
    "        #backward compatability with old cartopy\n",
    "        for z in zones:\n",
    "            z._geometry = z.geometry[0]\n",
    "\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    if transform is None:\n",
    "        tf = lambda x: x\n",
    "        itf = lambda x: x\n",
    "    elif transform == 'log':\n",
    "        tf = np.log10\n",
    "        itf = lambda x: np.power(10, x)\n",
    "    if key is not None:\n",
    "        values = [r.attributes[key] for r in zones]\n",
    "        if transform == 'log':\n",
    "            # Drop zeros for log transform\n",
    "            values = [v for v in values if v > 0]\n",
    "        count_max, count_min = tf(max(values)), tf(min(values))\n",
    "        count_span = count_max - count_min\n",
    "    patches = []\n",
    "    cm = plt.get_cmap(cmap)\n",
    "    for record in zones:\n",
    "        if key is not None:\n",
    "            try:\n",
    "                raw_value = record.attributes[key]\n",
    "                if transform == 'log' and raw_value == 0:\n",
    "                    # Skip 0 values for log\n",
    "                    continue\n",
    "                value = tf(raw_value)\n",
    "            except KeyError:\n",
    "                if transform == 'log':\n",
    "                    continue\n",
    "                else:\n",
    "                    value = 0\n",
    "            color = cm((value - count_min) / count_span)\n",
    "        else:\n",
    "            color = \"#efefef\"\n",
    "        patch = PolygonPatch(record._geometry, fc=color, ec='#555555', lw=0.2, alpha=1, zorder=1)\n",
    "        patches.append(patch)\n",
    "\n",
    "    mp = MultiPolygon([\n",
    "        shape(record._geometry)\n",
    "        for record in zones])\n",
    "    minx, miny, maxx, maxy = mp.bounds\n",
    "    w, h = maxx - minx, maxy - miny\n",
    "    aspect = w / h\n",
    "\n",
    "    fig = plt.figure(figsize=(8 * aspect / 0.9, 8))\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.set_xlim(minx - 0.2 * w, maxx + 0.2 * w)\n",
    "    ax.set_ylim(miny - 0.2 * h, maxy + 0.2 * h)\n",
    "    ax.set_aspect(1)\n",
    "    ax.add_collection(PatchCollection(patches, match_original=True))\n",
    "    \n",
    "    if key is not None:\n",
    "        # Add colorbar\n",
    "        plt.subplots_adjust(right=0.9)\n",
    "        cax = plt.axes([0.95, 0.1, 0.035, 0.8])\n",
    "        norm = mpl.colors.Normalize(vmin=count_min, vmax=count_max)\n",
    "        cb1 = mpl.colorbar.ColorbarBase(cax, cmap=cm,\n",
    "                                        norm=norm,\n",
    "                                        orientation='vertical')\n",
    "        if label is not None:\n",
    "            cb1.set_label(label)\n",
    "        else:\n",
    "            cb1.set_label(key)\n",
    "        if transform == 'log':\n",
    "            ticks, labels = log_ticks(count_min, count_max)\n",
    "            cb1.set_ticks(ticks)\n",
    "            cb1.set_ticklabels(labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reader = shpreader.Reader('taxi/taxi_zones/taxi_zones-latlong.shp')\n",
    "taxi_geo = list(reader.records())\n",
    "plot_geo_values(taxi_geo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Load Census Tracts\n",
    "- here we load data from the US census on geographical areas, called \"census tracts.\" __Note that census tracts are different than taxi zones__. In general, every data source that does geographic analysis or mapping uses its own way of dividing up space, which makes working with geography especially hard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reader = shpreader.Reader('tracts/cb_2017_36_tract_500k.shp')\n",
    "tract_geo = list(reader.records())\n",
    "\n",
    "# Keep tracts in NYC\n",
    "geoids = set(df_census['FIPS'])\n",
    "tract_geo = [\n",
    "    record for record in tract_geo\n",
    "    if int(record.attributes['GEOID']) in geoids]\n",
    "\n",
    "plot_geo_values(tract_geo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 4. Visualize Data\n",
    "#### Count Pick-Ups and Drop-Offs\n",
    "\n",
    "Let's start by visualizing the data. First, we use a helper function to count the number of pick-ups and drop-offs in each taxi zone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "def annotate_records(records, data, key):\n",
    "    for r in records:\n",
    "        geoid = int(r.attributes['GEOID'])\n",
    "        r.attributes[key] = data[geoid]\n",
    "\n",
    "def update_counts(df, records):\n",
    "    print('Updating counts')\n",
    "    counts_pickup = df.groupby('PULocationID').count()['tpep_pickup_datetime']\n",
    "    counts_dropoff = df.groupby('DOLocationID').count()['tpep_pickup_datetime']\n",
    "    zone_pu_counts = {}\n",
    "    zone_do_counts = {}\n",
    "    for zone, count in dict(counts_pickup).items():\n",
    "        zone_pu_counts[zone] = count\n",
    "    for zone, count in dict(counts_dropoff).items():\n",
    "        zone_do_counts[zone] = count\n",
    "    for r in records:\n",
    "        geoid = r.attributes['LocationID']\n",
    "        r.attributes['pickups'] = zone_pu_counts.get(geoid, 0)    \n",
    "        r.attributes['dropoffs'] = zone_do_counts.get(geoid, 0)    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "update_counts(df_taxi, taxi_geo)\n",
    "taxi_geo[0].attributes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualize Pick-Ups\n",
    "The following cell plots a map of New York taxi zones, with darker zones representing more taxi pick-ups. The number of pick-ups exhibits a heavy tail, so the colors are chosen based on the logarithm of the counts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_geo_values(taxi_geo, 'pickups', 'Pick-Ups', transform='log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Remove Airports\n",
    "In the above visualization, the taxi zones representing airports stand out as outliers. We can remove them from the analysis and update trip counts using the following code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify airport zones\n",
    "airports = [r for r in taxi_geo if 'Airport' in r.attributes['zone']]\n",
    "airport_zones = [r.attributes['LocationID'] for r in airports]\n",
    "\n",
    "# Remove airports zones and trips\n",
    "taxi_geo = [r for r in taxi_geo if 'Airport' not in r.attributes['zone']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_taxi = df_taxi[\n",
    "    (~df_taxi['DOLocationID'].isin(airport_zones))\n",
    "    & (~df_taxi['PULocationID'].isin(airport_zones))]\n",
    "\n",
    "# Update counts and plot\n",
    "update_counts(df_taxi, taxi_geo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_geo_values(taxi_geo, 'pickups', 'Pick-Ups', transform='log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\">\n",
    "    \n",
    "#### Short Answer 1\n",
    "Whether or not we want to include airports depends on the type of questions we want to answer. Most people traveling to and from an airport do not live or work at that airport. For each of the following two examples, explain why you would or would not include trips to and from airports in an analysis:\n",
    "\n",
    "1. Understanding wear and tear on roads.\n",
    "2. Understanding transportation between residential communities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualize income\n",
    "Now let's visualize household income. This data is organized into census tracts, which are typically smaller than taxi zones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "# Get geoid->values from df for each element in geo\n",
    "def get_geo_map(df, key, geo):\n",
    "    result = {}\n",
    "    for i, r in enumerate(geo):\n",
    "        geoid = int(r.attributes['GEOID'])\n",
    "        result[geoid] = df.loc[geoid, key]\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tract_income_map = get_geo_map(df_census, 'income', tract_geo)\n",
    "annotate_records(tract_geo, tract_income_map, 'income')\n",
    "plot_geo_values(tract_geo, 'income', transform='log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "#### Short Answer 2\n",
    "1. Do the areas with **high income** tend to have a higher than average, near average, or below average number of taxi pick-ups? \n",
    "- (Hint: compare the income map you just made with the trip maps you made earlier.)\n",
    "\n",
    "2. What is one way income level could influence the number of taxi trips to or from a neighborhood? \n",
    "    \n",
    "3. Can you think of a third variable that may be confounding the relationship between income and the number of trips? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pairwise Trip Data\n",
    "Rather than individual zone pick-ups, we're really interested in connections _between_ zones. In other words, we want to analyze the pairwise trip data. In urban planning, this kind of data is known as an __Origin-Destination Matrix (OD Matrix)__. Traditionally, this type of data is collected using surveys and used for planning travel capacity (i.e. roads), routing and zoning in urban areas. \n",
    "\n",
    "#### Create data frame for pairwise data\n",
    "The first step is to create a new data frame. The entries of this data frame are indexed by pairs of zones rather than single zones. __Since we are interested in travel across traffic zones, we ignore any trips that may have originated and ended within the same zone.__ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "def make_pair_df(geo):\n",
    "    index = set([\n",
    "        (r.attributes['LocationID'], s.attributes['LocationID'])\n",
    "        for r in geo\n",
    "            for s in geo\n",
    "                if r.attributes['LocationID'] != s.attributes['LocationID']\n",
    "    ])\n",
    "    df_pairs = pd.DataFrame(index=list(index))\n",
    "    return df_pairs\n",
    "\n",
    "def plot_pairs(\n",
    "    geo,\n",
    "    df_pairs,\n",
    "    key,\n",
    "    both_dir=True,\n",
    "    label=None,\n",
    "    bounds=None,\n",
    "    threshold=0,\n",
    "    transform=None,\n",
    "    linewidth=0.1,\n",
    "    alpha=0.3\n",
    "):\n",
    "    points = {}\n",
    "    patches = []\n",
    "    print('Finding representative points')\n",
    "    sys.stdout.flush()\n",
    "    for record in tqdm(geo):\n",
    "        patch = PolygonPatch(record._geometry, fc='#efefef', ec='#555555', lw=0.2, alpha=1, zorder=1)\n",
    "        patches.append(patch)\n",
    "        record_id = record.attributes['LocationID']\n",
    "        points[record_id] = record._geometry.representative_point().coords[0]\n",
    "\n",
    "    mp = MultiPolygon([\n",
    "        shape(record._geometry)\n",
    "        for record in geo])\n",
    "    minx, miny, maxx, maxy = mp.bounds\n",
    "    w, h = maxx - minx, maxy - miny\n",
    "\n",
    "    if bounds is None:\n",
    "        xlim = (minx - 0.2 * w, maxx + 0.2 * w)\n",
    "        ylim = (miny - 0.2 * h, maxy + 0.2 * h)\n",
    "    else:\n",
    "        xlim = (bounds[0], bounds[1])\n",
    "        ylim = (bounds[2], bounds[3])\n",
    "    ar = (xlim[1] - xlim[0]) / (ylim[1] - ylim[0])\n",
    "    \n",
    "    fig = plt.figure(figsize=(8 * ar, 8))\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.set_xlim(xlim)\n",
    "    ax.set_ylim(ylim)\n",
    "    ax.add_collection(PatchCollection(patches, match_original=True))\n",
    "    \n",
    "    if transform is None:\n",
    "        tf = lambda x: x\n",
    "        itf = lambda x: x\n",
    "    if transform == 'log':\n",
    "        tf = lambda x: np.where(x>0, np.log10(x), 0)\n",
    "        itf = lambda x: np.power(10, x)\n",
    "        \n",
    "               \n",
    "    trip_max = tf(df_pairs['trips'].max())\n",
    "    \n",
    "    cm = plt.get_cmap('afmhot')\n",
    "    sys.stdout.flush()\n",
    "    print('Drawing lines')\n",
    "    sys.stdout.flush()\n",
    "    for a_id, b_id in tqdm(df_pairs.index):\n",
    "        # Combine both directions for plot\n",
    "        if a_id >= b_id:\n",
    "            continue\n",
    "        a_point = points[a_id]\n",
    "        b_point = points[b_id]\n",
    "        try:\n",
    "            if both_dir:\n",
    "                raw_count = (\n",
    "                    df_pairs['trips'][(a_id, b_id)]\n",
    "                    + df_pairs['trips'][(b_id, a_id)])\n",
    "            else:\n",
    "                raw_count =  df_pairs['trips'][(a_id, b_id)]\n",
    "            try:\n",
    "                if raw_count < threshold:\n",
    "                    continue\n",
    "            except ValueError:\n",
    "                print(\"raw\", raw_count, 't', threshold, 'a', a_id, 'b', b_id)\n",
    "                raise\n",
    "            with np.errstate(divide='ignore'):\n",
    "                count = tf(raw_count)\n",
    "            f = count / trip_max\n",
    "            color = \"#000000\"\n",
    "            if both_dir:\n",
    "                l = mlines.Line2D(\n",
    "                    [a_point[0], b_point[0]],\n",
    "                    [a_point[1], b_point[1]],\n",
    "                    color='#000000',\n",
    "                    linewidth=linewidth,\n",
    "                    alpha=f*alpha)\n",
    "                ax.add_line(l)\n",
    "            else:\n",
    "                ax.arrow(x=a_point[0],\n",
    "                         y=a_point[1], \n",
    "                         dx=b_point[0] - a_point[0], \n",
    "                         dy=b_point[1] - a_point[1],\n",
    "                         linewidth=linewidth,\n",
    "                         color='#000000',\n",
    "                         alpha=f*alpha,\n",
    "                         length_includes_head=True)   \n",
    "        except RuntimeError:\n",
    "            pass\n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pairs = make_pair_df(taxi_geo)\n",
    "df_pairs.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Count number of trips\n",
    "Now we can count the number of trips between each pair of zones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Helper function\n",
    "\n",
    "# Calculate number of trips\n",
    "def count_trips(df_taxi, df_pairs):\n",
    "    print('Counting trips')\n",
    "    index = df_pairs.index\n",
    "    df_count = df_taxi.groupby(['PULocationID', 'DOLocationID']).count()['tpep_pickup_datetime']\n",
    "    df_pairs['trips'] = np.zeros(len(index))\n",
    "    for a_id, b_id in tqdm(df_pairs.index):\n",
    "        # Sum trips in both directions\n",
    "        count = df_count.get((a_id, b_id), 0)\n",
    "        df_pairs.trips[(a_id, b_id)] = count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count_trips(df_taxi, df_pairs)\n",
    "df_pairs.sort_values('trips', ascending=False).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualize Trips\n",
    "The following visualization draws a line between each pair of zones, with darker lines representing more traffic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "plot_pairs(\n",
    "    taxi_geo, df_pairs, 'trips',\n",
    "    label='Trips',\n",
    "    bounds=(-74.07, -73.83, 40.62, 40.93),\n",
    "    transform='log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Adding a Threshold\n",
    "The above visualization is a little difficult to interpret because there are so many lines. One way to reduce clutter is to add a threshold, only drawing lines between zones with a minimum number of trips."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "threshold = 5000\n",
    "plot_pairs(\n",
    "    taxi_geo, df_pairs, 'trips',\n",
    "    label='Trips',\n",
    "    bounds=(-74.07, -73.83, 40.62, 40.93),\n",
    "    threshold=threshold,\n",
    "    linewidth=1,\n",
    "    alpha=0.5,\n",
    "    transform='log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "#### Short Answer 3\n",
    "\n",
    "Try setting the minimum trip threshold in the cell above to different values.\n",
    "\n",
    "1. What do the resulting networks of taxi trips tell us about the relationships among neighborhoods? State two observations you found by changing the thresholds.\n",
    "    \n",
    "2. Pairwise trip data can be used to answer different questions regarding travel within the city. \n",
    "\n",
    "    (1) For what kind of questions would applying a minimum trip threshold be useful?\n",
    "    \n",
    "    (2) When would it be inappropriate?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Most important trips by zone\n",
    "\n",
    "While it is important to identify where in a city the majority of the traffic is concentrated, sometimes it is also useful to know how well different zones are inter-connected. For example, given any particular zone, we often want to know from which other zones (**generators**) people come to it often and which other zones (**attractors**) people from this zone may be traveling often to. \n",
    "\n",
    "In the following figure, we show only the lines corresponding to the top **n** trips going into (out of) a the dropoff (pickup) zone. Try different values for the type of focal zone and **n** in  the cell below. **Arrows always point from the pickup zone to the dropoff zone.**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_top_pairs = df_pairs[df_pairs.trips > 0].reset_index()\n",
    "df_top_pairs[\"pickup\"] = df_top_pairs.apply(lambda x: x[\"index\"][0],axis=1)\n",
    "df_top_pairs[\"dropoff\"] = df_top_pairs.apply(lambda x: x[\"index\"][1],axis=1)\n",
    "df_top_pairs.index = df_top_pairs[\"index\"]\n",
    "\n",
    "zone_type = \"pickup\" #values: [pickup, dropoff]\n",
    "top_n = 2 # values: 1,2,3,...\n",
    "df_top_pairs[\"rank\"] = df_top_pairs.groupby(zone_type)[\"trips\"].rank(method=\"min\",ascending=False)\n",
    "df_top_pairs = df_top_pairs[df_top_pairs[\"rank\"] <= top_n]\n",
    "plot_pairs(\n",
    "    taxi_geo, df_top_pairs, 'trips',\n",
    "    both_dir=False,\n",
    "    label='Trips',\n",
    "    bounds=(-74.07, -73.83, 40.62, 40.93),\n",
    "    threshold=5,\n",
    "    linewidth=1.0,\n",
    "    alpha=0.3,\n",
    "    transform='log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate Average Trip Distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def get_mean_distance(df_taxi, df_pairs):\n",
    "    print('Finding mean distance')\n",
    "    # Calculate mean trip distance\n",
    "    df_mean = df_taxi.groupby(['PULocationID', 'DOLocationID']).mean()\n",
    "    df_trip = df_mean['trip_distance']\n",
    "    df_pairs['trip_distance'] = np.zeros(len(df_pairs.index))\n",
    "    for a_id, b_id in tqdm(df_pairs.index):\n",
    "        # Get counts in each direction\n",
    "        ab_count = df_pairs['trips'].get((a_id, b_id), 0)\n",
    "        ba_count = df_pairs['trips'].get((b_id, a_id), 0)\n",
    "        # If no trips, average distance is undefined\n",
    "        if (ab_count + ba_count == 0):\n",
    "            df_pairs.at[(a_id, b_id), 'trip_distance'] = float('nan')\n",
    "        # Otherwise take a weighted average of each direction\n",
    "        else:\n",
    "            df_pairs.at[(a_id, b_id), 'trip_distance'] = (\n",
    "               ab_count * df_trip.get((a_id, b_id), 0)\n",
    "                + ba_count * df_trip.get((b_id, a_id), 0)) / (ab_count + ba_count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_mean_distance(df_taxi, df_pairs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(1, 2, 1);\n",
    "df_pairs['trip_distance'].hist(bins=30)\n",
    "plt.xlabel(\"Miles\")\n",
    "plt.ylabel(\"Number of Zone Pairs\")\n",
    "plt.subplot(1, 2, 2);\n",
    "plt.hist(df_taxi[df_taxi['trip_distance'] < 40]['trip_distance'], bins=30)\n",
    "plt.xlabel(\"Miles\")\n",
    "plt.ylabel(\"Number of Trips\")\n",
    "plt.grid()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "#### Short Answer 4\n",
    "\n",
    "Compare the distribution of distances between all possible pairs of traffic zones to the distribution of trip distances. What can you infer about the nature of taxi trips based on this comparison?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compute a map from census tract values to taxi zone\n",
    "- here is where we match up the census tracts with the taxi zones. You don't need to worry how we do this for the lab, but this code may be a useful reference if you want to do similar things in the future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "# Create a function mapping census tract vectors to taxi zone vectors\n",
    "def get_tract_to_zone(taxi_geo, tract_geo):\n",
    "    # Get fractions of zone overlapping each tract\n",
    "    weight_zone_tract, all_zones = get_weight_zone_tract(taxi_geo, tract_geo)\n",
    "    # Create a function to map census data to taxi zones\n",
    "    def tract_to_zone(tract_data):\n",
    "        zone_data = pd.Series(index=list(all_zones))\n",
    "        for zone_id in tqdm(all_zones):\n",
    "            zone_data[zone_id] = 0\n",
    "            for tract_id, d in tract_data.items():\n",
    "                f = weight_zone_tract[zone_id][tract_id]\n",
    "                zone_data[zone_id] += f * tract_data[tract_id]\n",
    "        return zone_data\n",
    "    # Return function\n",
    "    return tract_to_zone\n",
    "\n",
    "# Create a function mapping census tract vectors to taxi zone vectors\n",
    "def get_weight_zone_tract(taxi_geo, tract_geo):\n",
    "    \n",
    "    # Calculate the fraction of each zone that overlaps with a given tract\n",
    "    weight_zone_tract = {}\n",
    "    empty_zones = set()\n",
    "    all_zones = set()\n",
    "    for zone in tqdm(taxi_geo):\n",
    "        zone_id = zone.attributes['LocationID']\n",
    "        all_zones.add(zone_id)\n",
    "        weight_tract = {}\n",
    "        for tract in tract_geo:\n",
    "            tract_id = int(tract.attributes['GEOID'])\n",
    "            # This is a rough approximation that treats lat/long as cartesian\n",
    "            f = zone._geometry.intersection(tract._geometry).area / zone._geometry.area\n",
    "            weight_tract[tract_id] = f\n",
    "        # Exclude taxi zones without enough data\n",
    "        zone_total = sum(weight_tract.values())\n",
    "        if zone_total < 0.5:\n",
    "            empty_zones.add(zone_id)\n",
    "        else:\n",
    "            # Normalize\n",
    "            for tract_id in weight_tract.keys():\n",
    "                weight_tract[tract_id] /= zone_total\n",
    "        weight_zone_tract[zone_id] = weight_tract\n",
    "    return weight_zone_tract, all_zones\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tract_to_zone = get_tract_to_zone(taxi_geo, tract_geo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Taxi Zone Income\n",
    "- Once we have mapped the census tracts to taxi zones, we can estimate average income for people living in each taxi zone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "zone_incomes = tract_to_zone(df_census['income'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "for zone in tqdm(taxi_geo): \n",
    "    zone_id = zone.attributes['LocationID']\n",
    "    zone.attributes['income'] = zone_incomes[zone_id]\n",
    "plot_geo_values(taxi_geo, 'income', label=\"Income\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following histogram shows the distribution of incomes for each taxi zone in the city. Note that some neighborhoods have very high average incomes, while others have very low incomes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(zone_incomes, bins=24)\n",
    "plt.xlabel('Income')\n",
    "plt.ylabel('Number of Zones')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Income Difference\n",
    "- Here we compare the incomes of neighborhoods where taxi trips start with the incomes of the neighborhoods where they end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "# Calculate trips vs income difference\n",
    "def get_log_diff(df_pairs, attribute, values, min_value=0):\n",
    "    print('Finding attribute difference')\n",
    "    key = 'log_{}_diff'.format(attribute)\n",
    "    df_pairs[key] = np.zeros(len(df_pairs.index))\n",
    "    for i, j in tqdm(df_pairs.index):\n",
    "        a = values[i]\n",
    "        b = values[j]\n",
    "        # Ignore zero case and partial zones\n",
    "        if a > min_value and b > min_value:\n",
    "            df_pairs.at[(i, j), key] = np.log10(b) - np.log10(a)\n",
    "        else:\n",
    "            df_pairs.at[(i, j), key] = float('nan')\n",
    "\n",
    "def get_diff(df_pairs, attribute, values, min_value=0):\n",
    "    print('Finding attribute difference')\n",
    "    key = '{}_diff'.format(attribute)\n",
    "    df_pairs[key] = np.zeros(len(df_pairs.index))\n",
    "    for i, j in tqdm(df_pairs.index):\n",
    "        a = values[i]\n",
    "        b = values[j]\n",
    "        if a > min_value and b > min_value:\n",
    "            df_pairs.at[(i, j), key] = b - a\n",
    "        else:\n",
    "            df_pairs.at[(i, j), key] = float('nan')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min(df_census['income'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pairs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_diff(df_pairs, 'income', zone_incomes, min_value=12300)\n",
    "df_pairs['income_diff'].hist(bins=30)\n",
    "plt.xlabel('Difference in income (dollars)')\n",
    "plt.ylabel('Number of pairs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "#### Short Answer 5\n",
    "\n",
    "1. Based on the above figure, what can we say about taxi trips in terms of the incomes of the pickup and dropoff zones?\n",
    "2. Note that in the above figure, we have considered the absolute difference in income between the pickup and destination traffic zones. This means that a case where the dropoff zone has an income **\\$X** higher than the pickup zone and a case where the dropoff zone has an income **\\$X** lower than the pickup zone would both add to the frequency of pairs at the x-axis value **X** in the figure. Remove the code chunk **.abs()** in line 2 of the above cell to change the X axis to difference in income of the dropoff location in comparison to the pickup location (i.e. income at dropoff - income at pickup). What can you say about taxi trips from high income to low income zones vs. trips from low income to high income zones?    \n",
    "3. Give one reason that may explain your observation in 2.\n",
    "    \n",
    "    \n",
    "Note: the Y-axis here represent the number of pairs. A pair (A, B) exists if there's at least one trip from A to B in the dataset. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparing income and connectivity\n",
    "\n",
    "Let's do some systematic analysis to figure out if our intuition from the previous figure that travel is concentrated between zones with similar income levels holds true."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "def plot_binned(x, y, alpha=1, bins=30, log=None, xlabel=None, ylabel=None, title=None):\n",
    "    # Remove nan\n",
    "    isnan = np.isnan(x) | np.isnan(y)\n",
    "    x = np.array(x)[~isnan]\n",
    "    y = np.array(y)[~isnan]\n",
    "    num_bins = bins\n",
    "    low = min(x)\n",
    "    high = max(x)\n",
    "    delta = (high - low) / num_bins\n",
    "    bin_x = [[] for i in range(num_bins)]\n",
    "    bin_y = [[] for i in range(num_bins)]\n",
    "    bin_cx = [low + delta * (i + 1/2) for i in range(num_bins)]\n",
    "    for i in range(len(x)):\n",
    "        try:\n",
    "            b = int(math.floor((x[i] - low) / delta))\n",
    "        except ValueError:\n",
    "            print(x[i] == 'nan', x[i], low, delta)\n",
    "            raise\n",
    "        if b == num_bins:\n",
    "            # Include edge case on top bin\n",
    "            b -= 1\n",
    "        bin_x[b].append(x[i])\n",
    "        bin_y[b].append(y[i])\n",
    "    means = np.array([np.mean(by) for by in bin_y])\n",
    "    std = np.array([np.std(by) for by in bin_y])\n",
    "    count = np.sqrt(np.array([len(by) for by in bin_y]))\n",
    "    se = std / count\n",
    "    low95 = means - 2 * se\n",
    "    high95 = means + 2 * se\n",
    "    if log == 'y' or log == 'xy':\n",
    "        low95[low95 <= 0] = 0.1\n",
    "        high95[high95 <= 0] = 0.1\n",
    "    plt.fill_between(bin_cx, low95, high95, color='#dfdfdf', zorder=2, alpha=alpha)\n",
    "    if log == 'x':\n",
    "        plt.semilogx(x, y, '.', zorder=1, alpha=0.1)\n",
    "        plt.semilogx(bin_cx, means)\n",
    "        plt.semilogx(bin_cx, means, '--', color='#999999')\n",
    "        plt.semilogx(bin_cx, means + 2 * se, '--', color='#999999')\n",
    "    elif log == 'y':\n",
    "        plt.semilogy(x, y, '.', zorder=1, alpha=0.1)\n",
    "        plt.semilogy(bin_cx, means)\n",
    "        plt.semilogy(bin_cx, low95, '--', color='#999999')\n",
    "        plt.semilogy(bin_cx, high95, '--', color='#999999')\n",
    "    elif log == 'xy':\n",
    "        plt.loglog(x, y, '.', zorder=1, alpha=0.1)\n",
    "        plt.loglog(bin_cx, means)\n",
    "        plt.loglog(bin_cx, low95, '--', color='#999999')\n",
    "        plt.loglog(bin_cx, high95, '--', color='#999999')\n",
    "    else:\n",
    "        plt.plot(x, y, '.', zorder=1, alpha=0.1)\n",
    "        plt.plot(bin_cx, means)\n",
    "        plt.plot(bin_cx, means - 2 * se, '--', color='#999999')\n",
    "        plt.plot(bin_cx, means + 2 * se, '--', color='#999999')\n",
    "    if xlabel is not None:\n",
    "        plt.xlabel(xlabel)\n",
    "    if ylabel is not None:\n",
    "        plt.ylabel(ylabel)\n",
    "    if title is not None:\n",
    "        plt.title(title)\n",
    "\n",
    "def filter_borough(df_pairs, borough):\n",
    "    borough_ids = set()\n",
    "    for r in taxi_geo:\n",
    "        if r.attributes['borough'] == 'Manhattan':\n",
    "             borough_ids.add(int(r.attributes['LocationID']))\n",
    "    in_borough = list()\n",
    "    for s, t in df_pairs.index:\n",
    "        if s in borough_ids and t in borough_ids:\n",
    "            in_borough.append(True)\n",
    "        else:\n",
    "            in_borough.append(False)\n",
    "    df_pairs_borough = df_pairs[in_borough]\n",
    "    return df_pairs_borough"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pairs_nyc = filter_borough(df_pairs, 'Manhattan')\n",
    "print(\"Pairs in full set:\", len(df_pairs))\n",
    "print(\"Pairs in Manhattan:\", len(df_pairs_nyc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### See the number of trips between neighborhoods according to the difference in income between the neighborhoods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 8))\n",
    "plt.ylim([0,10000])\n",
    "plot_binned(\n",
    "    df_pairs_nyc['income_diff'], df_pairs_nyc['trips'],\n",
    "    #log='y', \n",
    "    alpha=0.6,\n",
    "    xlabel='Difference in income',\n",
    "    ylabel='Number of Trips')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  See the number of trips between neighborhoods according to the distance between the neighborhoods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate trip distance vs income difference\n",
    "plt.figure(figsize=(8, 8))\n",
    "plot_binned(\n",
    "    df_pairs_nyc['trip_distance'], df_pairs_nyc['trips'],\n",
    "    #log='xy', \n",
    "    alpha=0.8, bins=10,\n",
    "    xlabel='Mean ride distance between zones (miles)',\n",
    "    ylabel='Number of trips between zones')\n",
    "plt.ylim([0, 10000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### But wait, maybe distance and income are related.\n",
    "- Neighborhoods that are close together may have more similar income than neighborhoods that are far apart. \n",
    "- Maybe the relationship we saw between trips and income was just a side effect of the distance between neighborhoods, and not really a result of their incomes. \n",
    "- We can test this by only looking at trips of a specific length, say 1 mile or 7 miles, and seeing whether income influences where people go even when they're going the same distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "def filter_zones(df_pairs, distance, tolerance):\n",
    "    df = df_pairs[(df_pairs.trip_distance - distance).abs() < tolerance]\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 4))\n",
    "for i, miles in enumerate([1, 4, 7]):\n",
    "    df_pairs_miles = filter_zones(df_pairs_nyc, miles, 1.5)\n",
    "    plt.subplot(1, 3, i + 1)\n",
    "    plot_binned(\n",
    "        df_pairs_miles['income_diff'],\n",
    "        df_pairs_miles['trips'],\n",
    "        bins=10, #log='y',\n",
    "        alpha=0.6,\n",
    "        xlabel='Difference in income',\n",
    "        ylabel='Number of Trips',\n",
    "        title='{} Miles'.format(miles))\n",
    "    #plt.ylim([10, 10000])\n",
    "    #plt.xlim([-1, 1])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "#### Short Answer 6\n",
    "\n",
    "The above figures show that taxi trips tend to be concentrated between zones with similar income even after controlling for the distance between the pickup and dropoff. \n",
    "    \n",
    "- Provide one possible mechanism to explain why income is able to explain travel between zones to some extent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that once we control for the distance between the pickup and dropoff zones, travel from high-income zones to low-income zones is more frequent than vice versa. This is especially true at higher distances. Jane Jacobs, one of the most influential writers of modern American urban planning, provided one reason for this phenomenon in her seminal work \"The Death and Life of Great American Cities\"[JJ1961]. Low-income areas in a city are inherently more affordable, particularly, in terms of housing. This, in turn, leads to more diverse residents, unique and niche businesses and culture. This cultural richness attracts visitors from other parts of the city. Unfortunately, in the long-term, this dynamic makes these areas vulnrable to gentrification, increase in cost of living, and the loss of the very cultural capital that made them attractive. \n",
    "\n",
    "Can you think of any other explanations for why travel from high-income zones is more frequent than vice versa? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Time of Day\n",
    "- Another factor affects taxi trips: time of day. 7am taxi trips are probably different than 10pm taxi trips, for example. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\">\n",
    "    \n",
    "#### Short answer 7:\n",
    "- Give three examples of how different times of day might affect who is riding taxis and where they are going. Explain how each example might look in the maps of taxi traffic we made earlier. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "def get_trip_datetime(df):\n",
    "    print('Finding trip times')\n",
    "    \n",
    "    date = pd.to_datetime(df.tpep_pickup_datetime, ) #format='%Y-%m-%d %H:%M:%S'\n",
    "    \n",
    "    df['weekday'] = date.dt.weekday\n",
    "    df['hour'] = date.dt.hour\n",
    "    df['minute'] = date.dt.minute\n",
    "    return\n",
    "\n",
    "def get_time_pairs(df, taxi_geo, zone_incomes):\n",
    "    time_periods = {\n",
    "        'AM Rush': [(5, 10)],\n",
    "        'Midday': [(10, 16)],\n",
    "        'PM Rush': [(16, 20)],\n",
    "        'Night': [(0, 5), (20, 24)]\n",
    "    }\n",
    "    pairs = {}\n",
    "    for label, periods in time_periods.items():\n",
    "        selection = pd.Series(False, index=df.index)\n",
    "        # Check hour of day\n",
    "        for start, stop in periods:\n",
    "            selection = selection | (\n",
    "                (df.hour >= start)\n",
    "                & (df.hour < stop))\n",
    "            print(len(df[selection]))\n",
    "        # M-F only\n",
    "        selection = selection & (df.weekday < 5)\n",
    "        # Select data\n",
    "        df_period = df[selection]\n",
    "        pairs[label] = make_pair_df(taxi_geo)\n",
    "        count_trips(df_period, pairs[label])\n",
    "        get_mean_distance(df_period, pairs[label])\n",
    "        get_log_diff(pairs[label], 'income', zone_incomes, min_value=12300)\n",
    "    return pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "def get_time_pairs(df, taxi_geo, zone_incomes):\n",
    "    time_periods = {\n",
    "        'AM Rush': [(5, 10)],\n",
    "        'Midday': [(10, 16)],\n",
    "        'PM Rush': [(16, 20)],\n",
    "        'Night': [(0, 5), (20, 24)]\n",
    "    }\n",
    "    pairs = {}\n",
    "    for label, periods in time_periods.items():\n",
    "        selection = pd.Series(False, index=df.index)\n",
    "        # Check hour of day\n",
    "        for start, stop in periods:\n",
    "            selection = selection | (\n",
    "                (df.tpep_pickup_datetime.dt.hour >= start)\n",
    "                & (df.tpep_pickup_datetime.dt.hour < stop))\n",
    "            print(len(df[selection]))\n",
    "        # M-F only\n",
    "        selection = selection & (df.tpep_pickup_datetime.dt.weekday < 5)\n",
    "        # Select data\n",
    "        df_period = df[selection]\n",
    "        pairs[label] = make_pair_df(taxi_geo)\n",
    "        count_trips(df_period, pairs[label])\n",
    "        get_mean_distance(df_period, pairs[label])\n",
    "        get_diff(pairs[label], 'income', zone_incomes, min_value=12300)\n",
    "    return pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df = df_taxi[~df_taxi.index.duplicated(keep='first')]\n",
    "#get_trip_datetime(df)\n",
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "time_pairs = get_time_pairs(df, taxi_geo, zone_incomes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 8))\n",
    "for i, (label, pairs) in enumerate(time_pairs.items()):\n",
    "    plt.subplot(2,2,1 + i)\n",
    "    plot_binned(\n",
    "        pairs['income_diff'], pairs['trips'],\n",
    "        #log='y', \n",
    "        alpha=0.8,\n",
    "        xlabel='Difference in income',\n",
    "        ylabel='Number of Trips',\n",
    "        title=label)\n",
    "    #plt.ylim([0.5, 100])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_weekend = df[df.tpep_pickup_datetime.dt.weekday >= 5]\n",
    "df_pairs_weekend = make_pair_df(taxi_geo)\n",
    "count_trips(df_weekend, df_pairs_weekend)\n",
    "get_mean_distance(df_weekend, df_pairs_weekend)\n",
    "get_diff(df_pairs_weekend, 'income', zone_incomes, min_value=12300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 8))\n",
    "\n",
    "for i, miles in enumerate([1, 4, 7]):\n",
    "    df_pairs_miles = filter_zones(df_pairs_weekend, miles, 1.5)\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    plot_binned(\n",
    "        df_pairs_miles['income_diff'],\n",
    "        df_pairs_miles['trips'],\n",
    "        bins=15, #log='y', \n",
    "        alpha=0.8,\n",
    "        xlabel='Difference in income',\n",
    "        ylabel='Number of Trips',\n",
    "        title='{} Miles (weekend)'.format(miles))\n",
    "    #plt.ylim([0.9, 10000])\n",
    "    #plt.xlim([-1, 1])\n",
    "    \n",
    "for i, miles in enumerate([1, 4, 7]):\n",
    "    df_pairs_miles = filter_zones(time_pairs['AM Rush'], miles, 1.5)\n",
    "    plt.subplot(2, 3, i + 4)\n",
    "    plot_binned(\n",
    "        df_pairs_miles['income_diff'],\n",
    "        df_pairs_miles['trips'],\n",
    "        bins=15, #log='y', \n",
    "        alpha=0.8,\n",
    "        xlabel='Difference in income',\n",
    "        ylabel='Number of Trips',\n",
    "        title='{} Miles (rush hour)'.format(miles))\n",
    "    #plt.ylim([0.9, 10000])\n",
    "    #plt.xlim([-1, 1])\n",
    "    \n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Examining individual pairs\n",
    "- To understand these patterns more, let's look at some specific neighborhoods and the trips between them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "def print_zone(geo, location_id):\n",
    "    for r in geo:\n",
    "        if r.attributes['LocationID'] == location_id:\n",
    "            print('{}: {}'.format(\n",
    "                r.attributes['borough'],\n",
    "                r.attributes['zone']))\n",
    "            print('Median income:', r.attributes['income'])\n",
    "            break\n",
    "            \n",
    "def trip_weekday_hist(df, geo, source, dest, ymax):\n",
    "    r_source = taxi_geo_dict[source]\n",
    "    r_dest = taxi_geo_dict[dest]\n",
    "    bins = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]\n",
    "    df[(df.PULocationID == source) & \n",
    "       (df.DOLocationID == dest)\n",
    "      ].tpep_pickup_datetime.dt.weekday.hist(bins=bins)\n",
    "    plt.ylim([0, ymax])\n",
    "    plt.title('{} \\n-> {}'.format(\n",
    "        r_source.attributes['zone'],\n",
    "        r_dest.attributes['zone']))\n",
    "    plt.xticks(\n",
    "        range(7),\n",
    "        ['M', 'Tu', 'W', 'Th', 'F', 'Sa', 'Su'])\n",
    "    plt.xlim([-0.5,6.5])\n",
    "\n",
    "def trip_time_hist(df, geo, source, dest, ymax):\n",
    "    r_source = taxi_geo_dict[source]\n",
    "    r_dest = taxi_geo_dict[dest]\n",
    "    df[(df.PULocationID == source) & \n",
    "       (df.DOLocationID == dest)\n",
    "      ].tpep_pickup_datetime.dt.hour.hist(bins=24)\n",
    "    plt.ylim([0, ymax])\n",
    "    plt.title('{} \\n-> {}'.format(\n",
    "        r_source.attributes['zone'],\n",
    "        r_dest.attributes['zone']))\n",
    "    \n",
    "taxi_geo_dict = {}\n",
    "for r in taxi_geo:\n",
    "    taxi_geo_dict[r.attributes['LocationID']] = r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### neighborhoods with similar income\n",
    "- This code gets a list of neighborhoods with **simiar** incomes that have the most trips between them.\n",
    "- The neighborhood ID numbers appear on the left on the list in parentheses, e.g. `(262, 162)` means all trips between neighborhood with ID 262 and neighborhood with ID 162."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sel = (((time_pairs['AM Rush'].income_diff).abs() < 41339)\n",
    "       & (time_pairs['AM Rush'].trip_distance > 2))\n",
    "\n",
    "time_pairs['AM Rush'][sel].sort_values('trips', ascending=False).head().round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- This code lets you look up a neighborhood's name and income based on its ID number. \n",
    "- Try putting in some of the neighborhood ID numbers from above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_zone(taxi_geo, 148)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_zone(taxi_geo, 162)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- This code gets a list of neighborhoods with very **different** incomes that have the most trips between them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sel = (((time_pairs['AM Rush'].income_diff).abs() > 41339)\n",
    "       & (time_pairs['AM Rush'].trip_distance > 2))\n",
    "\n",
    "time_pairs['AM Rush'][sel].sort_values('trips', ascending=False).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_zone(taxi_geo, 229)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_zone(taxi_geo, 75)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (Just for fun!) Examine trips between two specific neighborhoods by time of day and day of week\n",
    "- Change the neighborhood ID numbers in the cell below to a pair of neighborhoods that you are interested in. \n",
    "- Then run the code below them to see which direction people travel between them, and when."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neighborhood_a = 75 # you can change this value\n",
    "neighborhood_b = 229 # you can change this value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(1,2,1)\n",
    "trip_weekday_hist(df, taxi_geo, source=neighborhood_a, dest=neighborhood_b, ymax=300)  \n",
    "plt.ylabel('trips')\n",
    "plt.subplot(1,2,2)\n",
    "trip_weekday_hist(df, taxi_geo, source=neighborhood_b, dest=neighborhood_a, ymax=300)    \n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(1,2,1)\n",
    "trip_time_hist(df, taxi_geo, source=neighborhood_a, dest=neighborhood_b, ymax=130)\n",
    "plt.ylabel('Trips')\n",
    "plt.xlabel('Hour of the day')\n",
    "plt.subplot(1,2,2)\n",
    "trip_time_hist(df, taxi_geo, source=neighborhood_b, dest=neighborhood_a, ymax=130)\n",
    "plt.xlabel('Hour of the day')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Try It Yourself: race and ethnicity\n",
    "- Income is not the only information we have about neighborhoods. \n",
    "- In the next section, you'll do the same analysis as before, but this time with data about the proportion of a neighborhood that identifies as white and not Hispanic. That is, instead of asking mow many taxi trips connect neighborhoods with different income levels, we will ash how many taxi trips connect neighborhoods with different proprotions of white people and people of color. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\">\n",
    "    \n",
    "#### Short answer 8:\n",
    "- What do you expect the relationship between neighborhoods' racial composition and taxi trips to look like? Explain why."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### First, a map of NYC showing the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zone_race = tract_to_zone(df_census['non_hispanic_white_frac'])\n",
    "for zone in taxi_geo: \n",
    "    zone_id = zone.attributes['LocationID']\n",
    "    zone.attributes['non_hispanic_white_frac'] = zone_race[zone_id]\n",
    "plot_geo_values(taxi_geo, 'non_hispanic_white_frac', \n",
    "                label=\"Fraction Non-Hispanic White\", \n",
    "                cmap='spring')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_diff(df_pairs, 'non_hispanic_white_frac', zone_race)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### trips according to race/ethnicity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 8))\n",
    "plot_binned(\n",
    "    df_pairs['non_hispanic_white_frac_diff'], df_pairs['trips'],\n",
    "    #log='y',\n",
    "    alpha=0.8,\n",
    "    xlabel='Difference in Fraction Non-Hispanic White',\n",
    "    ylabel='Number of Trips')\n",
    "plt.ylim([0,2000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "#### Short answer 9:\n",
    "- Was your prediction correct? Explain why or why not. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Reflection Questions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-warning\">\n",
    "\n",
    "#### Reflection Question 1\n",
    "\n",
    "We showed a number of factors that all influence taxi trips, including income, race, distance, time of day, and day of week. Name one other factor that might influence taxi trips. \n",
    "    \n",
    "- Explain how you would measure that using this data. Then explain how it is (or is not) related to the factors we looked at in the lab.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-warning\">\n",
    "\n",
    "#### Reflection Question 2\n",
    "Pick a city or town that you know that is not NYC. \n",
    "- How might the taxi transport patterns(x,y,z etc.) look different there? Explain why.\n",
    "- In what ways might it look similar? Explain why.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-warning\">\n",
    "\n",
    "#### Reflection Question 3\n",
    "    \n",
    "Taxis are one way to get around a city, but New Yorkers get around many ways including subways, ferries, walking, biking, driving, busses, and even (for the very wealthy) helecopters. \n",
    "    \n",
    "- What is one research question that taxi data would not be very helpful for answering? Why? \n",
    "- What is one research question that taxi data would be useful for answering? Why? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[JJ1961] Jacobs, J. The Death and Life of Great American Cities. Random House, New York."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}