Das Ergebnis

Das Video zeigt den Übergang vom realen Verlauf der Straßenbahnlinien zu dem Liniennetzplan der BSAG.

Motivation

Auf Reddit gibt es eine Sammlung von Animationen, die den Unterschied zwischen schematischen ÖPNV-Karten und der realen Geographie zeigen:

Ich möchte so eine Animation für das Straßenbahnnetz der BSAG in Bremen erstellen.

Dafür brauche ich die Pfade der Linien

aus dem Liniennetzplan
aus den Kartendaten von OpenStreetMaps (OSM)

Diese Pfade muss ich dann in ein gemeinsames Koordinatensystem überführen und dann den Übergang animieren.

Koordinaten der Linienführung

aus dem Liniennetzplan

Um den genauen Linienverlauf auf dem Netzplan nachbilden zu können, exportiere ich den Linienverlauf aus der PDF.

Vorbereitung

Inkscape installieren
Das Plugin ExportXY installieren
- Dazu die *.py und die *.inx Dateien in den Inkscape-Extension-Ordner kopieren
- Unter Windows ist der Pfad für den Extension-Ordner %AppData%\inkscape
Den Netzplan in Inkscape öffnen

Export der Koordinaten

Für jede Linie:

Eine neue Textdatei erstellen
Eine Straßenbahnlinie auswählen
Jetzt sollte die bounding box des Pfades ausgewählt sein. In der Werkzeugleiste werden dann die X und Y Koordinaten, Breite und Höhe der bounding box angezeigt. Beispiel:
Diese Werte als erste Zeile in die Textdatei eintragen, die Werte getrennt mit Tab
Jetzt das Pfad-Werkzeug auswählen
Dem Pfad Knoten hinzufügen mit Erweiterungen > Pfad modifizieren > Knoten hinzufügen
Die Koordinaten exportieren mit Erweiterungen > Exportieren > ExportXY
Die angezeigten Koordinaten kopieren und in die Textdatei einfügen

Die Datei von Linie 1 fängt dann z. B. so an:

574.500	1213.225	2108.693	688.133
0.0	0.0
0.710136	0.70221
1.41883	1.40711
2.1275	2.11266
...

Pfade aus OpenStreetMaps (OSM)

Die Linienführung aus den realen Karten kann aus OSM exportiert werden. Dazu benutze ich die Overpass API. Auf der Seite

Den gewünschten Kartenausschnitt wählen
Oben auf Wizard klicken
Als Suchbegriff route=tram eingeben und ausführen mit build and run query
Jetzt sollte die Linienführung auf der Karte sichtbar sein
Exportieren als GeoJSON

Daten zusammenführen

Um die Linien aus dem Netzplan und OSM zusammen darzustellen müssen sie in ein gemeinsames Koordinatensystem überführt werden. Die Koordinaten aus dem Liniennetzplan sind in relativen Bildkoordinaten. Um aus ihnen die absoluten Koordinaten zu berechnen habe ich die Maße der bounding box gespeichert. In der GeoJSON aus OSM sind die Linien als GPS-Koordinaten gespeichert.

Dafür benötige ich folgende Python Pakete

import pandas as pd
import geopandas as gpd
import numpy as np
from scipy import optimize as so

Straßenbahnlinien

Als erstes habe ich die Namen der Linien aus dem GeoJSON und die Farben aus dem Liniennetzplan herausgesucht

line = {
  1: {'name': 'Tram 1: Roland-Center => Bf Mahndorf',
      'color': '#00a650ff'},
  2: {'name': 'Tram 2: Sebaldsbrück => Gröpelingen',
      'color': '#1c63b7ff'},
  3: {'name': 'Tram 3: Gröpelingen => Weserwehr',
      'color': '#00adefff'},
  4: {'name': 'Tram 4: Lilienthal => Arsten',
      'color': '#ed1c24ff'},
  5: {'name': 'Tram 5: Bürgerpark => Bürgerpark',
      'color': '#0aabbaff'},
  6: {'name': 'Tram 6: Universität-Nord => Flughafen-Süd',
      'color': '#fbc707ff'},
  8: {'name': 'Tram 8: Kulenkampffallee => Roland-Center',
      'color': '#80cc28ff'},
  10: {'name': 'Tram 10: Sebaldsbrück => Gröpelingen',
      'color': '#2e3192ff'}
}

GPS Koordinaten laden

Die heruntergeladene GeoJSON habe ich in das Script-Verzeichnis als export.geojson gespeichert. Die Datei lade ich mit

routes = gpd.read_file('export.geojson')
for i in line.keys():
  route = routes[routes.name.str.match(line[i]['name'], na=False)]
  line[i]['gps'] = route.geometry.iloc[0].geoms[0]

und suche nach den minimalen und maximalen Koordinaten in allen Richtungen

