<jupyter_start><jupyter_text>NYC Yellow Cab Efficiency---> The goal of this code is to explore the impact of shared rides on efficiency by using aggregation strategies on NYC Yellow Cab trip data. The deliverables are estimates of gained efficiency if we switch from traditional cab rides to aggressive ride sharing strategies. My approach starts by defining an efficiency metric and an aggregation strategy and then I compare the distribution of efficiency of yellow cabs and ride shares across a range of spatial and temporal gradients in Manhattan.  Contents1. Definining an efficiency metric     * Formal definition.    * Assumptions and simplifications.    * Setup environment.2. Data processing pipeline to calculate efficiency     * Intro to aggregation strategy.    * Get trips in Manhattan.    * Cleanup data and calculate efficiency.3. Aggregation strategy     * Aggregation function.    * Find optimal waiting time threshold. 4. Compare efficiency for modes of transportation in Manhattan     * NYC Yellow cabs vs. Aggregated rides.    * Walking and biking baselines. 5. Visualizations for spatial and temporal trends in efficiency     * Neighborhood network maps.    * Efficiency vs. time of day and time of the week.6. Conclusions     * Take home points.    * Next steps.  Defining an efficiency metric for NYC transport: $$ \textbf{Efficiency} = \frac{\textit{Passenger displacement}}{\textit{cost}} $$$$ \textbf{Passenger displacement} = \mu_{\textit{distance covered}} * \mu_{\textit{number of passengers}} $$$$ \textbf{Cost} = \textit{(fuel cost + social carbon cost) * miles + idle time} + \textit{(drive salary + oppotunity cost) time} $$ Passenger displacement:* Distance travelled by each passenger.* Adds up when multiple passengers share a ride. Cost components:1. Driver salary     * $\$$20 / hour fixed salary.2. Fuel cost    * Assumes 13 mpg fuel efficiency.    * Assumes 25 mph average crusing speed to determine idle time and 0.5 gallon per hour consumption when idle.    * Based on a $\$$2.385 gallon and a $\$$0.453 social carbon cost per gallon emissions.3. Opportunity cost     * $\$$15 / hour opportunity cost per passenger for commuting. Assumptions:* Trips beginning and ending in Manhattan have twice the weight of those that only begin or end in Manhattan when determing efficiency. * Pickup and dropoff walking distances are are set to the Manhattan distance between pairs of pickup and dropoff points.* Shared rides can take up to 6 passengers at once.* Passengers can only aggregate at the beginning of a trip (i.e. no pickups or dropoffs mid-route).* Traffic speed, fuel cost and driver salary are the same for non-aggregated and aggregated rides. Data pipeline:* 1st full week of June 2016 $\Rightarrow$ training (find optimal waiting time).* 2nd full week of June 2016 $\Rightarrow$ validation (compare aggregated vs. yellow cab efficiency).  Load packages and read input data<jupyter_code>import pandas as pd
import numpy as np
from datetime import datetime, timedelta
from functools import partial
from collections import Counter, defaultdict
from multiprocessing import Pool
from sklearn import gaussian_process as gp
import plotnine as p9
import time

# gis imports
import geopandas as gpd
import folium
from shapely.geometry import Point, LineString, Polygon
from shapely.strtree import STRtree
from shapely.ops import nearest_points
from scipy.spatial import cKDTree
from osmnx import quadrat_cut_geometry

# set seed for reproducibility
np.random.seed(4)

