Wo befindet sich die "längste Gerade" im österreichischen Schienennetz?

Posted on So 28 Februar 2021 in Blog

Im Folder für den Streckenausbau Wien – Bratislava steht folgendes:

Beim vorliegendem Projekt handelt es sich um den Ausbau einer Bestandsstrecke mit der „längsten Geraden“ (rund 32,5 km) im österreichischen Eisenbahnnetz.

Daher dachte ich mir, das wäre ja nett einer Überprüfung zu unterziehen ;-) Mit ein wenig SQL und den Daten der OpenStreetMap ist eine Überprüfung auch möglich.

Inhalt:

Idee

In der OpenStreetMap ist ist das Liniennetz der ÖBB sehr gut erfasst. Daher sollte es doch eine Möglichkeit geben zu überprüfen, wo sich die "längste Gerade" befindet und wie lang diese ist.

Daher wurden die Daten auch in ein Schema at_osm in der Datenbank geladen. Erledigt wurde das mit osm2pgsql (siehe auch Routing von OSM-Daten mit OSRM und Postgres).

Um das Berechnung durchzuführen werden zuerst die entsprechenden Relationen aus den Daten der OpenStreetMap genommen. Danach die einzelnen Punkte extrahiert und der Winkel zueinander berechnet. Im nächsten Schritt wird ermittelt, ob der Winkel zwischen den einzelnen Zeilen gleich ist. Falls dem nicht so ist, ist das der Beginn einer neuen Gerade und für die Berechnung der Länge relevant. Von diesen Punkten ausgehend wird dann rekursiv die Ausdehnung berechnet. Danach wird noch die maximale Ausdehnung herangezogen und mittels des geography Datentyps die Länge ermittelt.

SQL

Um relativ leicht zu zusammenhängenden Streckenabschnitten zu kommen, habe ich mich entschieden Relationen, deren Tag route die Werte train oder railway hat, zu nehmen:

with recursive get_lines as (
    select osm_id, way, row_number() over() as tmp_id from at_osm.planet_osm_line 
    where route in ('train', 'railway')
),

Im nächsten Schritt werden die einzelnen Punkte einer Linie ermitlelt

dump_points as (
    select tmp_id, osm_id, path[1] as path_id, geom from get_lines as a inner join lateral st_dumppoints(way) as b on true 
),

Daraufhin wird die Lage bzw. der Winkel der Punkte zueinander ermittelt

calculate_azimuth as (
    select a.tmp_id, a.osm_id, a.path_id, 
    round(st_azimuth(a.geom, b.geom)::numeric, 1) as azimuth,a.geom
    from dump_points as a inner join dump_points as b on (a.osm_id=b.osm_id and a.path_id+1=b.path_id and a.tmp_id=b.tmp_id)
)

Jetzt wird überprüft, ob der Winkel davor gleich ist, wie der Aktuelle. Falls das nicht der Fall ist, dann haben wir den Beginn eines neuen gerade Streckenabschnitts

calculate_lag as (
    select tmp_id, osm_id, azimuth, geom, path_id,
    case when (azimuth IS DISTINCT FROM lag(azimuth) over (partition by tmp_id, osm_id order by path_id)) then true end as 
    calculate_recursion
    from calculate_azimuth
),

Jetzt wird rekursiv, ausgehend vom Beginn einer jeden "Gerade", die einzelnen Ausdehnungschritte berechnet

calculate_longest_line as (
    select tmp_id, osm_id, path_id as start_path_id, path_id, geom, geom as end_point, azimuth from calculate_lag where calculate_recursion=true
    union all
    select a.tmp_id, a.osm_id, b.start_path_id, a.path_id, b.geom, a.geom, a.azimuth from calculate_lag as a, calculate_longest_line as b
    where a.osm_id=b.osm_id and a.tmp_id=b.tmp_id and a.path_id=b.path_id+1 and a.calculate_recursion is null
),

Nun müssen wir noch ermitteln, was die maximale Ausdehnung ist

calculate_path_range as (
    select tmp_id, osm_id, start_path_id, max(path_id) as end_path_id from calculate_longest_line group by 1,2,3
),

Daraufhin kann die Länge berechnet werden. Hierfür wird einerseits der geography Datentyp herangezogen, um die tatsächliche Länge zu bekommen und um den Unterschied, der durch die Kartenprojektion entsteht, wird auch die Länge in der Projektion der OSM-Daten ermittelt.

calculate_length as (
    select a.tmp_id, a.osm_id, a.start_path_id, end_path_id, azimuth, 
    ST_Distance(ST_Transform(geom, 4326)::geography, ST_Transform(end_point, 4326)::geography) as length_geography, 
    ST_Distance(geom, end_point) as length_geom
    from calculate_longest_line as a,
    calculate_path_range as b where a.osm_id=b.osm_id and a.tmp_id=b.tmp_id and a.start_path_id=b.start_path_id
    and b.end_path_id=a.path_id
),

Jetzt wird noch die maximale Länge pro OSM-Relation berechnet. 

get_max_length as (
    select tmp_id, osm_id, max(length_geography) as max_length_geography,
    max(length_geom) as max_length_geom from calculate_length
    group by 1,2
)

Und zum Abschluss kann man sich Resultate sortiert ausgeben

select * from get_max_length order by max_length_geography desc

Zusammenfassung

Auch wenn der SQL-Code am Ende nicht übermäßig lang ist, so war der Weg dorthin durchaus ein wenig länger, da es doch einiges schwieriger war, als am Anfang gedacht. Und es ist auch noch nicht perfekt, da nicht jede Streck einmal vorkommt, sondern jede Linie, d.h. beispielsweise beinhalten die Relationen Graz-Prag und Ljubljana-Wien beide den Streckenabschnitt zwischen Graz Hauptbahnhof und Wien Hauptbahnhof.

Auch kann man natürlich darüber diskutieren, ob es sinnvoll ist auch tracks miteinzubeziehen (Beschreibung im OSM-Wiki). Generell macht, so denke ich, die Abfrage nach Relationen schon Sinn, da man so die gesamte Strecke bekommt. Ansonsten hat man oftmals nur Segmente und müsste diese dann verbinden. Ein weiterer Punkt, der noch genauer betrachtet kann, welche Unterschiede man in der Azimuth-Berechnung erlaubt.

Jedenfalls kommt man mit der durchgeführten Berechnung auf eine Länge von 31,9 km für den Streckenabschnitt zwischen Wien und Bratislava, was meiner Meinung nach schon ein sehr guter Wert ist und sehr nahe an dem, was die ÖBB in ihrem Folder angibt. Beziehungsweise auf 47,8 km für die Berechnung in der Web Mercator projection der OpenStreetMap. Das Explorer Magazin hat auch einen Artikel über die Feinheiten der Web Mercator Projection: "Web Mercator": Ein Blick hinter die Kartendienste im Web

Das gesamte Script hab ich als Github Gist zur Verfügung gestellt: osm-laengste-gerade.sql.

Postgres OpenStreetMap