Bus Mapping and Basic GIS with GeoPandas¶

# !pip install geopandas 
# !pip install rtree

# an rtree dependency 
# !brew install Spatialindex ( For the mac )

import pandas as pd 
import numpy as np

# geospatial tools
import geopandas
import rtree

# plotting tools
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline

# Source: DC open data portal http://opendata.dc.gov/datasets files about WMTA bus lines and stops 

# http://opendata.dc.gov/datasets/metro-bus-stops
bus_stops = geopandas.read_file('data/Metro_Bus_Stops/')

# http://opendata.dc.gov/datasets/metro-bus-lines
bus_routes = geopandas.read_file('data/Metro_Bus_Lines/')

bus_stops.crs

bus_routes.crs

# Source: US Census ACS 2015 5-year estimates

acs = geopandas.read_file('data/ACS_2015_5YR_BG_11_DISTRICT_OF_COLUMBIA.gdb')

acs.crs

acs = acs.to_crs({'init': 'epsg:4326'})
acs.crs

bus_stops.geometry[0:4]

# The Geopandas dataframe is aware of the POINT object within the column
# and knows how to plot it

bus_stops.geometry[0:4].plot()

pd.DataFrame(bus_stops.geometry[0:4]).geometry

# This won't work, this will throw TypeError
# pd.DataFrame(bus_stops.geometry[0:4]).plot()

bus_routes.geometry.head()

bus_routes.sample(2).plot(color='gray')

acs.geometry.head()

acs.plot(color = 'red', alpha=0.3)

# Create the plotting area 
fig, ax = plt.subplots(1, figsize=(7,7))

# Plot bus routes
base = bus_routes.plot(ax=ax, color='gray')

# Plot bus stops
bus_stops.plot(ax=base, marker="o", alpha=0.5)

# Plot bus stops
acs.plot(ax=base, alpha=0.2, color = 'red')
_ = ax.axis('off')

crosswalk_stops = geopandas.tools.sjoin(
    geopandas.GeoDataFrame(bus_stops[['OBJECTID','geometry']]),
    geopandas.GeoDataFrame(acs[['GEOID','geometry']]), 
                                        # Operation we are using is 'within'
                           how="inner", op='within')

crosswalk_stops.head()

bus_stops = bus_stops.merge(crosswalk_stops[['OBJECTID','GEOID']], on='OBJECTID', how='left')

sample_bg = acs.GEOID.sample().tolist()[0]

fig, ax = plt.subplots(1, figsize=(10,10))

# Plot bus stops within the block group
base = geopandas.GeoDataFrame(bus_stops[bus_stops.GEOID == sample_bg ]).\
            plot(ax=ax ,alpha=0.7)

# Plot the block group polygon
acs[acs.GEOID == sample_bg].plot(ax=base, alpha=0.5, color = 'gray')

#fig.suptitle(acs[acs.GEOID == sample_bg].NAMELSAD.tolist()[0])
x, y = acs[acs.GEOID == sample_bg][['TRACTCE','NAMELSAD']].to_dict('records')[0].values()
fig.suptitle("Tract: {0}\n{1}".format(x,y))

_ = ax.axis('off')

bus_stops['geom2']= bus_stops.geometry.map(lambda x: x.buffer(0.0005))
# For now lets not think too hard about the size of the buffer, since we would need to delve into the coordinate system

geopandas.GeoSeries(bus_stops[0:2].geom2).plot(color ='black')

geopandas.GeoSeries(bus_stops.geom2).plot()

# Plot only the stops that have a GEOID ( ACS block-group in DC )
geopandas.GeoSeries(bus_stops[~bus_stops.GEOID.isnull()].geom2).plot()

# Start by looking at one bus route
sample_route = bus_routes.ROUTE.sample(1).tolist()[0]
oneroute = bus_routes[bus_routes.ROUTE == sample_route].sample(1)
oneroute

# Perform the join just using the one bus route

sample_crosswalk = geopandas.tools.sjoin(
    geopandas.GeoDataFrame(oneroute[['OBJECTID','geometry']]),
    geopandas.GeoDataFrame(bus_stops[['OBJECTID','geom2']].rename(columns={'geom2': 'geometry'})), 
                           how="inner", op='intersects', 
                           lsuffix='_route',
                        rsuffix = '_stop')

# Join the stops and routes using the crosswalk

stops_on_route = sample_crosswalk[['OBJECTID__route','OBJECTID__stop']].\
                        merge(oneroute[['OBJECTID']], 
                                 left_on = 'OBJECTID__route',
                                 right_on = 'OBJECTID',
                                 how = 'inner')


fig, ax = plt.subplots(1, figsize=(11,8))

# Plot bus stops along the route      
base = geopandas.GeoDataFrame(
            bus_stops.merge(stops_on_route,
                         left_on = 'OBJECTID',
                         right_on = 'OBJECTID__stop',
                         how = 'inner'),
            geometry = 'geom2').\
            plot(ax=ax ,alpha=0.7)    
    