gps_bounds = np.array([ r['gps'].bounds for r in line.values() ])
gps_xmin = np.min(gps_bounds[:, (0, 2)])
gps_xmax = np.max(gps_bounds[:, (0, 2)])
gps_ymin = np.min(gps_bounds[:, (1, 3)])
gps_ymax = np.max(gps_bounds[:, (1, 3)])

Linien aus Netzplan laden

Die aus Inkscape exportierten Koordinaten habe ich in einem Unterordner lines und benannt nach der Liniennummer, z. B. 1.txt, gespeichert. Die Dateien lade ich und transformiere die Koordinaten von relativen zu absoluten Bildkoordinaten. Dazu passe ich erst die Skalierung an die bounding box an und verschiebe sie dann an die richtige Position

for i in line.keys():
  path_bb = pd.read_csv(f'linien/{i}.txt', names=['x', 'y', 'w', 'h'],
                                index_col=False, sep='\t', nrows=1)
  path_rel = pd.read_csv(f'linien/{i}.txt', names=['x', 'y'],
                                index_col=False, sep='\t', skiprows=1)
  
  # Convert relative to absolute coordinates
  # Fix scaling of path
  _xmin = path_rel['x'].min()
  _xmax = path_rel['x'].max()
  path_abs_x = path_rel['x'].to_numpy() * path_bb['w'].iloc[0] / (_xmax - _xmin)
  
  _ymin = path_rel['y'].min()
  _ymax = path_rel['y'].max()
  path_abs_y = -path_rel['y'].to_numpy() * path_bb['h'].iloc[0] / (_ymax - _ymin)
  
  # Add offset from bounding box
  path_abs_x += path_bb['x'].iloc[0] - np.min(path_abs_x)
  path_abs_y += path_bb['y'].iloc[0] - np.min(path_abs_y)
  
  line[i]['path'] = np.stack((path_abs_x, path_abs_y)).T

Da die y-Achse in Bildkoordinaten in die falsche Richtung zeigt muss ich sie noch spiegeln.

path_ymax = np.min([ r['path']['y'].min() for r in line.values() ] )
for i in line.keys():
  line[i]['path']['y'] = path_ymax - line[i]['path']['y']

Wie bei den GPS-Koordinaten suche ich nach den minimalen und maximalen Koordinaten in allen Richtungen

path_bounds = []
for i in line.keys():
  _x = line[i]['path']['x']
  _y = line[i]['path']['y']
  path_bounds.append((_y.min(), _x.min(), _y.max(), _x.max()))

path_bounds = np.array(path_bounds)
path_xmin = np.min(path_bounds[:, (1, 3)])
path_xmax = np.max(path_bounds[:, (1, 3)])
path_ymin = np.min(path_bounds[:, (0, 2)])
path_ymax = np.max(path_bounds[:, (0, 2)])

Korrektur von Linie 5

Der Pfad von Linie 5 aus dem Netzplan ist nicht vollständig, da er eine Schleife enthält. Ich suche nahc dem Punkt, an dem sich die Schleife schließt und füge eine Kopie des restlichen Pfads an.

loop_idx = np.argmin(np.hypot(line[5]['path'][1:, 0] - line[5]['path'][0, 0],
                              line[5]['path'][1:, 1] - line[5]['path'][0, 1]))
line[5]['path'] = np.concatenate((line[5]['path'][:loop_idx:-1, :], line[5]['path']))[::-1]

Interpolation der GPS-Koordinaten

Um einfacher mit den GPS-Koordinaten arbeiten zu können, passe ich ihre Anzahl an die Anzahl der Bildkoordinaten an. Ich interpoliere die GPS-Koordinaten so, dass die neuen Koordinaten äquidistant sind.

for i in line.keys():
  steps = np.linspace(0, 1, line[i]['path'].shape[0])
  gps_interp = [ line[i]['gps'].interpolate(s, normalized=True).coords[:][0] for s in steps ]
  line[i]['gps_interp'] = np.array( gps_interp )

Koordinaten aus dem Netzplan transformieren

Um die Daten zusammen darzustellen transformiere ich die Bildkoordinaten aus dem Netzplan in den Wertebereich dere GPS-Koordinaten. Dafür gibt es zwei Möglichkeiten.

Übereinstimmende Kanten

Um die Koordinaten zu transfomieren kann ich annehmen, dass die nördlichsten, östlichsten, südlichsten und westlichsten Koordinaten übereinstimmen.

path_scale_x = (gps_xmax - gps_xmin)/(path_xmax - path_xmin)
path_scale_y = (gps_ymax - gps_ymin)/(path_ymax - path_ymin)

for i in line.keys():
  path_gps_x = (line[i]['path'][:, 0] - path_xmin) * path_scale_x + gps_xmin
  path_gps_y = (line[i]['path'][:, 1] - path_ymin) * path_scale_y + gps_ymin
  line[i]['path_gps'] = np.stack((path_gps_x, path_gps_y)).T

In einer Animation würden sich in dem Fall die Punkte an den Kanten recht wenig bewegen. Ein Großteil der bewegungen würde in der Mitte des Bildes stattfinden.