# read input data
pd.set_option('display.float_format', lambda x: '%.5f' % x)
boroughs = gpd.read_file('Borough_Boundaries/geo_export_acd45f17-302a-4908-b0d7-223e9510dc04.shp')
trips_all = pd.read_csv('yellow_tripdata_2016-06.csv')
manhattan = boroughs.loc[boroughs.boro_name == 'Manhattan', 'geometry'].values[0]
boroughs['city'] = ['NYC'] * len(boroughs)
all_boroughs = boroughs.dissolve(by='city').geometry.values[0]
trips_all.describe()<jupyter_output><empty_output><jupyter_text>Calculating efficiency ---* Create function to get trips starting or ending in Manhattan.* Create function to get efficiency and emissions for those trips.* Calculate efficiency for Yellow cabs in Manhattan. Function to get Manhattan trips<jupyter_code>def get_trips_polygon(df: pd.DataFrame, polygon: Polygon, bounds: Polygon=None) -> pd.DataFrame:
    """
    Helper function to subset trips to those starting and/or ending inside a georefrenced polygon
    Inputs:
    df: pandas DataFrame with trips
    polygon: shapely polygon
    """
    # create columns for point geometries
    points_pick = [Point(*ele) for ele in zip(df['pickup_longitude'], df['pickup_latitude'])]
    points_drop = [Point(*ele) for ele in zip(df['dropoff_longitude'], df['dropoff_latitude'])]
    pick_id = [id(ele) for ele in points_pick]
    drop_id= [id(ele) for ele in points_drop]
               
    # create a dictionary from id to index
    id_to_idx = {point_id: df.index[idx] for idx, point_id in enumerate(pick_id)}
    for idx, point_id in enumerate(drop_id):
        id_to_idx[point_id] = df.index[idx]
    
    # create R-tree for points to query 
    pick_tree = STRtree(points_pick)
    drop_tree = STRtree(points_drop)
    
    # filter to stay within 4 borough boundaries
    if bounds:
        query_pick = pick_tree.query(bounds)
        quert_drop = drop_tree.query(bounds)
        ids_pick = set([id_to_idx[id(ele)] for ele in query_pick])
        ids_drop = set([id_to_idx[id(ele)] for ele in quert_drop])
        df = df.loc[list(ids_pick.intersection(ids_drop))]
        
        # update trees
        pick_tree = STRtree(query_pick)
        drop_tree = STRtree(quert_drop)

    # chop polygon into quadrats for R-Tree searching
    pol_cut = quadrat_cut_geometry(polygon, 0.025)
    
    # loop through quadrats storing points that start or end in each quadrant
    ids_pick = set()
    ids_drop = set()
    for quadrat in pol_cut:
        drop_quadrat = set([id_to_idx[id(ele)] for ele in drop_tree.query(quadrat)])
        pick_quadrat = set([id_to_idx[id(ele)] for ele in pick_tree.query(quadrat)])
        ids_drop = ids_drop.union(drop_quadrat)
        ids_pick = ids_pick.union(pick_quadrat)
    
    # keep points that start or end in Manhattan, give extra weight to those that do both
    weight2_trips = list(ids_drop.intersection(ids_pick))
    df = df.loc[df.index.isin(list(ids_pick.union(ids_drop)))]
    df = df.assign(weight=[1] * len(df))
    df.loc[weight2_trips, 'weight'] = 2

    # change datetime columns to datetime and return# change time columns to datetime
    for col in ['tpep_pickup_datetime', 'tpep_dropoff_datetime']:
        df[col] = df[col].apply(lambda x: datetime.strptime(x, '%Y-%m-%d %H:%M:%S'))
        
    # return df
    return df<jupyter_output><empty_output><jupyter_text>Get all Manhattan trips during the second week of June 2016<jupyter_code># subset trips by date
trips = trips_all.loc[[12 < int(ele[8:10]) < 20 for ele in trips_all.tpep_pickup_datetime]]

# subset by location and add geometries and weights
tic = time.time()
chunks = (trips.iloc[idx:idx + 10000] for idx in range(0, len(trips), 10000))
pool = Pool(24)
results = pool.map_async(partial(get_trips_polygon, polygon=manhattan,
                                bounds=all_boroughs), chunks)
trips = results.get()
pool.close()
pool.join()
print(f'Finished processing trips in {time.time() - tic} seconds')

