Open Data mit Postgres und Postgis (Teil 2)
Posted on Mo 23 Dezember 2024 in Blog
In Teil 1 wurden schon erste Verschneidungen vorgenommen und alle Rastereinheiten, die innerhalb einer Gemeinde liegen, zugeordnet. Im zweiten Teil soll nun eine Möglichkeit gesucht werden, wie man die nicht überschneidungsfreien Rastereinheiten Gemeinden zuordnet.
Zusätzlich sieht man, wie Postgis sich automatisch kleiner Optimierungstricks behilft, indem automatisch der entsprechende Index herangezogen wird.
Inhalt:
Zuordnung der Rastereinheiten nach Anteil zu Gemeinden
Eine Idee ist es, die Rastereinheiten zu teilen und die Bevölkerung dann der Gemeinde zuzuordnen, die den größten Anteil an der Rastereinheit hat. Man wird wohl nie eine perfekte Aufteilung, im Sinne von "die Aufteilung entspricht 1:1 der Bevölkerung zum 1.1.2024" hinbekommen, da man hinsichtlich der Rastereinheiten, die sich auf mehrere Gemeinden aufteilen, immer Annahmen treffen muss.
Sicherlich könnte man noch zusätzliche Datensätze hinzuziehen, um eine bessere Aufteilung vornehmen zu können. Beispiel dafür wären eventuell Informationen zu Gebäuden aus der OpenStreetMap oder auch aus dem Adressregister des BEV. Aber selbst wenn man Informationen zu den Gebäuden hat, weiß man nicht, wie sich die Bevölkerung in dieser Rasterzelle auf die Gebäude verteilt. Ein Beispiel dafür ist die Grenze zwischen Bad Gleichenberg und Kapfenstein: https://www.openstreetmap.org/#map=19/46.891132/15.936136.
Um die Zellen aufzuteilen, wird man mit st_intersects, st_area und st_difference arbeiten.
explain analyze
with order_cells as (
select
rr.name as grid_name
, rr.id as grid_id
, gag.reg_code
, row_number() over (partition by rr.name order by st_area(st_difference(rr.geom, st_transform(gag.geom, 3035))) / st_area(rr.geom) desc) as grid_id
from dwh_oesterreich.gemeinden_aktueller_gebietsstand gag,
dwh_inspire.at_bev_rastereinheiten_2024 rr
where st_intersects(rr.geom, st_transform(gag.geom, 3035))
)
select * from order_cells where grid_id = 1
Mittels window function werden die Zellen nach Größe geordnet und es wird dann jeweils die Rastereinheit der Gemeinde zugeordnet, die den größten Anteil an der Zelle hat. Das ganze geht dank Index auch wieder ziemlich flott, wie der Output zeigt:
Subquery Scan on order_cells (cost=23058077.01..54645857.39 rows=6218 width=45) (actual time=708465.380..710897.689 rows=587945 loops=1)
Filter: (order_cells.grid_id = 1)
Rows Removed by Filter: 24482
-> WindowAgg (cost=23058077.01..54630313.59 rows=1243504 width=53) (actual time=708465.378..710849.335 rows=612427 loops=1)
-> Gather Merge (cost=23058077.01..23206967.51 rows=1243504 width=45) (actual time=708465.355..710543.201 rows=612427 loops=1)
Workers Planned: 4
Workers Launched: 4
-> Sort (cost=23057076.95..23057854.14 rows=310876 width=45) (actual time=707921.248..707937.767 rows=122485 loops=5)
Sort Key: rr.name, ((st_area(st_difference(rr.geom, st_transform(g.geom, 3035), '-1'::double precision)) / st_area(rr.geom))) DESC
Sort Method: quicksort Memory: 12641kB
Worker 0: Sort Method: quicksort Memory: 12590kB
Worker 1: Sort Method: quicksort Memory: 12643kB
Worker 2: Sort Method: quicksort Memory: 12666kB
Worker 3: Sort Method: quicksort Memory: 12668kB
-> Nested Loop (cost=297822.53..23028715.77 rows=310876 width=45) (actual time=2931.919..707095.721 rows=122485 loops=5)
-> Parallel Hash Join (cost=297822.39..372034.61 rows=146986 width=687) (actual time=2869.211..3364.099 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..90.923 rows=117589 loops=5)
-> Parallel Hash (cost=238842.95..238842.95 rows=1686995 width=159) (actual time=1736.383..1736.384 rows=1686991 loops=5)
Buckets: 2097152 Batches: 8 Memory Usage: 215136kB
-> Parallel Seq Scan on regionalstatistische_rastereinheiten_000100_2016 rr (cost=0.00..238842.95 rows=1686995 width=159) (actual time=0.026..613.457 rows=1686991 loops=5)
-> Index Scan using gemeinden_2024_geom_transform_idx on gemeinden_2024 g (cost=0.14..25.16 rows=1 width=28963) (actual time=2.826..3.789 rows=1 loops=587945)
Index Cond: (st_transform(geom, 3035) && rr.geom)
Filter: st_intersects(rr.geom, st_transform(geom, 3035))
Rows Removed by Filter: 1
Planning Time: 0.403 ms
Execution Time: 710919.148 ms
Was fällt am Query-Plan auf? Obwohl nur st_intersects verwendet wird, nutzt die Abfrage den Index. Zuerst wird im Index überprüft, ob die bounding boxen sich schneiden (siehe: Index Cond: (st_transform(geom, 3035) && rr.geom) und dann wird nach den Kombinationen gefiltert, wo sich die beiden Geometrien schneiden.
Jetzt kann man die Abfrage aus Teil 1 mit der obigen Abfrage kombinieren und so einen Gesamtdatensatz für Österreich erstellen.
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))
),
order_cells as (
select
rr.name as grid_name
, rr.id as grid_id
, gag.reg_code
, rr.value
, row_number() over (partition by rr.name order by st_area(st_difference(rr.geom, st_transform(gag.geom, 3035))) / st_area(rr.geom) desc) as grid_id
from dwh_oesterreich.gemeinden_aktueller_gebietsstand gag,
dwh_inspire.at_bev_rastereinheiten_2024 rr
where st_intersects(rr.geom, st_transform(gag.geom, 3035))
)
select grid_name, grid_id, reg_code, value as bev from select_matches
union
select grid_name, grid_id, reg_code, value as bev from order_cells where grid_id = 1
Natürlich passen die Resultate nicht 1:1 mit den Einwohnerzahlen zusammen, die von der STATISTIK AUSTRIA veröffentlicht werden, da eine Aufteilung der Rastereinheiten nicht eindeutig sein kann. Wie so eine Abweichung aussehen kann, sieht man am Beispiel von Graz. Mit der gewählten Methode der Verschneidungen kommt man für Graz auf 302.539 Einwohner:innen. Laut offiziellen Zahlen hat Graz zum 1.1.2024 eine Bevölkerung von 302.749 Personen. Wie ich finde, durchaus eine gute Näherung.
Fazit und Ausblick
Für ganz Österreich kommen wir auf die Summe von 9.158.750 Einwohner:innen. Das schaut mir sehr stark nach der Zahl aus, die wir uns in Teil 1 erwartet hätten.
Hier sieht man, wie man sich mit den Standardfunktionen von Postgis sehr gut behelfen kann, um räumliche Probleme zu lösen. Im nächsten Teil wird man dann sehen, wie man die Informationen aus dem Breitbandatlas nutzen kann, um über die verfügbaren Internetgeschwindigkeiten in einer Gemeinde Auskunft zu geben.