Geringste Distanz

Alternativ kann ich eine Metrik einführen, wie gut die beiden Koordinaten nach der Transformation übereinander liegen. Dann kann ich mit der Funktion interpolate() aus dem Paket scipy.optimize kann die optimalen Parameter finden, die diese Metrik minimieren.

Die Metrik berechne ich als Summe der Distanzen zwischen den Punkten der jeweiligen Koordinatensystemen.

def dist(a, line):
    dist = 0
    for i in line.keys():
      path_gps_x = line[i]['path'][:, 0] * a[0] + a[1]
      path_gps_y = line[i]['path'][:, 1] * a[2] + a[3]
      
      gps_x = line[i]['gps_interp'][:, 0]
      gps_y = line[i]['gps_interp'][:, 1]
      
      point_dists = np.hypot(gps_x - path_gps_x, gps_y - path_gps_y)
      dist += np.sum(point_dists)
      
    return dist

Diese Metrik minimiere ich

  opt = so.minimize(dist, (1, 0, 1, 0), args=(line))

und transformiere dann die Koordinaten mit den ermittelten Parametern

  for i in line.keys():
    path_gps_x = line[i]['path'][:, 0] * opt.x[0] + opt.x[1]
    path_gps_y = line[i]['path'][:, 1] * opt.x[2] + opt.x[3]
    line[i]['path_gps'] = np.stack((path_gps_x, path_gps_y)).T

Richtung der Linien

Als letztes muss ich noch sicherstellen, dass die Daten alle den richtigen Start- und Endpunkt haben. Auf dem Bild ist die Linienführung im Netzplan und OSM dargestellt, jeweils mit den Startpunkten in rot und den Endpunkten in blau. Wie man sieht, stimmen diese bei einigen Linien nicht überein. Ich drehe die Reihenfolge der Koordinaten um, wenn Start- und Endpunkt näher aneinander sind als die jeweiligen Startpunkte.

for i in lines.keys():
  dist_first = np.hypot(lines[i]['gps_interp'][0, 0] - lines[i]['path_gps'][0, 0],
                        lines[i]['gps_interp'][0, 1] - lines[i]['gps_interp'][0, 1])
  dist_last = np.hypot(lines[i]['gps_interp'][-1, 0] - lines[i]['path_gps'][0, 0],
                        lines[i]['gps_interp'][-1, 1] - lines[i]['gps_interp'][0, 1])
  if dist_last < dist_first:
    lines[i]['path_gps'][:, :] = lines[i]['path_gps'][::-1, :]

Linienführung im Netzplan und OSM nach Korrektur, jeweils mit dem Startpunkt in rot und dem Endpunkt in blau. Nach der Korrektur sind die Start- und Endpunkte in der richtigen Reihenfolge. In der Animation erstelle ich jetzt einen Übergang vom rechten in des linke Bild.

Animation

Interpolation des Übergangs

Ich habe jetzt für jede Linie zwei Arrays: Koordinaten im Netzplan und in OSM. Beide Arrays haben jeweils die gleiche Anzahl an äquidistanten Punkten. Ich kann jetzt durch lineare Interpolation für jeden Punkt im Übergang Punkt in OSM -> Punkt in Netzplan die Positionen berechnen.

frames = 150
delay = 50
for i in line.keys():
  anim_path_x = np.linspace(line[i]['gps_interp'][:, 0], line[i]['path_gps'][:, 0], frames)
  anim_path_y = np.linspace(line[i]['gps_interp'][:, 1], line[i]['path_gps'][:, 1], frames)
  
  path_anim = np.stack((anim_path_x, anim_path_y)).T
  start_delay = np.repeat(line[i]['gps_interp'][:, None, :], delay, axis=1)
  end_delay = np.repeat(line[i]['path_gps'][:, None, :], delay, axis=1)
  line[i]['path_anim'] = np.concatenate((start_delay, path_anim, end_delay), axis=1)

Hier sind frames die Anzahl der Bilder im Übergang und delay die Anzahl Bilder, in denen am Anfang und am Ende ein Standbild stehen bleibt.

Plot

Die Animation erstelle ich dann mit Matplotlib. Für jeden Frame tausche ich die Daten der Linien im Plot.

from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots(figsize=(8, 8))
line_p = []
for i in line.keys():
  line_p.append(ax.plot(line[i]['path_anim'][:, 0, 0], line[i]['path_anim'][:, 0, 1],
                        color=line[i]['color'])[0])
ax.axis('off')

def animate(frame):
  global line, delay
  for i, ln in enumerate(line_p):
    ln.set_data(list(line.values())[i]['path_anim'][:, frame, 0], 
                list(line.values())[i]['path_anim'][:, frame, 1])
  return line_p

plt.tight_layout()
anim = FuncAnimation(fig, animate, frames=frames+2*delay, interval=20, blit=False)
plt.show()

Das gesamte Skript und die Dateien sind in diesem Repository zu finden.

BSAG Netzkarte Animation