# subset trips and explore resulting data
trips = pd.concat(trips)
trips = trips.assign(total_wait=[0] * len(trips)) 
trips.describe()<jupyter_output>Finished processing trips in 18.447881937026978 seconds<jupyter_text>Cleanup trip data* There is no missing data on columns of interest.* Some trips have a passenger count of zero.* Some trips have zero distance.* Some trips have an unrealistic distance (> 500 mi).* Some trips have unrealistic duration given their distance (maybe the driver forgot to turn off the meter?). Data imputation strategy:* Substitute passenger counts of zero with a sample from passenger counts.* Substitute missing / unrealistic distances with Manhattan distances.<jupyter_code># get number of trips with zero distance, zero passengers and unrealisticly high distances
print(f'Number of trips with zero distance: {len(trips.loc[trips.trip_distance == 0])}')
print(f'Number of trips with zero passengers: {len(trips.loc[trips.passenger_count == 0])}')
print(f'Number of trips over 50 miles:  {len(trips.loc[trips.trip_distance > 50])}')

# Helper function top get manhattan distances between vectors of lat, lon
def manhattan_dist(lat1, lat2, lon1, lon2):
    dlat = abs(lat1 - lat2)
    dlon = abs(lon1 - lon2)
    return 69 * dlat + 69 * dlon * np.cos(np.radians(lon2)) 

# function to clean trips
def cleanup_trips(trips: pd.DataFrame, minimum_speed=3.4):
    
    # make sure pickup date comes before dropoff date and duration is reasonable for trip distance
    trips = trips.loc[trips.tpep_pickup_datetime < trips.tpep_dropoff_datetime]
    trips = trips.loc[(trips.trip_distance / ((trips['tpep_dropoff_datetime'] - \
                                              trips['tpep_pickup_datetime']).astype('timedelta64[s]') / 3600)) > minimum_speed]
    
    # fill zero passenger count trips sampling from observed passenger counts
    passenger_counts = list(trips.loc[trips.passenger_count > 0].passenger_count)
    trips.loc[trips.passenger_count == 0, 'passenger_count'] = np.random.choice(passenger_counts, 
                         size=len(trips.loc[trips.passenger_count == 0]))

    # fill trips with unrealistic or zero distance using the Manhattan distance 
    dist_trips = trips.loc[(trips.trip_distance < 50) & (trips.trip_distance != 0)]
    trips.loc[(trips.trip_distance < 50) | (trips.trip_distance != 0), 'trip_distance'] =  \
        manhattan_dist(dist_trips.pickup_latitude, dist_trips.dropoff_latitude, 
                       dist_trips.pickup_longitude, dist_trips.dropoff_longitude)

    return trips

# clean trips and display new values
trips = cleanup_trips(trips)
trips.describe()<jupyter_output>Number of trips with zero distance: 6504
Number of trips with zero passengers: 45
Number of trips over 50 miles:  26<jupyter_text>Efficiency function<jupyter_code># function to get efficiency for trip
def get_efficiency(trips, gallon_price: float=2.385, mpg: float=13.0,
                   social_carbon_cost: float=0.453, idle_cost: float=0.5, 
                   avg_speed: float=25.0, driver_salary: float=20.0, 
                   opportunity_cost: float=15.0) -> pd.DataFrame:
    """
    Determines the efficiency of a trip as passenger displacement / trip cost
    
    gallon_price (float): price of a gallon of fuel in US$
    mpg (float): miles per gallon for vehicles
    idle_cost: gallons consumed per hour when car is idle
    avg_speed: average speed in mph when not suck in traffic
    driver_salary: driver hourly pay in US$
    opportunity_cost: hourly ocpportunity cost for passengers in US$
    
    returns: efficiency for trip, doubled if it starts and ends within Manhattan 
    """
    # add extra weight to manhattan trips
    trips = trips.append(trips.loc[trips.weight == 2], ignore_index=True)
    
    # get displacement
    displacement = trips['trip_distance'] * trips['passenger_count']
    
    # get duration in hours
    duration = (trips['tpep_dropoff_datetime'] - trips['tpep_pickup_datetime']).astype('timedelta64[s]') / 3600
    
    # get fuel and driver cost
    idle_time = np.clip(duration - trips['trip_distance'] / avg_speed, 0, 10E6)
    total_fuel = trips['trip_distance'] / mpg + idle_cost * idle_time
    cost_fuel = total_fuel * (gallon_price + social_carbon_cost)
    emissions = total_fuel * 9.07185 / 1000
    cost_salary = duration * driver_salary + ((trips['total_wait'] + duration)) * opportunity_cost
    
    # get efficiency 
    total_cost = cost_fuel + cost_salary
    efficiency = displacement / total_cost
    
    # return metrics
    metrics = pd.DataFrame({'efficiency': efficiency,
                            'duration': (trips['total_wait'] + duration) * 60,
                            'emissions': emissions})
    return metrics<jupyter_output><empty_output><jupyter_text>Get efficiency for NYC yellow cabs<jupyter_code># get mean efficiency for yellow cabs
