Open Data mit Postgres und Postgis (Teil 1)

Posted on Mo 16 Dezember 2024 in Blog

In einer mehrteiligen Serie will ich Schritt für Schritt dem Ziel mehrere Datensätze miteinander zu verschneiden, um neue Informationen zu generieren, näher kommen.

Das große Ziel ist es, die Bevölkerung auf 100x100 Meter Rastereinheiten Gemeinden zuordnen, um herauszufinden, wie viel Prozent der Bevölkerung theoretisch 1000 Mbit-Internet zur Verfügung haben.

Vor einiger Zeit habe ich schon einmal einen Beitrag zu den Inspire Bevölkerungsdaten geschrieben. Jetzt geht es darum, unter anderem diesen Datensatz zu nutzen, um neue Fragestellungen beantworten zu können.

Inhalt:

Vorbereitung

Bei den Daten, die für diese Blogserie verwendet werden, handelt es sich um vier Open-Data-Datensätze:

Die Daten habe ich in die Datenbank eingespielt und die Tabellen folgendermaßen genannt:

  • data_gv_at.breitbandatlas
  • inspire.at_bev_rastereinheiten_2024
  • oesterreich.gemeinden_2024
  • oesterreich.regionalstatistische_rastereinheiten_000100_2016

Da die Daten in unterschiedlichen CRS vorliegen, wird bei der Tabelle oesterreich.gemeinden_2024 auch noch ein Index erstellt, der die Daten umprojiziert:

CREATE INDEX gemeinden_2024_geom_transform_idx ON oesterreich.gemeinden_2024 USING gist (st_transform(geom, 3035));

Zusätzlich habe ich noch in zwei Schemata, die mittels dbt bespielt werden, zwei Views erstellt.

Viewdefinition dwh_oesterreich.gemeinden_aktueller_gebietsstand

CREATE OR REPLACE VIEW dwh_oesterreich.gemeinden_aktueller_gebietsstand
AS SELECT g.g_id AS reg_code,
    g.g_name AS reg_name,
    g.geom
   FROM oesterreich.gemeinden_2024 g;

Viewdefinition dwh_inspire.at_bev_rastereinheiten_2024

CREATE OR REPLACE VIEW dwh_inspire.at_bev_rastereinheiten_2024
AS WITH unnest_data AS (
         SELECT s.statisticalvalue_value AS value,
            replace(s.statisticalvalue_dimensions_dimensions_spatial_href::text, 'https://data.inspire.gv.at/77679c2c-302c-11e3-beb4-0000c1ab0db6/su.StatisticalGridCell/AT_'::text, ''::text) AS spatial_id
           FROM inspire.at_bev_rastereinheiten_2024 s
        )
 SELECT rr.gid,
    rr.id,
    rr.name,
    rr.geom,
    ud.value
   FROM oesterreich.regionalstatistische_rastereinheiten_000100_2016 rr,
    unnest_data ud
  WHERE rr.name::text = ud.spatial_id;

So weit ist alles einmal vorbereitet, daher kann man sich an die Umsetzung machen.

Erste Zuordnung der Rastereinheiten zu Gemeinden

In einem ersten Versuch, die Daten zu verschneiden, sollen in einem ersten Schritt die Bevölkerung den Gemeinden zugeordnet werden. Ein erster Ansatz dafür ist den @ Operator zu verwenden:

explain analyze
select
    rr.name as grid_name
    , rr.id as grid_id
    , gag.reg_code
from dwh_oesterreich.gemeinden_aktueller_gebietsstand gag,
dwh_inspire.at_bev_rastereinheiten_2024 rr
where
rr.geom @ st_transform(gag.geom, 3035);

Das Ergebnis der Ausführung:

Nested Loop  (cost=298822.53..539489.08 rows=1243504 width=37) (actual time=3911.313..14447.485 rows=1081508 loops=1)
  ->  Gather  (cost=298822.39..431829.11 rows=587945 width=687) (actual time=3910.324..4979.164 rows=587945 loops=1)
        Workers Planned: 4
        Workers Launched: 4
        ->  Parallel Hash Join  (cost=297822.39..372034.61 rows=146986 width=687) (actual time=3622.743..4001.742 rows=117589 loops=5)
              Hash Cond: (replace((s.statisticalvalue_dimensions_dimensions_spatial_href)::text, 'https://data.inspire.gv.at/77679c2c-302c-11e3-beb4-0000c1ab0db6/su.StatisticalGridCell/AT_'::text, ''::text) = (rr.name)::text)
              ->  Parallel Seq Scan on at_bev_rastereinheiten_2024 s  (cost=0.00..29467.86 rows=146986 width=121) (actual time=0.293..167.364 rows=117589 loops=5)
              ->  Parallel Hash  (cost=238842.95..238842.95 rows=1686995 width=159) (actual time=2400.825..2400.826 rows=1686991 loops=5)
                    Buckets: 2097152  Batches: 8  Memory Usage: 215104kB
                    ->  Parallel Seq Scan on regionalstatistische_rastereinheiten_000100_2016 rr  (cost=0.00..238842.95 rows=1686995 width=159) (actual time=0.149..1306.069 rows=1686991 loops=5)
  ->  Index Scan using gemeinden_2024_geom_transform_idx on gemeinden_2024 g  (cost=0.14..0.16 rows=1 width=28963) (actual time=0.014..0.016 rows=2 loops=587945)
        Index Cond: (st_transform(geom, 3035) ~ rr.geom)
Planning Time: 48.913 ms
Execution Time: 14485.325 ms

Ganz flott unterwegs. Wie gehofft wird auch der oben angelegte Index für eine schnellere Ausführung genutzt.

Allerdings ist das Resultat nicht wirklich brauchbar, denn wenn man sich die Summe der Bevölkerung ansieht, dann erkennt man den Fehler:

with bev as (
    select
        rr.name as grid_name
        , rr.id as grid_id
        , gag.reg_code
        , rr.value
    from dwh_oesterreich.gemeinden_aktueller_gebietsstand gag,
    dwh_inspire.at_bev_rastereinheiten_2024 rr
    where
    rr.geom @ st_transform(gag.geom, 3035)
)
select sum(value) from bev;

sum
16.456.577

Österreich ist entweder sehr stark gewachsen oder es gibt Probleme mit der Zuordnung. Und in diesem Fall liegt es an der Zuordnung und am @ - Operator. Dieser vergleicht Bounding Boxen miteinander und wenn man sich die Dokumentation ansieht, dann erkennt man das Problem:

The @ operator returns TRUE if the bounding box of geometry A is completely contained by the bounding box of geometry B.

Bei den Raster-Einheiten sind die Bounding Boxen noch sehr einfach, aber bei komplexeren Geometrien, wie es bei den Gemeindegebieten der Fall ist, handelt es sich bei Bounding Boxen um eine Annäherung, die dazu führt, dass viele Rastereinheiten, die nicht im Gemeindegebiet sind, ausgewählt werden (siehe auch die Erklärung im GIS Garden: Bounding Box).

Durch die Nutzung des Index ist die Abfrage allerdings sehr flott. Aber wie funktioniert ein räumlicher Index? Ein Blick in die Dokumentation von Postgis zahlt sich hier aus:

Spatial indexes are a little different – they are unable to index the geometric features themselves and instead index the bounding boxes of the features.

Daher müssen wir die Abfrage anpassen. Allerdings nehmen wir die Bounding Boxen, um den Index trotzdem nutzen zu können und die Abfrage zu beschleunigen.

Danach filtern wir die Resultate noch nach den Zellen, die wirklich innerhalb der Gemeinde liegen. Dafür fragen wir ab, ob die Geometrie innerhalb der Gemeinde liegt.

explain analyze
select
    rr.name as grid_name
    , rr.id as grid_id
    , gag.reg_code
from dwh_oesterreich.gemeinden_aktueller_gebietsstand gag,
dwh_inspire.at_bev_rastereinheiten_2024 rr
where
rr.geom @ st_transform(gag.geom, 3035) and st_within(rr.geom, st_transform(gag.geom, 3035));

Wie man an der Ausgabe von explain analyze sieht, ist die Dauer der Abfrage etwas länger: 14 Sekunden vs. über 7 Minuten. Hier sollte man aber nicht das Augenmerk auf die Dauer alleine legen, sondern auch auf das Resultat.

Zumal die Dauer natürlich auch immer abhängig ist von der verwendeten Postgres-Version, der eingesetzten Hardware usw. und soll daher hier nur zeigen, bei der Dauer kommt es auch drauf an, welche Genauigkeit man für eine Aufgabe benötigt.

