Datos de Automated Fare Collection del sistema de transporte público de Bahía Blanca

Este documento contiene código fuente y notas preliminares para un estudio de algunos aspectos de la utilización del sistema de transporte público en Bahía Blanca. En esta etapa, nos proponemos una exploración de la información disponible que informe futuros análisis.

Los posibles trabajos sobre esta información incluyen, pero no se limitan a:

  • Impacto del precio del boleto en la economía de los pasajeros (vinculando origen y destino de los viajes con indicadores socioeconómicos georeferenciados)
  • Caracterización de los perfiles de usuarios del sistema.
  • Eficiencia del sistema (ociocidad de las unidades, etc)
  • Impacto ambiental
  • Adherencia a los horarios planificados.
  • Planificación de horarios, servicios y rutas.
  • Medidas de demanda, oferta y performance de los servicios existentes.
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import psycopg2
import pandas as pd
from mpl_toolkits.basemap import Basemap
from matplotlib.collections import LineCollection
import fiona
from shapely.geometry import Point, LineString
from shapely.wkt import loads as wkt_load
from itertools import izip
import math

db = psycopg2.connect('dbname=bondis_bahia user=manuel')
colors = ['#1f77b4','#ff7f0e','#2ca02c','#d62728','#9467bd','#8c564b','#e377c2','#7f7f7f','#bcbd22','#17becf']

Pre-procesamiento de los datos

Los datos provistos por Gobierno Abierto de la Municipalidad de Bahía Blanca se procesaron para guardarlos en una tabla Postgres/PostGIS:

CREATE TABLE omnibus_usuarios (
    id integer NOT NULL,
    empresa text,
    tipo_pasaje text,
    linea text,
    recorrido integer,
    id_tarjeta integer,
    interno integer,
    id_colectivo integer,
    fechahora timestamp without time zone,
    pasajeros integer,
    locacion geometry(Point,4326)
);

Perfil temporal de un usuario

Consideramos a un usuario como una observación de 7 * 24 (días de la semana y horas del día) variables, sobre el período disponible. Lo caracterizamos a través de su perfil semanal, construído como la suma de los viajes realizados en un día y hora. Esta construcción será útil para posibles caracterizaciones de tipos o grupos de usuarios.

El concepto de perfil temporal está basado en la publicación Understanding Passenger Patterns in Public Transit Through Smart Card and Socioeconomic Data (2014).

El perfil se construye mediante la siguiente función SQL:

CREATE OR REPLACE FUNCTION get_user_temporal_profile(id_tarjeta integer, round_interval interval)
RETURNS table(id_tarjeta integer, day_of_week integer, hour double precision, trips integer)
AS $$ 

SELECT id_tarjeta,
       cast(extract(dow
                    FROM fechahora) AS integer) AS day_of_week,
       (date_part('hour', date_round(fechahora, $2)) * 100 + date_part('minute', date_round(fechahora, $2))) AS hour,
       cast(count(*) AS integer) AS trips
 FROM omnibus_usuarios
WHERE id_tarjeta = $1
GROUP BY id_tarjeta,
         hour,
         day_of_week
ORDER BY day_of_week;
$$
LANGUAGE sql;

La siguiente función en lenguaje Python grafica un perfil de usuario:

In [197]:
def plot_temporal_profile(weekly_profile, xsize=15):
    a = np.zeros((7,24))
    plt.figure(figsize=(12,24/7.))
    ax = plt.gca()
    plt.xticks(range(0,24))
    t = ax.transData
    canvas = ax.figure.canvas
    plt.title("Perfil del usuario: %s" % weekly_profile['id_tarjeta'][0])
    for dow, hour, trips in zip(weekly_profile['day_of_week'], weekly_profile['hour'] / 100, weekly_profile['trips']):
        a[dow][int(hour)] = trips
        text = ax.text(int(hour) - 0.25, dow + 0.1, trips, transform=t, color='black')
        text.draw(canvas.get_renderer())
    plt.imshow(a,interpolation='nearest', cmap='Blues')

Observamos, a modo de ejemplo, los perfiles temporales para dos usuarios elegidos aleatoriamente.

El eje X contiene las horas del día, de 0 a 23. El eje Y muestra los días de la semana, con Domingo = 0 y Sábado = 6. La intensidad del color representa la cantidad de viajes efectuados por el usuario en un día de la semana y hora determinada, dentro del período estudiado (15 de julio de 2014 hasta 15 de septiembre del mismo año).

En el primero, notamos un patrón uniforme en el que presumiblemente el pasajero usa el transporte público alrededor de las 7 de la mañana y retorna a su casa alrededor de las 18. Hay uso ocasional durante los fines de semana.

El segundo comparte la característica del viaje matutino, pero se observa un uso consistente durante el mediodía, probablemente a causa del horario comercial "cortado" que se acostumbra en Bahía Blanca.