get_efficiency(trips).describe()<jupyter_output><empty_output><jupyter_text>Trip aggregation strategy ---1. Merges rides with similar pickup locations and times goint to a similar dropoff location respecting a vehicle capacity of 6 passengers.2. Convert distance of nearby pickups and dropoffs to time assuming a fixed walking speed of 3.6 mph.3. Aggregates passengers from trip $\textit{A}$ to those of trip $\textit{B}$ if doing so does not increase trip duration by more than a given threshold nor exceeds vehicle capacity.4. Waiting time threshold is capped at twice the original trip duration to prevent inneficient aggregations on short trips.5. Choose an optimal waiting time threshold by computing mean efficiency for a random subset of trips across a range of reasonable waiting times.         Get first week of June 2016 to find the most efficient waiting time threshold<jupyter_code># subset training trips by date
trips_train = trips_all.loc[[5 < int(ele[8:10]) < 13 for ele in trips_all.tpep_pickup_datetime]]

# subset by location and add geometries and weights
tic = time.time()
chunks = (trips_train.iloc[idx:idx + 10000] for idx in range(0, len(trips_train), 10000))
pool = Pool(24)
results = pool.map_async(partial(get_trips_polygon, polygon=manhattan,
                                bounds=all_boroughs), chunks)
trips_train = results.get()
pool.close()
pool.join()
trips_train = pd.concat(trips_train)

