Naive räumliche Interpolation in Postgres

Posted on So 10 Jänner 2021 in Blog

Im Artikel Routing von OSM-Daten mit OSRM und Postgres habe ich beschrieben, wie Routing auf Basis von Daten der OpenStreetMap funktionieren kann. Bei der Karte Erreichbarkeit von Oberzentren habe ich die TIN-Interpolation in QGIS verwendet. Jetzt habe ich mir gedacht, ich könnte das ganze auch mit einer ziemlich naiven Interpolation in SQL versuchen.

Inhalt:

Idee

Die Idee dieser naiven Interpolation ist es, die Regionalstatistischen Rastereinheiten der STATISTIK AUSTRIA heranzuziehen und diese Einheiten dann als Basis für die Interpolation zu verwenden. Denn auch bei der Interpolation in QGIS werden die Resultate als Rasterdatei gespeichert und den einzelnen Pixeln Werte zugewiesen. Natürlich ist dort die Auswahlmöglichkeit der Algorithmen, die diese Interpolation durchführen, besser und man sollte sie bei ernsthaften Analysen immer dem nachfolgenden Vorschlag vorziehen.

Routing

Um ein dichteres Netz an Punkten zu haben, habe ich im Gegensatz zur Beschreibung in Routing von OSM-Daten mit OSRM und Postgres noch den Typ village hinzugenommen. Obwohl das Routing auch mit Village auf meinem Server ganz flott durchgelaufen ist, habe ich zusätzlich noch folgendes kleines Python-Script geschrieben, das mir die Anfragen parallelisiert. Für die Verbindung wird das Service- und das Connection-File verwendet (siehe auch Connection Service File und pgpass.

import psycopg2
from multiprocessing import Pool

def route(modulo):
    connection = psycopg2.connect("service=osm")
    cursor = connection.cursor()
    query_string = "insert into routing.routes " \
        "with select_edge as ( " \
        "select routing_id, ST_X(ST_Transform(b.way, 4326)) || ',' || ST_Y(ST_Transform(b.way, 4326)) || ';' || " \
        "ST_X(ST_Transform(c.way, 4326)) || ',' || ST_Y(ST_Transform(c.way, 4326)) as route from " \
        "routing.edges as a, public.planet_osm_point as b, public.planet_osm_point as c " \
        "where b.osm_id=a.id1 and c.osm_id=a.id2 and routing_id % 10 = " + str(modulo) + " ) select routing_id, " \
        "(content::json->'routes'->0->>'distance')::numeric as distance, " \
        "(content::json->'routes'->0->>'duration')::numeric as duration from select_edge 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"

    cursor.execute(query_string,)
    connection.commit()
    cursor.close()
    connection.close()
    return True


if __name__ == '__main__':
    with Pool(10) as p:
        print(p.map(route, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]))

Shapes

Um meinen Ansatz umzusetzen, habe ich mir noch die Regionalstatistische Rastereinheiten vom Open Data-Portal der STATISTIK AUSTRIA besorgt. Von Eurostat kamen zusätzlich noch die Urban Audit 2020 Daten. Einspielen kann mit Dateien in Postgis mittels shp2pgsql

shp2pgsql -s 3035 STATISTIK_AUSTRIA_L001000_LAEA.shp statistik.regionalstatistische_rastereinheiten_001000_2016 | psql service=osm
shp2pgsql -s 3857 URAU_RG_100K_2020_3857_CITIES.shp eurostat.urban_audit_2020_cities | psql service=osm

SQL-Spaß

Das Herzstück des ganzen Versuchs ist natürlich die SQL-Abfrage, die die Eigenschaften von Postgis und Postgres nutzt. Hier holen wir uns im ersten Schritt alle Punkte, die innerhalb der City-Klassifikation des Urban Audit 2020 liegen:

with select_geom_ways as (
    select osm_id, way from eurostat.urban_audit_2020_cities as a, public.planet_osm_point as b 
    where cntr_code = 'AT' and place in ('town', 'city', 'village') and st_within(way, geom)
),

Danach holen wir uns alle Routen, die diese Punkte betreffen:

select_routes as (
    select id1, id2, duration  from routing.edges as a, routing.routes as b, select_geom_ways as c where a.routing_id=b.routing_id and id1=c.osm_id
    union 
    select id2, id1, duration from routing.edges as a, routing.routes as b, select_geom_ways as c where a.routing_id=b.routing_id and id2=c.osm_id
),

Daraufhin ermitteln wir für alle Punkte, die nicht über die City-Klassifikation ausgewählt wurden, die kürzeste Fahrtdauer zu einem der im ersten Schritt ausgewählten Punkte. Zusätzlich wird die Dauer für alle Punkte, die der City-Klassifikation entsprechen, auf 0 gesetzt.

