TL;DR
Download der Rohdaten: hier
Osmosis-Befehl: ./osmosis/bin/osmosis --read-pbf germany-latest.osm.pbf --tf reject-relations --tf accept-ways railway=rail --used-node --write-xml filtered_data.osm
Download des Python-Skriptes: hier
Download des fertigen Shapefiles: Shapefile, Karte, PDF, Karte, PNG, QGIS-Datei (Stand 2022-12-05, © OpenStreetMap-Mitwirkende)
Nachtrag: Europa-Karte, PDF
Inspiration
Auch wenn es Websites wie die openrailwaymap mit ihrem Elektrifizierungslayer oder den DB-GeoViewer für die Gleisanzahl gibt: Beide bieten keine praktische Funktion, beide Metriken miteinander zu vereinen. Aus einem Thread auf Twitter inspiriert war das Ziel also, beides auf einer von üblichen Computern händelbaren Karte zu vereinen.
Also gings ran an die Arbeit. Da ein Export von den diversen Overpass-Servern in den (Flächen-)Größenordnungen außer Disussion stand, lag die Arbeit mit einem direkten Dump aus der OSM nahe. Zu Testzwecken nutze ich dafür gerne einen kleineren Auszug, nämlich den von Tübingen. Der ist im Vergleich zu den knapp 4GB des Deutschland-Dumps mit ~100MB etwas handlicher für die ersten Tests.
Da die Bahnspezifischen Daten nur einen kleinen Teil des Datensatzes ausmachen lag es nahe, zuersteinmal alles irrelevante herauszufiltern. Hierfür nutzte ich das Tool “Osmosis”, das sich mit wenigen Befehlen schon vorkompiliert auf dem Rechner “installieren” lässt.
Die eigentliche Filterung erfolgt wie folgt:
./osmosis/bin/osmosis \ # Aufrufen des Tools
--read-pbf germany-latest.osm.pbf \ # einlesen des OSM-Exportes im .pbf-Format
--tf reject-relations \ # Relationsdaten löschen
--tf accept-ways railway=rail \ # nur die "ways" übernehmen, die mit "railway=rail" getaggt sind
--used-node \ # alle nodes, die von den "ways" referenziert wurde auch mit übernehmen
--write-xml filtered_data.osm # die Ergebnisse in die Ausgabedatei ausgebenDer Prozess selbst dauert dann einen Moment. Für meine Tübingen-Testdatei waren es gut 15 Sekunden, das Runterbrechen des Deutschland-Exportes brauchte dagegen schon gut zehn Minuten. Diese Datei in Openstreetmap’s XML-Format lässt sich zwar schon z.B. in QGIS öffnen, ist aber gerade im großen Deutschland-Export noch arg laggy und unbequem auszuwerten.
Über die python-Library des Osmium Tools wird deshalb der gefilterte Export noch weiter von unnötigem “Informationsballast” befreit und per “pyshp”-Library zu einem kompakten Shapefile (~700kB für Tübingen, 40MB für Deutschland) eingedampft.
# installation of all dependencies
!pip3 install pyshp
!pip3 install osmium
import osmium
import shapefile
import timeOsmium arbeitet mit einer Handler-Klasse, die dann die geforderten Aktionen sehr effizient abarbeitet. In der Init-Funktion wird, für den späteren Export als Shapedatei, die pyshp-Bibliothek entsprechend initialisiert und die zwei Datenfelder erzeugt (“POWERLINE” und “TRACK_CNT”, leider gibt es dort eine Zeichenbegrenzung der Spaltennamen).
class Worker(osmium.SimpleHandler):
def __init__(self, output_filename):
osmium.SimpleHandler.__init__(self)
self.filename = output_filename.replace(".shp", "")
self.shapes = shapefile.Writer(self.filename + ".shp", shapeType=shapefile.POLYLINE)
self.shapes.autoBalance = 1
self.shapes.field("POWERLINE", "C")
self.shapes.field("TRACK_CNT", "C")Die way()-Funktion wird von der Bibliothek bei jedem geparsten way aufgerufen.
Diese liest die zwei gefragten Tags aus, extrahiert aus den assoziierten Nodes die Koordinaten und speichert diese zusammen mit den Elektrifizierungs und “Gleisigkeits”-Metadaten direkt in das Shapefile.
def way(self, w):
electrified = w.tags.get("electrified")
track_count = w.tags.get("passenger_lines")
coords = []
for node in w.nodes:
if node.location.valid:
coords.append([node.location.lon, node.location.lat])
self.shapes.line([coords])
self.shapes.record(electrified, track_count)Shapefiles sind per-se Koordinatensystemlos, das heißt sie müssen für die weitere Verarbeitung mit QGIS damit angereichert werden. Da die Openstreetmap EPSG-4326 benötigt, wird der entsprechende Passus in die Projektions (.prj)-Datei abgelegt.
def saveall(self):
# see https://spatialreference.org/ref/epsg/4326/prettywkt/ for coordinate-system
epsg = """GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.01745329251994328,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]]"""
with open(self.filename + ".prj", "w") as f:
f.write(epsg)
f.close()
self.shapes.close()Jetzt ist der Worker/Handler soweit fertig bestückt und kann ausgeführt werden. Der Vorgang braucht selbst für den (reduzierten) Deutschland-Datensatz nur um die fünfzehn Sekunden.
runner = Worker("shapefile_germany.shp")
runner.apply_file("filtered_data.osm", locations=True)
runner.saveall()Im entsprechenden Verzeichnis liegen nun die vier für die Shapedatei benötigten Files.
Visualisierung
Mit QGIS können nun die Daten relativ performant visualisiert werden. Hierfür wird das Shapefile über “Layer > Layer hinzufügen > Vektorlayer hinzufügen” und der entsprechenden Angabe des Pfades zur .shp-Datei hinzugefügt.
Per Rechtsklick auf die “shapefile”-Zeile im Layer-Fenster und die Auswahl von “Eigenschaften” kann die Symboldarstellung angepasst werden (siehe Screenshot)
Über die Auswahl von “Kategorisiert” und der Variable “POWERLINE” im Bereich “Darstellung” können die Farben für die Elektrifizierung eingerichtet werden. Per Klick auf die “Symbol”-Zeile kann in den Symboleinstellungen zusätzlich die Linienbreite in Abhängigkeit von der Variable “TRACK_CNT” gesetzt werden.
Ein klick auf “Anwenden” aktualisiert die Darstellung in der QGIS-Karte die ggfs. um einen Hintergrundlayer wie z.B. die Openstreetmap-Grundkarte erweitert werden kann.
Download der Karte als PDF
Nachtrag
Dank der Option, auf einen Rechner mit mehr ressourcen zugreifen zu können, entstand außerdem eine Europa-Karte aus demselben Datensatz, nur in etwas größer. Die Karte siehe hier