# cleanup trips
trips_train = cleanup_trips(trips_train)
print(f'Finished processing trips in {time.time() - tic} seconds')<jupyter_output>Finished processing trips in 20.521830320358276 seconds<jupyter_text>Ride aggregation function<jupyter_code># helper function to aggregate by waiting time
def agg_by_time(df: pd.DataFrame, walking_speed: float, max_wait: float,
                max_passengers: int) -> pd.DataFrame:
    
      
    # find distance buffer given walking speed
    dist_buffer = walking_speed * max_wait / 60 
    
    # add rides at capacity to final output
    agg = pd.DataFrame()
    at_capacity = df.passenger_count >= max_passengers
    agg = agg.append(df.loc[at_capacity])
    agg= agg.assign(total_wait=[0] * len(agg))
    df = df.loc[-at_capacity]
    
    # loop through trips
    while len(df) > 0:
        
        # get first row
        curr = df.iloc[0]
        lat_pick, lon_pick = curr.pickup_latitude, curr.pickup_longitude
        lat_drop, lon_drop = curr.dropoff_latitude, curr.dropoff_longitude
        pick_time = curr.tpep_pickup_datetime
        df = df.iloc[1:]
        
        # adjust max wait time for short trips to prevent too much waiting / travelling time
        wait_time = min(max_wait, (pick_time - curr.tpep_dropoff_datetime).seconds / 30)
        
        # filter by pickup time, latitude and dropoff latitude
        time_diffs = (df.loc[:, 'tpep_pickup_datetime'] - pick_time).astype('timedelta64[s]') / 3600
        filtered = df.loc[np.array(time_diffs) < wait_time]
        lat_diffs_pick = abs(filtered.pickup_latitude - lat_pick)
        filtered = filtered.loc[np.array(lat_diffs_pick) < (dist_buffer / 69)]
        lat_diffs_drop = abs(filtered.dropoff_latitude - lat_drop)
        filtered = filtered.loc[np.array(lat_diffs_drop) < (dist_buffer / 69)]
        
        # subset by distance to dropoff location
        distance_drop = manhattan_dist([lat_drop] * len(filtered), filtered.dropoff_latitude,
                                       [lon_drop] * len(filtered), filtered.dropoff_longitude)
        filtered = filtered.assign(wait_time=(distance_drop / walking_speed)) 
        filtered = filtered.loc[filtered.wait_time < wait_time]
        
        # subset by pickup location
        distance_pick = manhattan_dist([lat_pick] * len(filtered), filtered.pickup_latitude,
                                       [lon_pick] * len(filtered), filtered.pickup_longitude)
        filtered['wait_time'] += distance_pick / walking_speed
        
        # threshold by waiting time
        filtered = filtered.loc[filtered.wait_time < wait_time]
         
        # add best passenger groups 
        filtered = filtered.sort_values(by=['wait_time'])
        total_wait = 0
        idcs = set([])
        for idx, row in filtered.iterrows():
            if curr.passenger_count == max_passengers:
                break
            if (row.passenger_count + curr.passenger_count) <= max_passengers:
                total_wait += row.passenger_count * row.wait_time
                curr.loc['passenger_count'] += row.passenger_count
                idcs.add(idx)  
        curr['total_wait'] = total_wait
        agg = agg.append(curr)
        
        # remove aggregated rides from dataframe
        df = df.loc[~df.index.isin(idcs)]
    
    # return aggregated rides
    return agg<jupyter_output><empty_output><jupyter_text>Search for optimal waiting time threshold for aggregation *WARNING: due to quadratic time operations on aggregation function, getting efficiency for all thresholds can take > 25 minutes*<jupyter_code># create date ranges for chunking
start_date = min(trips_train.tpep_pickup_datetime)
end_date = max(trips_train.tpep_pickup_datetime)
total_minutes = ((end_date - start_date).days + 1) * 1440  
date_range = [[start_date + timedelta(minutes=ele), 
               start_date + timedelta(minutes=ele + 15)] for ele in range(0, total_minutes, 15)]
date_range = [date_range[idx] for idx in np.random.choice(range(len(date_range)), size=96, replace=False)]

# search for optimal waiting time 
efficiency_wait = {}
for max_wait in range(5, 26, 5):
    
    # generate chunks using date ranges
    chunks = (trips_train.loc[(trips_train.tpep_pickup_datetime >= date[0]) & 
                        (trips_train.tpep_pickup_datetime < date[1])]  
              for idx, date in enumerate(date_range))

    # process chunks
    tic = time.time()
    pool = Pool(24)
    results = pool.map_async(partial(agg_by_time, walking_speed=3.4, 
                                     max_wait=max_wait, max_passengers=6), chunks)
    agg_trips = results.get()
    pool.close()
    pool.join()
    agg_trips = pd.concat(agg_trips)
    print(f"finished aggregating with a {max_wait} minute threshold in {time.time() - tic} seconds.")
    print(f"number of aggregated trips: {len(agg_trips)}.")
    
    # add weight to Manhattan trips and get efficiency for shared rides 
    agg_trips = agg_trips.append(agg_trips.loc[agg_trips.weight == 2], ignore_index=True)
    efficiency_wait[max_wait] =  get_efficiency(agg_trips, opportunity_cost=15)
    print(f"mean efficiency: {efficiency_wait[max_wait].efficiency.mean()}\n")<jupyter_output>finished aggregating with a 5 minute threshold in 147.716317653656 seconds.
number of aggregated trips: 123240.
mean efficiency: 0.645684499676153

finished aggregating with a 10 minute threshold in 144.93352389335632 seconds.
number of aggregated trips: 112861.
mean efficiency: 0.666883032961444