select_overall as (
    select id2 as osm_id, way, min(duration / 60) as duration from select_routes as a, public.planet_osm_point as b where a.id2=b.osm_id and id2 not in (
    select osm_id from select_geom_ways
)
group by 1,2
union
select osm_id, way, 0::numeric as duration from select_geom_ways

Im nächsten Abschnitt transformieren wir die Rastereinheiten der STATISTIK AUSTRIA

select_raster as (
    select id, st_transform(geom, 3857) as geom from statistik.regionalstatistische_rastereinheiten_001000_2016
),

Danach wählen wir alle Polygone aus, die schon einen Punkt beinhalten und nehmen die kürzeste Dauer für diese Polygone an

select_osm_points as (
    select id, min(duration) as duration from select_raster as a,
    select_overall as b where ST_Within(b.way, a.geom)
    group by 1
)

Für alle Polygone der Rastereinheiten, die noch keine Punkte beinhalten, werden im Umkreis von 15.000 Metern Punkte ausgewählt. Dies geschieht jeweils vom Centroid der Polygone.

get_neighbours as (
    select id, geom as geom, ST_LENGTH(ST_MakeLine(ST_Centroid(geom), way)) as distance, duration, osm_id from
    select_raster as a join select_overall as b on (ST_DWithin(a.geom, b.way, 15000))
    where id not in (select id from select_osm_points)
)

Und im letzten Schritt wird dann die durchschnittliche Dauer bzw. der Median für jede Rastereinheit berechnet

select id, avg(duration) as duration, percentile_cont(0.5) WITHIN group (order by duration) as duration_median 
from get_neighbours group by id
UNION
SELECT id, duration, duration FROM select_osm_points;

Sehr viele der Parameter, die verwendet wurden, wurden mehr oder minder willkürlich gewählt. Beispielsweise hatte die Auswahl der 15.000 Meter den Hintergrund, die Anzahl der Leerfelder zu reduzieren. Um den Layer in QGIS zu laden, habe ich dann eine Ad-Hoc-Query verwendet. Das ist über die DB-Verwaltung in QGIS möglich:

select a.id, geom, duration, duration_median from statistik.regionalstatistische_rastereinheiten_001000_2016 as a, 
spielwiese.osm_ort_duration as b where a.id=b.id

Zusammenfassung

Die Umsetzung der Idee war nur durch Open Data möglich. So stellen Eurostat und die STATISTIK AUSTRIA sehr gute Datensätze zur Verfügung. Und das Straßennetz der OpenStreetMap ist sowieso sehr gut.

Ich habe mich beispielsweise für die 1000 Meter Rastereinheiten entschieden, da es ein guter Kompromiss zwischen Detailgrad und Rechenzeit war. Diese steigt natürlich, je kleiner die Rastereinheiten gewählt werden. Vom Routing her, wäre es natürlich auch möglich gewesen, die Centroide der Rastereinheiten zu den Städten des Urban Audits zu berechnen (mache ich unter Umständen noch, um zu sehen, wie sich die Resultate unterscheiden).

Ich persönlich sehe den Routing-Datensatz aber als Möglichkeit unterschiedliche Fragestellungen hinsichtlich Erreichbarkeiten mit nur einem Datensatz zu beantworten, ohne jedes Mal wieder einige Zeit in das Routing zu stecken. Denn so kann man mit einem Datensatz (sehr simpel) die relevanten Punkte auswählen.

Und wie man untenstehender Karte sieht, hat es trotzdem überraschenderweise gut funktioniert. Natürlich gibt es Rastereinheiten, die offensichtlich die falsche Farbe haben bzw. gar keine Farben, da kein Punkt des Routings innerhalb der 15.000 Meter erreichbar war.

Generell wirkt das Ergebnis nicht so verkehrt. Vor allem, wenn man die Einfachheit des Ansatzes bedenkt. Die ganze Abfrage ist auf meinem Server rund 10 Minuten gelaufen. Je nach Detailgrad der Interpolation in QGIS kann es wesentlich länger als die 10 Minuten dauern. Zusätzlich steigt dann auch die Dateigröße des generierten GeoTIFF, während die Größe der Resultate in der Datenbank durchaus überschaubar bleibt, da die räumliche Information in einer anderen Tabelle gespeichert wird und dadurch mehrfach verwendet werden kann.

Erreichbarkeiten von Städten in Österreich

Anhang

Da die SQL-Abfrage für die Erklärung in einzelne Stücke unterteilt wurde, gibt es sie hier noch als ganze Abfrage:

create table spielwiese.osm_ort_duration as
-- in OSM alle Punkte auswählen, die in die Urban Audit Klassifikation City fallen
with select_geom_ways as (
select osm_id, way from eurostat.urban_audit_2020_cities as a, public.planet_osm_point as b
where cntr_code = 'AT' and place in ('town', 'city', 'village') and st_within(way, geom)
),
-- alle routen von und zu diesen Punkten auswählen
select_routes as (
select id1, id2, duration  from routing.edges as a, routing.routes as b, select_geom_ways as c where a.routing_id=b.routing_id and id1=c.osm_id
union
select id2, id1, duration from routing.edges as a, routing.routes as b, select_geom_ways as c where a.routing_id=b.routing_id and id2=c.osm_id
),
-- für jeden Punkt die minimale Fahrtdauer auswählen, die Dauer in Minuten umwandlen und die "City"-Punkte auf 0 setzen
select_overall as (
select id2 as osm_id, way, min(duration / 60) as duration from select_routes as a, public.planet_osm_point as b where a.id2=b.osm_id and id2 not in (
select osm_id from select_geom_ways
)
group by 1,2
union
select osm_id, way, 0::numeric as duration from select_geom_ways
),
-- den Raster de Statistik Austria auswählen und nach EPSG:3857 transformieren
select_raster as (
select id, st_transform(geom, 3857) as geom from statistik.regionalstatistische_rastereinheiten_001000_2016
),
-- die Raster auswählen, die schon einen Punkt beinhalten
select_osm_points as (
select id, min(duration) as duration from select_raster as a,
select_overall as b where ST_Within(b.way, a.geom)
group by 1
),
-- für alle Raster, die keine Punkte beinhalten, die "Nachbarn" berechnen
get_neighbours as (
select id, geom as geom, ST_LENGTH(ST_MakeLine(ST_Centroid(geom), way)) as distance, duration, osm_id from
select_raster as a join select_overall as b on (ST_DWithin(a.geom, b.way, 15000))
where id not in (select id from select_osm_points)
)
select id, avg(duration) as duration, percentile_cont(0.5) WITHIN group (order by duration) as duration_median
from get_neighbours group by id
UNION
SELECT id, duration, duration FROM select_osm_points;