Routing mit einem dichten Punktnetz

Zum Thema Routing mit Daten der OpenStreetMap gab es in letzter Zeit schon zwei Posts:

• Routing von OSM-Daten mit OSRM und Postgres
• Naive räumliche Interpolation in Postgres

Bei beiden war das "Problem", dass zwischen den Punkten interpoliert wurde und es so aussieht, als ob jeder Punkt in Österreich mit dem Auto erreichbar wäre. Was ja nicht ganz den Tatsachen entspricht. Die Inspiration für den nächsten Versuch habe ich mir über das QGIS-Plugin QNEAT3 geholt.

Inhalt

• Idee
• Vorbereitungen
• SQL
• Zusammenfassung

Idee

Als Basis werden die Daten der OpenStreetMap und die Klassifikation des Urban Audit (2020) von Eurostat dienen. Im ersten Schritt werden die entsprechenden Straßentypen aus der OpenStreetMap ausgewählt, die für ein Routing mittels Auto in Frage kommen. Bei den Straßen handelt es sich um Multiline-Objekte. Daher kommt die Idee, die einzelnen Punkte für das Routing heranzuziehen. So entsteht ein sehr feines Netz in Regionen mit hoher Straßendichte, während es im alpinen Raum spärlicher sein wird.

Vorbereitungen

Als Daten dienen die Daten der OpenStreetMap, die bei mir im Schema at_osm gespeichert sind und mittels osm2pgsql eingespielt wurden. Die Daten des Urban Audit liegen im Schema at_routing. Sind die Punktdaten von Eurostat, die dann nur die österreichischen Orte umfassen.

SQL

Im ersten Schritt werden alle Straßen, die für das Routing in Frage kommen aus den Daten der OpenStreetMap extrahiert und mittels lateral join die einzelnen Punkte mittels ST_DumpPoints extrahiert. Danach werden alle extrahierten Punken mit allen Städten des Urban Audit verknüpft. Das Resultat der Abfrage wird auch gleich in einer Tabelle gespeichert. Disclaimer: aus Gründen der Faulheit und Einfachheit habe ich auch die Punkte dazugespeichert. Ist nicht unbedingt der ideale Ansatz, aber bei der überschaubaren Datenmenge ist es nicht das große Problem

create table at_routing.urban_edges as
with select_osm_points as (
    select osm_id, x.geom from 
    at_osm.planet_osm_line left join lateral ST_DumpPoints(st_points(way)) as x on true 
    where highway in ('motorway', 'trunk', 'primary', 'secondary', 'tertiary', 'unclassified', 'residential', 
    'motorway_link', 'trunk_link', 'primary_link', 'secondary_link', 'tertiary_link', 'living_street')
)
select id, b.geom, osm_id, a.geom as target_geom from select_osm_points as a, at_routing.cities_at_urban_audit as b;

Im nächsten Schritt wird noch eine routing_id erstellt

alter table at_routing.urban_edges add column routing_id int GENERATED BY DEFAULT AS IDENTITY;

Danach wird die Tabelle für die Resultate des Routings erstellt:

create table at_routing.urban_routes (
    routing_id int,
    distance numeric,
    duration numeric
);

Für das Routing werde ich wieder eine Kombination aus OSRM und dem PostgreSQL HTTP Client verwenden. Eine genauere Anleitung findet sich im Beitrag Routing von OSM-Daten mit OSRM und Postgres. Die Abfrage selbst bietet wieder wenig Überraschungen

insert into at_routing.urban_routes
with select_routes as (
    select routing_id, ST_X(ST_Transform(geom, 4326)) || ',' || ST_Y(ST_Transform(geom, 4326)) || ';' || ST_X(ST_Transform(target_geom, 4326)) || ',' || ST_Y(ST_Transform(target_geom, 4326)) as route
from at_routing.urban_edges
)
select routing_id, (content::json->'routes'->0->>'distance')::numeric as distance, (content::json->'routes'->0->>'duration')::numeric as duration
from select_routes left join lateral http_get('http://127.0.0.1:5000/route/v1/driving/' || route || '?steps=false&overview=simplified&annotations=false') as b on true;

Prinzipiell lässt sich das Routing auch wieder mit einem simplen Python-Script parallelisieren. Siehe dazu auch Naive räumliche Interpolation in Postgres. Nach erfolgreichem Routing geht es an die Auswertung der Daten. Hierfür werden die Straßen wieder in einzelne Punkte zerlegt und mit den Resultaten verknüpft. Danach werden die einzelnen Punkte wieder zu Linien zusammengefügt. Jede Wegstrecke bekommt die Summe der zwei Punkte, die dann durch zwei dividiert wird, als Fahrzeit zugerechnet. Ist ein etwas simpler Ansatz, aber für einen ersten Test durchaus ausreichend. Das ganze speichere ich dann als materialized view.

create materialized view at_osm.v_routing_street_results as
with select_osm_points as (
    select osm_id, x.geom, x.path[1] as path_number from 
    at_osm.planet_osm_line left join lateral ST_DumpPoints(st_points(way)) as x on true 
    where highway in ('motorway', 'trunk', 'primary', 'secondary', 'tertiary', 'unclassified', 'residential', 
    'motorway_link', 'trunk_link', 'primary_link', 'secondary_link', 'tertiary_link', 'living_street')
),
get_min_duration as (
    select osm_id, target_geom, min(duration) / 60 as min_duration from at_routing.urban_routes as a, at_routing.urban_edges as b 
    where a.routing_id=b.routing_id
    group by 1,2
),
get_min_distance_per_point as (
    select a.osm_id, geom, path_number, min_duration from select_osm_points as a inner join 
    get_min_duration as b on (a.osm_id=b.osm_id and geom=target_geom)
)
select a.osm_id, ST_MakeLine(a.geom, b.geom) as geom, (a.min_duration + b.min_duration) / 2.0 as duration from 
get_min_distance_per_point as a inner join get_min_distance_per_point as b on (a.osm_id=b.osm_id and a.path_number+1=b.path_number)

Zusammenfassung

Im Unterschied zu Daten, die interpoliert werden sieht man am Beispiel, der hier dargestellten Erreichbarkeiten, schön die Topografie und die damit einhergehenden Herausforderungen. Während es im außeralpinen Raum ein dichtes Verkehrsnetz gibt, gibt es inneralpin naturgemäß weniger Verkehrswege. Ob weniger Verkehrswege automatisch etwas Schlechtes sind, muss jede(r) für sich selbst entscheiden (Stichwort Flächenversiegelung). Hier soll es einfach darum gehen die Möglichkeiten im Zusammenspiel mit Daten der OpenStreetMap und Postgres bzw. QGIS zu zeigen.