finished aggregating with a 15 minute threshold in 152.78219413757324 seconds.
number of aggregated trips: 110021.
mean efficiency: 0.6685475643575135

finished aggregating with a 20 minute threshold in 161.4207980632782 seconds.
number of aggregated trips: 108763.
mean efficiency: 0.6676949441097717

finished aggregating with a 25 minute threshold in 170.6648828983307 seconds.
number of aggregated trips: 108147.
mean efficiency: 0.6669731114646525<jupyter_text>Compare modes of transportation in Manhattan ---* Compare using efficiency at the second full week of June 2016.* Use best waiting time threshold for aggregating rides in the previous week.* Obtain passenger mile / cost for each mode of transportation.* Assume walking speed of 3.4 mph and biking speed of 11.5 mph.*WARNING: due to quadratic time operations on aggregation function, getting efficiency for aggregated rides can take > 25 minutes*<jupyter_code># helper function to get walking and biking efficiency
def get_efficiency_walk(trips: pd.DataFrame, speed: float=3.4, opportunity_cost: float=15.0):
    
    # add extra weight to manhattan trips
    trips = trips.append(trips.loc[trips.weight==2], ignore_index=True)
    
    # get walking distances
    distances = manhattan_dist(trips.pickup_latitude, trips.dropoff_latitude, 
                               trips.pickup_longitude, trips.dropoff_longitude)
        
    # get efficiency 
    duration = distances / speed
    efficiency = distances / (duration * opportunity_cost)
    
    # return metrics
    metrics = pd.DataFrame({'efficiency': efficiency,
                            'duration': duration * 60,
                            'emissions': [0] * len(efficiency)})
    return metrics

# aggregate rides using a 15 minute cutoff
start_date = min(trips.tpep_pickup_datetime)
end_date = max(trips.tpep_pickup_datetime)
total_minutes = ((end_date - start_date).days + 1) * 1440  
date_range = [[start_date + timedelta(minutes=ele), 
               start_date + timedelta(minutes=ele + 15)] for ele in range(0, total_minutes, 15)]
chunks = (trips.loc[(trips.tpep_pickup_datetime >= date[0]) & 
                    (trips.tpep_pickup_datetime < date[1])]  
          for idx, date in enumerate(date_range))

pool = Pool(24)
results = pool.map_async(partial(agg_by_time, walking_speed=3.4, 
                                 max_wait=15, max_passengers=6), chunks)
agg_trips = results.get()
pool.close()
pool.join()
agg_trips = pd.concat(agg_trips)

# walking
walk_eff = get_efficiency_walk(trips)
print(f"mean efficiency for walking: {walk_eff.efficiency.mean()}\n")

# biking
bike_eff = get_efficiency_walk(trips, speed=11.5)
print(f"mean efficiency for biking: {bike_eff.efficiency.mean()}\n")

# yellow cab
cab_eff = get_efficiency(trips)
print(f"mean efficiency for yellow cabs: {cab_eff.efficiency.mean()}\n")

# ride share 
share_eff = get_efficiency(agg_trips)
print(f"mean efficiency for aggregated rides: {share_eff.efficiency.mean()}\n")<jupyter_output>mean efficiency for walking: 0.22666666666666657

mean efficiency for biking: 0.7666666666666668

mean efficiency for yellow cabs: 0.35677535293897467