In [199]:
plot_temporal_profile(pd.read_sql("select * from get_user_temporal_profile(11061406, '1 hour')", db))
plot_temporal_profile(pd.read_sql("select * from get_user_temporal_profile(99019579, '1 hour')", db))

Promedio de cantidad total de viajes por hora para el mes de agosto

Los siguientes gráficos muestran la distribución de la cantidad promedio de viajes por hora, considerando lunes a viernes. En primer lugar, se muestra dicha estadística para todos los tipos de pasaje. Luego, se muestran para los tipos Pesos, Frecuente, Libre y Viaje.

Se observa una distribución similar en todos los casos, con picos de actividad alrededor de las 7 de la mañana, el mediodía y el fin de la jornada laboral.

In [ ]:
avg_trips_per_hour = pd.read_sql("""with viajes_por_dia_y_hora as (
select date_trunc('day', fechahora) as dia, date_part('hour', fechahora) as hora, sum(pasajeros) as p
  from omnibus_usuarios
 where date_part('month', fechahora) = 8
   and extract(dow from fechahora) <> 0
   and extract(dow from fechahora) <> 6
 group by dia, hora)

select hora, avg(p) from viajes_por_dia_y_hora
group by hora
order by hora""", db)

avg_trips_per_hour_by_type = pd.read_sql("""select date_trunc('day', fechahora) as dia, date_part('hour', fechahora) as hora, tipo_pasaje, sum(pasajeros) as p
  from omnibus_usuarios
 where date_part('month', fechahora) = 8
   and extract(dow from fechahora) <> 0
   and extract(dow from fechahora) <> 6
 group by dia, hora, tipo_pasaje""", db)
In [6]:
plt.figure(figsize=(15,15))
plt.subplot2grid((3,2), (0,0), colspan=2)
plt.title('Todos los tipos de pasaje')
plt.ylim(0,10000)
avg_trips_per_hour['avg'].plot(kind='bar', color=colors[0])

tipos_de_pasaje = ('Pesos', 'Frecuente', 'Libre', 'Viaje')
for i, t in enumerate(tipos_de_pasaje):
    plt.subplot2grid((3,2), ((i / 2) + 1,i % 2))
    plt.ylim(0,10000)
    plt.xlim(0,23)
    plt.title('Tipo de pasaje: %s' % t)
    avg_trips_per_hour_by_type[avg_trips_per_hour_by_type['tipo_pasaje'] == t].groupby(['hora'])['p'].mean().plot(kind='bar', color=colors[i+1])
    
plt.savefig('avg_trips.png', format='png')

Viajes encadenados

Decimos que un pasajero p realizó dos viajes encadenados v1 y v2 si:

  • el intervalo de tiempo entre v1 y v2 es menor a 45 minutos
  • las línea de ómnibus de los viajes v1 y v2 son distintas (para descartar los viajes de retorno al mismo lugar)
  • v1 y v2 sucedieron en la misma fecha
In [2]:
encadenados_cte = """with encadenados as (select o1.fechahora as fechahora1, o1.locacion as locacion1, o1.linea as linea1, o2.fechahora as fechahora2, o2.locacion as locacion2, 
o2.linea as linea2, 
st_distance_sphere(o1.locacion, o2.locacion) as distancia
from omnibus_usuarios o1
inner join omnibus_usuarios o2
 on o2.id_tarjeta = o1.id_tarjeta
where abs(extract(epoch from o2.fechahora - o1.fechahora)) < 45*60 -- menos de 45m entre viajes consecutivos 
  and o2.fechahora > o1.fechahora
  and o1.linea <> o2.linea -- distintas lineas
  and date_trunc('day', o2.fechahora) = date_trunc('day', o1.fechahora) -- misma fecha
order by distancia desc)"""
In [ ]:
encadenados = pd.read_sql(encadenados_cte + 
"""select linea1, linea2, count(*) as c from encadenados
group by linea1, linea2
order by linea1, linea2""", db)

El siguiente gráfico representa la matriz de frecuencias de viajes encadenados para el período disponible (15/7 a 15/9). Se observa que la combinación de viajes en las líneas 514 y 517 sucede frecuentemente, sugiriendo que hay una población considerable de usuarios en las zonas de influencia de dichas líneas que deben combinar servicios para llegar a destino.

In [225]:
lineas = ['500','502','503','504','505','506','507','509','512','513','513E','514','516','517','518','519','519A']
pares = np.zeros((17, 17))
for linea1, linea2, viajes in zip(encadenados['linea1'], encadenados['linea2'], encadenados['c']):
    pares[lineas.index(linea1)][lineas.index(linea2)] = viajes