Aber die längere Dauer ist verständlich, da mehrere räumliche Checks durchgeführt werden müssen. In einem ersten Schritt wird nur die Menge der Geometrien, die verglichen werden müssen, verringert, aber es bleiben trotzdem noch viele Vergleiche über:

Gather  (cost=298822.53..15180913.84 rows=1244 width=37) (actual time=2639.986..435999.956 rows=563685 loops=1)
    Workers Planned: 4
    Workers Launched: 4
    ->  Nested Loop  (cost=297822.53..15179789.44 rows=311 width=37) (actual time=2626.555..435024.269 rows=112737 loops=5)
          ->  Parallel Hash Join  (cost=297822.39..372034.61 rows=146986 width=687) (actual time=2561.304..3016.059 rows=117589 loops=5)
                Hash Cond: (replace((s.statisticalvalue_dimensions_dimensions_spatial_href)::text, 'https://data.inspire.gv.at/77679c2c-302c-11e3-beb4-0000c1ab0db6/su.StatisticalGridCell/AT_'::text, ''::text) = (rr.name)::text)
                ->  Parallel Seq Scan on at_bev_rastereinheiten_2024 s  (cost=0.00..29467.86 rows=146986 width=121) (actual time=0.015..84.071 rows=117589 loops=5)
                ->  Parallel Hash  (cost=238842.95..238842.95 rows=1686995 width=159) (actual time=1648.108..1648.109 rows=1686991 loops=5)
                      Buckets: 2097152  Batches: 8  Memory Usage: 215168kB
                      ->  Parallel Seq Scan on regionalstatistische_rastereinheiten_000100_2016 rr  (cost=0.00..238842.95 rows=1686995 width=159) (actual time=0.022..586.139 rows=1686991 loops=5)
          ->  Index Scan using gemeinden_2024_geom_transform_idx on gemeinden_2024 g  (cost=0.14..25.17 rows=1 width=28963) (actual time=2.886..3.673 rows=1 loops=587945)
                Index Cond: ((st_transform(geom, 3035) ~ rr.geom) AND (st_transform(geom, 3035) ~ rr.geom))
                Filter: st_within(rr.geom, st_transform(geom, 3035))
                Rows Removed by Filter: 1
  Planning Time: 47.381 ms
  Execution Time: 436024.527 ms

Sind damit alle relevanten Rastereinheiten einer Gemeinde zugeordnet? Dafür sollten wir uns die Bevölkerung ansehen:

with select_matches as (
    select
        rr.name as grid_name
        , rr.id as grid_id
        , gag.reg_code
        , rr.value
    from dwh_oesterreich.gemeinden_aktueller_gebietsstand gag,
    dwh_inspire.at_bev_rastereinheiten_2024 rr
    where
    rr.geom @ st_transform(gag.geom, 3035) and st_within(rr.geom, st_transform(gag.geom, 3035))
)
select sum(value) from select_matches

sum
8.838.911

Da fehlen also noch Personen, denn die Bevölkerung nach 100m ETRS-LAEA-Raster beinhaltet 9.158.750 Personen, was der offiziellen Bevölkerung zum 1.1.2024 entspricht.

Woran liegt die Differenz?

Es gibt natürlich auch bei 100x100 Metern Rastereinheiten, durch die direkt die Gemeindegrenzen verlaufen, daher fehlen in diesem Schritt noch einige Zellen, die man sich noch gesondert ansehen muss.

Fazit und Ausblick

In einem ersten Teil hat man gesehen, wie man mit einer simplen Abfrage die Rastereinheiten den einzelnen Gemeinden zuordnen kann und was man dabei beachten muss. Man konnte schon einen Großteil der Bevölkerung den Gemeinden zuweisen, aber es gibt noch Verbesserungsmöglichkeiten.

Wie man diese Zuordnung noch verbessern kann und welche Herausforderungen dabei warten, sieht man dann im zweiten Teil der Serie.

In dritten und letzten Teil sind man dann, wie man die Ergebnisse kombinieren und mit anderen Datensätzen verknüpfen kann: Open Data mit Postgres und Postgis (Teil 3)

Postgres GIS Postgis Open Data 2024