mean efficiency for aggregated rides: 0.6918108208367088<jupyter_text>Create Visualizations ---* Neighborhood traffic map.* Efficiency by means of transportation.* Split by neighborhood: Downtown / Midtown / Uptown.* Split by time of day: Morning / Afternoon / Evening.* Split by weekdays vs. weekends. Get traffic between neighborhoods<jupyter_code># function to find traffic between neighborhoods
def get_traffic(trips: pd.DataFrame, neighborhoods: gpd.GeoDataFrame):
    
    # extract locations for neighborhoods, pickups and dropoffs
    pickup_points = list(zip(trips.pickup_latitude, trips.pickup_longitude))
    dropoff_points = list(zip(trips.dropoff_latitude, trips.dropoff_longitude))
    nb_points = list(zip(neighborhoods.geometry.y, neighborhoods.geometry.x))
    
    # get number of trips between each pair neighborhoods removing inter neighborhood trips
    nb_tree = cKDTree(nb_points)
    _, idx_pick = nb_tree.query(pickup_points, k=1)
    _, idx_drop = nb_tree.query(dropoff_points, k=1)
    pairs = [tuple(sorted(ele)) for ele in (zip(neighborhoods.loc[idx_pick, 'name'].values,
                                                neighborhoods.loc[idx_drop, 'name'].values))]
    trips = trips.assign(nb=pairs)
    traffic = Counter([ele for ele in pairs if ele[0] != ele[1]])
    
    # get the 200 most frequent trips
    traffic = traffic.most_common(200)
    
    # get locations and efficiency for each pair of neighborhoods
    traffic_df = pd.DataFrame()
    for ele, val in traffic:
        coords = [neighborhoods.loc[neighborhoods.name==ele[0], 'geometry'].y.values[0],
                  neighborhoods.loc[neighborhoods.name==ele[0], 'geometry'].x.values[0],
                  neighborhoods.loc[neighborhoods.name==ele[1], 'geometry'].y.values[0],
                  neighborhoods.loc[neighborhoods.name==ele[1], 'geometry'].x.values[0]]
        row = pd.Series({'n_trips': val, 
                         'efficiency': trips.loc[trips.nb == ele, 'efficiency'].mean(),
                         'emissions': trips.loc[trips.nb == ele, 'emissions'].sum(),
                         'passengers': trips.loc[trips.nb == ele, 'passenger_count'].mean(),
                         'coordinates': coords},
                       name=ele)
        traffic_df = traffic_df.append(row)
    return traffic_df
    

# get neighborhood locations and add efficiency metrics to trips
neighborhoods = gpd.read_file('Neighborhood_names/geo_export_ecd7b650-b6be-4ee3-827b-92a33f23d30f.shp')
trips = pd.concat([trips.reset_index(), get_efficiency(trips)], axis=1, join='inner')
agg_trips = pd.concat([agg_trips.reset_index(), get_efficiency(agg_trips)], axis=1, join='inner')

# get traffic 
traffic = get_traffic(trips, neighborhoods)
traffic_agg = get_traffic(agg_trips, neighborhoods)<jupyter_output><empty_output><jupyter_text>Traffic maps* Line width is proportional to the total number of trips along that route.* Hovering over lines provides tooltips with mean efficiency, number of trips, mean number of passengers and total emissions for a route.* Color scales from red to green with efficiency:    * efficiency $<$ 0.2            $\Rightarrow$ red    * 0.2 $\leq$ efficiency $<$ 0.4 $\Rightarrow$ orange    * 0.4 $\leq$ efficiency $<$ 0.6 $\Rightarrow$ blue    * efficiency $\geq$ 0.6 $\Rightarrow$ green Yellow cab traffic map<jupyter_code># start yellow cab map
cab_map = folium.Map(location=[40.746911, -73.9716688], zoom_start=13, tiles='Stamen Toner')

# helper function to color by efficiency
def get_color(num):
    if num < 0.2:
        return 'red'
    elif num < 0.4:
        return 'orange'
    elif num < 0.6:
        return 'blue'
    else:
        return 'green'

# add routes including tooltips for efficiency, emissions and number of trips
for idx, route in traffic.iterrows():
    coords = list(route.coordinates)
    points = [tuple(coords[0:2]), tuple(coords[2:])]
    tooltip = f"{idx[0]} <--> {idx[1]} <br> mean efficiency: {round(route.efficiency, 3)} miles / $<br>" + \
              f"total emissions: {round(route.emissions, 3)} tons CO<sup>2</sup> <br> total trips: {int(route.n_trips)}" + \
              f"<br> mean passengers: {round(route.passengers, 3)}"
    folium.PolyLine(points, color=get_color(route.efficiency), weight=route.n_trips / 2500, 
                    opacity=0.5, tooltip=tooltip).add_to(cab_map)