plt.figure(figsize=(10,10))
plt.xticks(range(0,17), lineas)
plt.yticks(range(0,17), lineas)
plt.imshow(pares,interpolation='nearest', cmap='Reds')
c = plt.colorbar()

Vale la pena visualizar la distribución espacial de las paradas de origen de los usuarios que realizan habitualmente la combinación 514-517. El mapa sugiere un trayecto frecuente, que no está conectado por niguna línea de ómnibus.

In [101]:
import itertools

def plot_encadenados(linea1, linea2):
    encadenados_l1_l2 = pd.read_sql("""
with encadenados as (select o1.fechahora as fechahora1, o1.locacion as locacion1, o1.linea as linea1, o2.fechahora as fechahora2, o2.locacion as locacion2, 
o2.linea as linea2, 
st_distance_sphere(o1.locacion, o2.locacion) as distancia
from omnibus_usuarios o1
inner join omnibus_usuarios o2
 on o2.id_tarjeta = o1.id_tarjeta
where abs(extract(epoch from o2.fechahora - o1.fechahora)) < 45*60 -- menos de 45m entre viajes consecutivos 
  and o2.fechahora > o1.fechahora
  and o1.linea = %s and o2.linea = %s -- distintas lineas
  and date_trunc('day', o2.fechahora) = date_trunc('day', o1.fechahora) -- misma fecha
order by distancia desc)
    select linea1, linea2, date_round(fechahora1, '30 minutes') as fechahora, count(*) as c, 
    st_x(locacion1) as longitud1, st_y(locacion1) as latitud1,
    st_x(locacion2) as longitud2, st_y(locacion2) as latitud2
    from encadenados
    where linea1 = %s and linea2 = %s
      and date_part('hour', date_round(fechahora1, '30 minutes')) < 9
    group by linea1, linea2, locacion1, locacion2, fechahora
    """, params=(linea1, linea2, linea1, linea2), con=db)
    
    shp = fiona.open('shps/bahia_4326_2.shp')
    bds = shp.bounds
    shp.close()
    extra = 0.01
    ll = (bds[0], bds[1])
    ur = (bds[2], bds[3])
    coords = list(itertools.chain(ll, ur))
    w, h = coords[2] - coords[0], coords[3] - coords[1]


    plt.clf()
    fig = plt.figure(figsize=(18,18))
    ax = fig.add_subplot(111, axisbg='w', frame_on=False)

    m = Basemap(
        projection='tmerc',
        lon_0=-62.2584,
        lat_0=-38.6879,
        ellps = 'WGS84',
        llcrnrlon=coords[0] - extra * w,
        llcrnrlat=coords[1] - extra + 0.01 * h,
        urcrnrlon=coords[2] + extra * w,
        urcrnrlat=coords[3] + extra + 0.01 * h,
        lat_ts=0,
        resolution='i',
        suppress_ticks=True)

    # transformar puntos a las coordenadas proyectadas del mapa
    map_points_514 = pd.Series(
        [Point(m(mapped_x, mapped_y)) 
         for mapped_x, mapped_y in zip(encadenados_l1_l2['longitud1'], encadenados_l1_l2['latitud1'])])

    map_points_517 = pd.Series(
        [Point(m(mapped_x, mapped_y)) 
         for mapped_x, mapped_y in zip(encadenados_l1_l2['longitud2'], encadenados_l1_l2['latitud2'])])

    m.scatter(
        [geom.x for geom in map_points_514],
        [geom.y for geom in map_points_514],
        s=[math.sqrt(trips) * 20 for trips in encadenados_l1_l2['c']], marker='o', lw=.55,
        facecolor='red', edgecolor='w',
        alpha=1, antialiased=True, 
        zorder=3)

    m.scatter(
        [geom.x for geom in map_points_517],
        [geom.y for geom in map_points_517],
        s=[math.sqrt(trips) * 15 for trips in encadenados_l1_l2['c']], marker='o', lw=.55,
        facecolor='blue', edgecolor='w',
        alpha=1, antialiased=True,
        zorder=3)

    x = pd.read_sql("select linea, st_astext(recorrido) as r from omnibus_recorridos where linea = %s or linea = %s", 
                    params=(linea1, linea2),
                    con=db)
    for r,linea in izip([wkt_load(g) for g in x[u'r']], x['linea']):
        col = LineCollection([[m(*c) 
                               for c in r.coords]])
        col.set_color((1,0,0,0.5) if linea == linea1 else (0,0,1,0.5))
        ax.add_collection(col)

    d = m.readshapefile(
        'shps/bahia_4326_2',
        'bahia_4326',
        color=(0.1,0.1,0.1,0.3),
        zorder=2)
    
In [102]:
plot_encadenados('517', '514')
# plt.savefig('mapa.pdf', format='pdf')
<matplotlib.figure.Figure at 0x110ce5890>