Code
import pandas as pd
import numpy as np
import glob
import geopandas as gpd1: Setup and Initial Data Processing
This notebook builds the core dataset used throughout the project. Raw data is cleaned, merged, and saved as a single processed dataset for downstream analysis and modeling
Data Sources
- CitiBike trip Data: I choose the months of April, May, and June of 2024, in hopes of getting a nice variance of weather values, as well as potentially fluctuating ridership counts from school being in vs out of session for NYC students
- NOAA Weather: Hourly temperature, precipitation, and wind speed from the Central Park Weather Station (WBAN 94728)
- NYC Spatial Data: Neighborhood Tabulation Area (NTA) boundaries for mapping, and school locations from the NYC Department of City Planning
1.1: Load and Aggregate CitiBike Data
Load all CSV files matching file_pattern from CitiBike data and aggregate it to the hourly station level. Track departures, arrivals, bike type breakdown, and trip duration. This will serve as the primary source of downstream analysis, and data additions will build on top of this
I chose 2 CSV’s from each month, which should represent a diverse enough range of ridership to be representative for this project. You can call the load_and_aggregate_month function on whatever months you choose to include, as long as you have the corresponding weather data, as it not matching will cause downstream merging issues. You will also have to adjust the bleed row cutoff for the first month in 1.1b
Code
def load_and_aggregate_month(file_pattern):
cols = ['rideable_type', 'started_at', 'ended_at', 'start_station_name', 'start_station_id', 'end_station_name',
'end_station_id', 'start_lat', 'start_lng', 'end_lat', 'end_lng']
files = glob.glob(file_pattern)
departure_chunks = []
arrival_chunks = []
for f in files:
# create chunk from csv file, add new features
chunk = pd.read_csv(f, usecols=cols, low_memory=False)
chunk['started_at'] = pd.to_datetime(chunk['started_at'])
chunk['ended_at'] = pd.to_datetime(chunk['ended_at'])
chunk['start_hour'] = chunk['started_at'].dt.floor('h')
chunk['end_hour'] = chunk['ended_at'].dt.floor('h')
# clean any null entries just in case
chunk['start_lat'] = chunk['start_lat'].replace(0, pd.NA)
chunk['start_lng'] = chunk['start_lng'].replace(0, pd.NA)
chunk['end_lat'] = chunk['end_lat'].replace(0, pd.NA)
chunk['end_lng'] = chunk['end_lng'].replace(0, pd.NA)
chunk = chunk.dropna(subset=['start_station_id', 'start_lat', 'start_lng', 'end_station_id'])
# compute trip duration, clear any outlier / invalid durations
chunk['duration_min'] = (chunk['ended_at'] - chunk['started_at']).dt.total_seconds() / 60
chunk = chunk[(chunk['duration_min'] > 0) & (chunk['duration_min'] < 1440)]
# aggregate departure chunk
departures = (chunk.groupby(['start_hour', 'start_station_id', 'start_station_name', 'start_lat', 'start_lng', 'rideable_type'])
.agg(
departures=('started_at', 'count'),
avg_duration_min=('duration_min', 'mean')
)
.reset_index())
departure_chunks.append(departures)
# aggregate arrival chunk
arrivals = (chunk.groupby(['end_hour', 'end_station_id'])
.agg(arrivals=('started_at', 'count'))
.reset_index()
.rename(columns={
'end_hour': 'start_hour',
'end_station_id': 'start_station_id'
}))
arrival_chunks.append(arrivals)
# merge departures, arrivals across files
departure_aggregates = (pd.concat(departure_chunks)
.groupby(['start_hour', 'start_station_id', 'start_station_name',
'start_lat', 'start_lng', 'rideable_type'])
.agg(
departures=('departures', 'sum'),
avg_duration_min=('avg_duration_min', 'mean')
)
.reset_index())
arrival_aggregates = (pd.concat(arrival_chunks)
.groupby(['start_hour', 'start_station_id'])['arrivals']
.sum()
.reset_index())
print("departure columns:", departure_aggregates.columns.tolist())
print("arrival columns:", arrival_aggregates.columns.tolist())
# merge arrivals and departures, add total_traffic feature
month_agg = departure_aggregates.merge(arrival_aggregates, on=['start_hour', 'start_station_id'], how='left')
month_agg['arrivals'] = month_agg['arrivals'].fillna(0).astype(int)
month_agg['total_traffic'] = month_agg['departures'] + month_agg['arrivals']
return month_agg
# this is the place if you want to change which months you choose. aggregates using all files
# that match the given month header
april = load_and_aggregate_month('../data/raw/citibike/202404*.csv')
may = load_and_aggregate_month('../data/raw/citibike/202405*.csv')
june = load_and_aggregate_month('../data/raw/citibike/202406*.csv')
df = pd.concat([april, may, june], ignore_index=True)
print(df.head())departure columns: ['start_hour', 'start_station_id', 'start_station_name', 'start_lat', 'start_lng', 'rideable_type', 'departures', 'avg_duration_min']
arrival columns: ['start_hour', 'start_station_id', 'arrivals']
departure columns: ['start_hour', 'start_station_id', 'start_station_name', 'start_lat', 'start_lng', 'rideable_type', 'departures', 'avg_duration_min']
arrival columns: ['start_hour', 'start_station_id', 'arrivals']
departure columns: ['start_hour', 'start_station_id', 'start_station_name', 'start_lat', 'start_lng', 'rideable_type', 'departures', 'avg_duration_min']
arrival columns: ['start_hour', 'start_station_id', 'arrivals']
start_hour start_station_id start_station_name \
0 2024-03-31 11:00:00 5721.07 Thompson St & Bleecker St
1 2024-03-31 11:00:00 6224.05 W 20 St & 8 Ave
2 2024-03-31 12:00:00 5335.03 Bialystoker Pl & Delancey St
3 2024-03-31 13:00:00 6224.05 W 20 St & 8 Ave
4 2024-03-31 13:00:00 6224.06 8 Ave & W 24 St
start_lat start_lng rideable_type departures avg_duration_min arrivals \
0 40.728401 -73.999688 classic_bike 1 756.402067 0
1 40.743453 -74.000040 classic_bike 1 877.048350 0
2 40.716226 -73.982612 classic_bike 1 827.761483 0
3 40.743453 -74.000040 classic_bike 1 1128.324350 0
4 40.745911 -73.998071 classic_bike 1 743.863017 0
total_traffic
0 1
1 1
2 1
3 1
4 1 1.1b: Collapse Rideable Type for Modeling
Create a collapsed version for modeling where Y = total departures per station per hour. Leave a separate table intact including splits by bike type (electric vs classic), although that is not the primary dataset that will be used
Code
df_model = (df.groupby(['start_hour', 'start_station_id', 'start_station_name', 'start_lat', 'start_lng'])
.agg(
ride_count=('departures', 'sum'),
arrivals=('arrivals', 'sum'),
total_traffic=('total_traffic', 'sum'),
avg_duration_min=('avg_duration_min', 'mean')
)
.reset_index())
# drop march bleed rows
df = df[df['start_hour'] >= '2024-04-01'].copy()
df_model = df_model[df_model['start_hour'] >= '2024-04-01'].copy()
print(df_model.shape)(1892070, 9)1.2: Spatial Joining
Take the latitude and longitude points of stations, and map them to a neighborhood via an NTA shapefile. Also join with school data. Stations that are on water boundaries are assigned via a nearest-neighbor fallback. Some feature engineering on the school data is also done ahead of time.
Code
# form dataframe of unique stations
stations = df_model[['start_station_id', 'start_station_name', 'start_lat', 'start_lng']].drop_duplicates()
# initialize a GeoDataFrame with coordinates from each station
gdf = gpd.GeoDataFrame(stations,
geometry=gpd.points_from_xy(stations.start_lng, stations.start_lat),
crs='EPSG:4326')
# load neighborhood tabulation areas
nta = gpd.read_file('../data/raw/spatial/nta.shp').to_crs('EPSG:4326')
# left join the GeoDataFrame and the nta shape file to map each station to an NTA polygon
stations_nta = gpd.sjoin(gdf, nta[['ntaname', 'boroname', 'geometry']], how='left', predicate='within')
# verify any stations that may have not been matched to a polygon (e.g. ones that exist on a boundary)
unmatched = stations_nta['ntaname'].isnull().sum()
print(f"{unmatched} unmatched stations out of {len(stations_nta)}")
# join unmatched stations to nearest neighbor, merge back with original stations_nta dataframe
if unmatched > 0:
matched = stations_nta[stations_nta['ntaname'].notna()]
unmatched_gdf = gdf[~gdf['start_station_id'].isin(matched['start_station_id'])]
nearest = gpd.sjoin_nearest(unmatched_gdf.to_crs('EPSG:3857'),
nta[['ntaname', 'boroname', 'geometry']].to_crs('EPSG:3857'),
how='left').to_crs('EPSG:4326')
stations_nta = pd.concat([matched, nearest], ignore_index=True)
# load school location data and assign it to a station
schools = gpd.read_file('../data/raw/spatial/schools.shp')
print(schools.columns.tolist())
print(schools.crs)
print(schools.shape)
stations_proj = gdf.to_crs('EPSG:3857')
stations_school = gpd.sjoin_nearest(
stations_proj,
schools[['geometry']],
how='left',
distance_col='dist_to_nearest_school_m'
)
# drop duplicates to maintain one row per stations
stations_school = (stations_school[['start_station_id','dist_to_nearest_school_m']]
.drop_duplicates(subset='start_station_id'))
# binary feature to denote whether a station is within 300 meters of a school - makes for interesting insight down the line
stations_school['near_school'] = (
stations_school['dist_to_nearest_school_m'] <= 300
).astype(int)
# merge back
df_model = df_model.merge(
stations_nta[['start_station_id', 'ntaname', 'boroname']],
on='start_station_id', how='left')
df = df.merge(
stations_nta[['start_station_id', 'ntaname', 'boroname']],
on='start_station_id', how='left')36 unmatched stations out of 2289
['ATS', 'Building_C', 'Location_C', 'Name', 'Geographic', 'Latitude', 'Longitude', 'geometry']
EPSG:3857
(1950, 8)1.3: Weather Merge
Merge CitiBike data with NOAA hourly readings by timestamp
Code
# initialize weather dataframe from NOAA local climatological data for April-June 2024
weather = pd.read_csv('../data/raw/weather/noaa_nyc_lcd_aprmayjun_2024.csv')
weather['DATE'] = pd.to_datetime(weather['DATE']).dt.floor('h')
weather = weather[['DATE', 'HourlyDryBulbTemperature', 'HourlyPrecipitation', 'HourlyWindSpeed']]
# replace any trace precipitation (labeled 'T' in NOAA data) with a small nonzero value
weather['HourlyPrecipitation'] = weather['HourlyPrecipitation'].replace('T', 0.001)
# convert all weather data entries to numeric values
weather_cols = ['HourlyDryBulbTemperature', 'HourlyPrecipitation', 'HourlyWindSpeed']
for col in weather_cols:
weather[col] = pd.to_numeric(weather[col], errors="coerce")
# drop duplicates from weather data, as the same readings may exist, but under different report types (NOAA data quirk)
weather = weather.dropna(subset=['DATE']).drop_duplicates(subset=['DATE'])
# merge into model dataframe, fill gaps in weather data
df_model = df_model.merge(weather, left_on='start_hour', right_on='DATE', how='left')
for col in weather_cols:
df_model[col] = df_model[col].ffill(limit=3)
df_model['HourlyPrecipitation'] = df_model['HourlyPrecipitation'].fillna(0)
df_model['HourlyWindSpeed'] = df_model['HourlyWindSpeed'].fillna(df_model['HourlyWindSpeed'].median())1.4: Additional Time Features & Final Checks
Additional time features derived for downstream analysis. Analyze the final dataframes that will be saved
Code
# add additional time features to the data
df_model['hour_of_day'] = df_model['start_hour'].dt.hour
df_model['day_of_week'] = df_model['start_hour'].dt.dayofweek
df_model['month'] = df_model['start_hour'].dt.month
df_model['is_weekend'] = df_model['day_of_week'].isin([5,6]).astype(int)
df_model['is_rush_hour'] = df_model['hour_of_day'].isin([7,8,9,16,17,18,19]).astype(int)
print("\n── df_model (used for modeling + most EDA) ──")
print(df_model.shape)
print(df_model.dtypes)
print(df_model.isnull().sum())
print(df_model.describe())
print("\n── df (rideable_type split, used for bike type EDA) ──")
print(df.shape)
print(df.head())
── df_model (used for modeling + most EDA) ──
(1893958, 20)
start_hour datetime64[us]
start_station_id str
start_station_name str
start_lat float64
start_lng float64
ride_count int64
arrivals int64
total_traffic int64
avg_duration_min float64
ntaname str
boroname str
DATE datetime64[us]
HourlyDryBulbTemperature float64
HourlyPrecipitation float64
HourlyWindSpeed float64
hour_of_day int32
day_of_week int32
month int32
is_weekend int64
is_rush_hour int64
dtype: object
start_hour 0
start_station_id 0
start_station_name 0
start_lat 0
start_lng 0
ride_count 0
arrivals 0
total_traffic 0
avg_duration_min 0
ntaname 0
boroname 0
DATE 1502
HourlyDryBulbTemperature 3490
HourlyPrecipitation 0
HourlyWindSpeed 0
hour_of_day 0
day_of_week 0
month 0
is_weekend 0
is_rush_hour 0
dtype: int64
start_hour start_lat start_lng ride_count \
count 1893958 1.893958e+06 1.893958e+06 1.893958e+06
mean 2024-05-15 21:55:14.187115 4.073949e+01 -7.395968e+01 3.048675e+00
min 2024-04-01 00:00:00 4.063338e+01 -7.407106e+01 1.000000e+00
25% 2024-04-23 12:00:00 4.070167e+01 -7.398752e+01 1.000000e+00
50% 2024-05-15 14:00:00 4.073544e+01 -7.396316e+01 2.000000e+00
75% 2024-06-07 09:00:00 4.076767e+01 -7.393641e+01 4.000000e+00
max 2024-06-30 23:00:00 4.088630e+01 -7.384672e+01 8.300000e+01
std NaN 5.078752e-02 3.411769e-02 3.297049e+00
arrivals total_traffic avg_duration_min \
count 1.893958e+06 1.893958e+06 1.893958e+06
mean 4.296156e+00 7.344831e+00 1.325942e+01
min 0.000000e+00 1.000000e+00 1.002100e+00
25% 0.000000e+00 1.000000e+00 6.529794e+00
50% 0.000000e+00 3.000000e+00 1.031778e+01
75% 4.000000e+00 8.000000e+00 1.578089e+01
max 2.500000e+02 3.090000e+02 1.437528e+03
std 9.606297e+00 1.174319e+01 1.747202e+01
DATE HourlyDryBulbTemperature \
count 1892456 1.890468e+06
mean 2024-05-15 22:11:28.278935 6.592976e+01
min 2024-04-01 00:00:00 3.700000e+01
25% 2024-04-23 12:00:00 5.700000e+01
50% 2024-05-15 15:00:00 6.700000e+01
75% 2024-06-07 10:00:00 7.500000e+01
max 2024-06-30 23:00:00 9.300000e+01
std NaN 1.137510e+01
HourlyPrecipitation HourlyWindSpeed hour_of_day day_of_week \
count 1.893958e+06 1.893958e+06 1.893958e+06 1.893958e+06
mean 1.823546e-03 4.790189e+00 1.351341e+01 2.973247e+00
min 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
25% 0.000000e+00 3.000000e+00 9.000000e+00 1.000000e+00
50% 0.000000e+00 5.000000e+00 1.400000e+01 3.000000e+00
75% 0.000000e+00 7.000000e+00 1.800000e+01 5.000000e+00
max 5.600000e-01 2.900000e+01 2.300000e+01 6.000000e+00
std 1.730785e-02 3.343479e+00 5.758814e+00 1.979552e+00
month is_weekend is_rush_hour
count 1.893958e+06 1.893958e+06 1.893958e+06
mean 4.992520e+00 2.757606e-01 3.929227e-01
min 4.000000e+00 0.000000e+00 0.000000e+00
25% 4.000000e+00 0.000000e+00 0.000000e+00
50% 5.000000e+00 0.000000e+00 0.000000e+00
75% 6.000000e+00 1.000000e+00 1.000000e+00
max 6.000000e+00 1.000000e+00 1.000000e+00
std 8.095207e-01 4.468969e-01 4.884000e-01
── df (rideable_type split, used for bike type EDA) ──
(2580089, 12)
start_hour start_station_id start_station_name start_lat \
0 2024-04-01 2883.03 3 Ave & Wakeman Pl 40.638246
1 2024-04-01 3046.03 E 3 St & Church Ave 40.643817
2 2024-04-01 3056.05 12 Ave & 36 St 40.643546
3 2024-04-01 3169.07 53 St & 4 Ave 40.644862
4 2024-04-01 3172.01 Tilden Ave & Nostrand Ave 40.646830
start_lng rideable_type departures avg_duration_min arrivals \
0 -74.024714 electric_bike 2 8.777042 0
1 -73.977433 electric_bike 1 5.181850 0
2 -73.986418 electric_bike 2 17.163958 0
3 -74.014531 electric_bike 1 2.992950 0
4 -73.949330 electric_bike 1 14.203033 0
total_traffic ntaname boroname
0 2 Sunset Park (West) Brooklyn
1 1 Kensington Brooklyn
2 2 Kensington Brooklyn
3 1 Sunset Park (West) Brooklyn
4 1 East Flatbush-Erasmus Brooklyn 1.5: Data Saving
Save all the relevant dataframes to csv files. Outputs should be
data/processed/hourly_merged.csv- one row per station per hourdata/processed/station_lookup.csv- each individual station with spatial attributesdata/processed/hourly_by_rideable_type.csv- same as hourly_merged, split by classic vs electric bike type
Code
# save main working dataset
df_model.to_csv('../data/processed/hourly_merged.csv', index=False)
# save rideable type split for classic vs electric bike explorations
df.to_csv('../data/processed/hourly_by_rideable_type.csv', index=False)
# save station/geo-spatial data
station_lookup = (df_model[['start_station_id','start_station_name',
'start_lat','start_lng','ntaname','boroname']]
.drop_duplicates()
.merge(stations_school[['start_station_id', 'dist_to_nearest_school_m', 'near_school']],
on='start_station_id', how='left'))
station_lookup.to_csv('../data/processed/station_lookup.csv', index=False)