cab_map<jupyter_output><empty_output><jupyter_text>Shared ride traffic map* line width scaling is 4x steeper to deal with the dramatically smaller number of trips.<jupyter_code># start shared ride map
share_map = folium.Map(location=[40.746911, -73.9716688], zoom_start=13, tiles='Stamen Toner')

# add routes
for idx, route in traffic_agg.iterrows():
    coords = list(route.coordinates)
    points = [tuple(coords[0:2]), tuple(coords[2:])]
    tooltip = f"{idx[0]} <--> {idx[1]} <br> mean efficiency: {round(route.efficiency, 3)} miles / $<br>" + \
              f"total emissions: {round(route.emissions, 3)} tons CO<sup>2</sup> <br> total trips: {int(route.n_trips)}" + \
              f"<br> mean passengers: {round(route.passengers, 3)}"
    folium.PolyLine(points, color=get_color(route.efficiency), weight=route.n_trips / 650, 
                    opacity=0.5, tooltip=tooltip).add_to(share_map)

share_map<jupyter_output><empty_output><jupyter_text>Plot efficiency across time of the day * Efficiency vs. hour.* Weekdays vs. weekends.* Gaussian Process Regression (GPR) curves with Bayesian posterior distributions for mean efficiency by time.     * Line: predicted mean.    * Swat: 95% credible interval for prediction. * Warning: the GPR implementation I am using is slow even with a modest sample size, expect ~ 5 minutes*<jupyter_code># add time columns to trips
def add_time_cols(trips: pd.DataFrame):
    
    # add time of the day
    trips = trips.assign(hour=trips.tpep_pickup_datetime.dt.hour)
    
    # add weekday vs. weekend
    trips = trips.assign(weekday=['weekday'] * len(trips))
    trips.loc[trips.tpep_pickup_datetime.dt.weekday > 4, 'weekday'] = 'weekend'
    
    return trips

# format trips for plotting
n_samples = 5000
trips_plot = pd.concat([add_time_cols(trips.sample(n_samples, random_state=4)),
                        add_time_cols(agg_trips.sample(n_samples, random_state=4))])
trips_plot = trips_plot.assign(trip_type = ['Yellow Cab'] * n_samples + ['Aggregated'] * n_samples)

# define gaussian kernel for GPR
kernel = gp.kernels.ConstantKernel(1.0, (1e-1, 20)) * gp.kernels.RBF(10.0, (1e-3, 5))

# plot efficiency for NYC yellow cabs and aggregated rides across the day
eff_plot = (p9.ggplot(data=trips_plot, mapping=p9.aes(x='hour', y='efficiency', 
                                                     fill='trip_type')) +
           p9.theme_light() +
           p9.geom_smooth(alpha=0.5, method='gpr', method_args={'kernel': kernel,
                                                                'n_restarts_optimizer': 10,
                                                                'alpha': 0.25,
                                                                'normalize_y': True}) +
           p9.scale_fill_manual(values=["green", "#FF4500"]) +
           p9.facet_wrap('weekday') +
           p9.labels.labs(x='Time of day (hours)', y='Efficiency (passenger mile / $)') +
           p9.ggtitle('Posterior for Efficiency'))
eff_plot<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09648268216481339, 'fraction_numerical': 0.02664023919105049, 'mean_word_length': 3.5398810431599816, 'pattern_counts': {'":': 0, '<': 87, '<?xml version=': 0, '>': 81, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '212730', 'license': 'MIT', 'n_tokens_mistral': 9604, 'n_tokens_neox': 8796, 'n_words': 2659, 'path': 'NYC_YellowCab/Trip_Aggregation.ipynb', 'repo_name': 'bentocg/DataChallenges', 'size': '29767'}