# Plot the single bus route
oneroute.plot(ax=base, alpha=0.5, color = 'green')

name = oneroute.LINENAME.tolist()[0]
fig.suptitle("Route Name: {0}".format(name))

_ = ax.axis('off')

crosswalk_acs = geopandas.tools.sjoin(
    geopandas.GeoDataFrame(oneroute[['OBJECTID','geometry']]),
    geopandas.GeoDataFrame(acs[['GEOID','geometry']]), 
                           how="inner", op='intersects')

crosswalk_acs.head()

fig, ax = plt.subplots(1, figsize=(10,10))


# The block groups that intersect that bus line
base = acs.merge(crosswalk_acs[['GEOID']], how = 'inner').geometry.\
    plot(ax = ax,  color='gray', alpha=0.3)
    
# The one bus line
oneroute.\
    plot(ax = base, color = 'green')
    
# Bus stops within the ACS block groups that intersect that bus line
geopandas.GeoSeries(
    bus_stops.merge(crosswalk_acs[['GEOID']], how = 'inner').geom2).\
    plot(ax = base,alpha=0.5)

name = oneroute.LINENAME.tolist()[0]
fig.suptitle("Route Name: {0}".format(name))
_ = ax.axis('off')

dc_shape = acs.unary_union
dc_shape

lines_in_dc = bus_routes.intersection(dc_shape)
lines_in_dc.plot()

# The difference shows  bus route boundaries that are outside of DC
line_diffs = bus_routes.difference(dc_shape)
line_diffs.plot()

# Percentage of the bus line mileage within DC boundaries.
100*round(lines_in_dc.length.sum()/bus_routes.geometry.length.sum(),2)

crosswalk_route_stop = geopandas.tools.sjoin(
    geopandas.GeoDataFrame(bus_routes[['OBJECTID','geometry']]),
    geopandas.GeoDataFrame(bus_stops[['OBJECTID','geom2']].rename(columns={'geom2': 'geometry'})), 
                           how="inner", op='intersects', 
                           lsuffix='_route',
                        rsuffix = '_stop')

crosswalk_route_stop[['OBJECTID__route', 'OBJECTID__stop']].sample(5)

crosswalk_stop_acs = geopandas.tools.sjoin(
    geopandas.GeoDataFrame(bus_stops[['OBJECTID','geometry']]),
    geopandas.GeoDataFrame(acs[['GEOID','geometry']]), 
                           how="inner", op='intersects')

crosswalk_stop_acs[['OBJECTID', 'GEOID']].head()

crosswalk_route_stop = crosswalk_route_stop.merge(crosswalk_stop_acs, 
                           left_on = 'OBJECTID__stop', 
                           right_on = 'OBJECTID',
                           # Some stops are in DC , and some are not
                           how = 'left')[['OBJECTID__route', 'OBJECTID__stop', 'GEOID']]

crosswalk_route_stop.sample(10)

# Write to CSV for use in another step
crosswalk_route_stop.to_csv('data/crosswalk_route_stop.csv', index=False)

# Each stop has multiple OBJECTID__route values associated with it
crosswalk_route_stop.groupby('OBJECTID__stop').\
        OBJECTID__route.\
        apply(lambda x: x.unique().tolist()).head()

# Each route has multiple OBJECTID__stop values associated with it
crosswalk_route_stop.groupby('OBJECTID__route').\
        OBJECTID__stop.\
        apply(lambda x: x.unique().tolist()).head()

bus_routes.merge(crosswalk_route_stop, left_on  ='OBJECTID', right_on = 'OBJECTID__route').\
    groupby('ROUTE').OBJECTID__stop.nunique().\
    agg({'Average Stops per Route':np.mean, 
         'Minimum Stops on a Route': np.min, 
         'Maximum Stops on a Route': np.max,
         'Median Stops': np.median}).apply(round).to_frame()

bus_routes.merge(crosswalk_route_stop, left_on  ='OBJECTID', right_on = 'OBJECTID__route').\
    groupby('GEOID').ROUTE.nunique().\
    agg({'Average Routes':np.mean, 
         'Minimum Routes': np.min, 
         'Maximum Routes': np.max,
         'Median Routes': np.median}).apply(round).to_frame()

fig, ax = plt.subplots(1, figsize=(11,10))


# The block groups that intersect that bus line
base = acs.merge(bus_routes.merge(crosswalk_route_stop, left_on  ='OBJECTID', right_on = 'OBJECTID__route').\
    groupby('GEOID').ROUTE.nunique().to_frame().reset_index())[['geometry','ROUTE']].\
    plot(column='ROUTE', ax=ax,
         cmap=plt.cm.BrBG, legend=True, k = 8)
    
fig.suptitle("Number of Bus routes that cross through the block-group")
_ = ax.axis('off')