Render archived MTA NYCT A Division ATS data

subviewer.pde

import processing.video.*;

MovieMaker mm;
OSMMercator om;
PFont theFont;


void setup() {
  size(1280, 720);
  smooth();
  fill(0);
  noStroke();
  background(255);
  mm = new MovieMaker(this, width, height, "drawing.mov", 
  30, MovieMaker.ANIMATION, MovieMaker.HIGH);

  om = new OSMMercator(40.756, -73.987, 1.05/width, width, height);

  theFont = createFont("HelveticaNeue", 48);
  textFont(theFont, 48);

  String lines[] = loadStrings("out.csv");

  String[] boroughs = loadStrings("shapeout");

  drawBG(boroughs);

  int lastFrame = 0;
  drawClock(lastFrame);

  for (int i=0; i < lines.length; i++) {
    String[] parts = splitTokens(lines[i], ",");
    int frame = int(parts[0]);
    float latitude = float(parts[1]);
    float longitude = float(parts[2]);
    String fillcolor = parts[3];
    int direction = int(parts[4]);

    if (frame > lastFrame) {
      for (int j=0; j < (frame - lastFrame); j++) {
        mm.addFrame();
        background(255);
        drawBG(boroughs);
        drawClock(lastFrame+j);
      }
    }

    color clr = color(unhex(fillcolor.substring(1, 3)), 
    unhex(fillcolor.substring(3, 5)), 
    unhex(fillcolor.substring(5, 7)), 255);  

    fill(clr);
    if (direction == 0) { //northbound, right
      arc((float)om.x(longitude), (float)om.y(latitude), 10, 10, radians(-90), radians(90));
    } 
    else { //southbound, left
      arc((float)om.x(longitude), (float)om.y(latitude), 10, 10, radians(90), radians(270));
    }
    fill(0);

    lastFrame = frame;
    /*if (frame > 500) {
     break;
     }*/
  }
  mm.addFrame();
  mm.finish();
  exit();
}

void drawClock(int frameCount) {
  int total_seconds = frameCount * 30;

  int hours = total_seconds / 3600;
  int minutes = (total_seconds % 3600) / 60;
  int seconds = (total_seconds % 3600) % 60;

  String out = nf(hours, 2) + ":" + nf(minutes, 2) + ":" + nf(seconds, 2);
  textAlign(RIGHT, CENTER);
  text(out, 1280, 720/2);
}

void drawBG(String[] boroughs) {
  for (int i=0; i < boroughs.length; i++) {
    String[] segments = splitTokens(boroughs[i], ";");


    for (int j=0; j < segments.length; j++) {
      String[] start_point = splitTokens(segments[j], ",");


      float s_x = (float)om.x(float(start_point[0]));
      float s_y = (float)om.y(float(start_point[1]));

      String end_point;
      if (j+2 > segments.length) {
        end_point = segments[0];
      } 
      else {
        end_point = segments[j+1];
      }

      String[] end_point_val = splitTokens(end_point, ",");

      float e_x = (float)om.x(float(end_point_val[0]));
      float e_y = (float)om.y(float(end_point_val[1]));

      stroke(0);
      line(s_x, s_y, e_x, e_y);
      noStroke();
    }
  }
}

subviewer.py

from csv import DictReader, DictWriter
from collections import defaultdict
import xml.etree.ElementTree as etree
from datetime import datetime
import math


import pytz

from gtfs import Schedule
from gtfs.entity import Stop, Route

eastern = pytz.timezone('US/Eastern')

def total_seconds(td):
    return int(math.floor((td.days * 3600 * 24) + td.seconds + (td.microseconds / 100000.0)))


dr = DictReader(open('data.csv'))
schedule = Schedule('google_transit.db')

trains = defaultdict(list)
timestamps = []

for row in dr:
    trains[row['train_id']].append(row)
    timestamps.append(eastern.localize(datetime.strptime(row['timestamp'], "%Y-%m-%d %H:%M:%S")))
    
print len(trains)

dr = DictReader(open('data.csv'))
min_ts =  min(timestamps)
max_ts =  max(timestamps)

points = []
frames = defaultdict(set)

for train, events in trains.iteritems():
    stops = defaultdict(list)
    try:
        color = '#' + Route.query.get(events[0]['route_id']).route_color
    except Exception:
        color = '#000000'

    direction = events[0]['direction_id']
        
    for event in events:
        stops[event['stop_id']].append(event)

    for stop, stop_events in stops.iteritems():
        assert len(stop_events) <= 2

        arr = dep = None
        for event in stop_events:
            if event['event_type'] == '1':
                arr = event['timestamp']
            if event['event_type'] == '2':
                dep = event['timestamp']
            if arr is None and dep is not None:
                arr = dep
            if dep is None and arr is not None:
                dep = arr

        dep = eastern.localize(datetime.strptime(dep, "%Y-%m-%d %H:%M:%S"))
        arr = eastern.localize(datetime.strptime(arr, "%Y-%m-%d %H:%M:%S"))

        stop_data = Stop.query.get(stop)

        start_frame = total_seconds(arr - min_ts)/30
        end_frame = total_seconds(dep - min_ts)/30

        for frame in xrange(start_frame, end_frame + 1):
            frames[frame].add((stop_data.stop_lat, stop_data.stop_lon, color, direction))
            """frames.append({'latitude': stop_data.stop_lat,
                           'longitude': stop_data.stop_lon,
                           'frame': frame,
                           'color': color})"""

dw = DictWriter(open('out.csv', 'w'), ['frame', 'latitude', 'longitude', 'color', 'direction'])

for frame, points in frames.iteritems():
    for point in points:
        dw.writerow({'frame': frame,
                     'latitude': point[0],
                     'longitude': point[1],
                     'color': point[2],
                     'direction': point[3]
                  })

transform.py

import osgeo.osr as osr
import shapefile

base = "nybb_10c_av/nybb"

prj = open(base+".prj").read()

srs_in = osr.SpatialReference()
srs_in.ImportFromWkt(prj)
srs_in.MorphFromESRI()

srs_out = osr.SpatialReference()
srs_out.SetWellKnownGeogCS("WGS84")

ct = osr.CoordinateTransformation(srs_in, srs_out)

sf = shapefile.Reader("nybb_10c_av/nybb")

with open("shapeout", "w") as outfile:
    for shape in sf.shapes():

        lastStartIndex = 0
        print shape.parts
        for partStartIndex in shape.parts:
            print lastStartIndex, partStartIndex
            points = shape.points[lastStartIndex:partStartIndex]
            lastStartIndex = partStartIndex
            if len(points) > 0:
                point_out = ct.TransformPoints(points)
                outfile.write(";".join(['%s,%s' % p[0:2] for p in point_out])+'\n')