Publicerad den Lämna en kommentar

Mera PostGIS ”st” magi

I förra veckan sneglade jag lite på detaljplaner och PostGIS och ett av de grundläggande problemen var att använda så få dataset till så mycket informationslager som möjligt, utan att informationen krockar och överlagras. Mycket information i planerna är linjer, men det är i många fall smidigare att registrera dessa som polygoner. I detta avsnitt tänkte jag mera systematiskt titta på vilka möjligheter vi har med olika ”st” kommandon i PostGIS.

För mina tester skapar jag en väldigt enkel polygon delad i ett rutnät, där tabellen i huvudsak har två numeriska fält (typ_a och typ_b).

Skärmbild från 2018-09-08 12-55-23.png

I QGIS delar jag upp värdet i de båda fälten så att jag kan göra indelningar i såväl typ_a som typ_b på olika sätt.

Skärmklipp från 2018-09-08 13:08:32.png

Om jag vill rita ut den omgivande kantlinjen runt alla ingående, men sammanhängande polygoner så går det inte med QGIS stilar. Men ett virtuellt lager skulle fungera, vilket använder sql, men beräkningen görs i QGIS. Genom att använda DB Manager så kan jag skicka ett sql kommando till PostGIS databasen där beräkningen görs och jag får bara resultatet. Det första är smidigt, speciellt med lokala data, men om det exempelvis är en begränsning i nätverkets överföringshastighet så är DB Manager bättre. Många kommandon stöds heller inte i virtuella lager och i många fall går det snabbare på en dedikerad databasserver. Jag kommer därför enbart att presentera kommandon på detta vis, men mycket av det som presenteras går som sagt att skapa även med virtuella lager direkt mot lokala lager i QGIS.

Skärmklipp från 2018-09-08 13:15:55

I bilden ovan har jag skapat ett enkelt SQL-kommando i DB Manager (klicka på andra knappen från vänster i DB Manager).

SELECT st_union(geom)
From "Detaljplan"."test_tabell";

Detta slår samman alla geometrier i tabellen ”test_tabell” i schemat ”Detaljplan”. Resultatet laddas in i QGIS som nytt lager och här behöver man först välja geometrikolumn och ge lagret ett namn innan det går att läsa in resultatet. Rutan stängs sedan inte, så nu kan jag dels spara frågan till senare tillfällen, eller modifiera den och köra igen.

Om jag inte vill ha en polygon, utan endast en gränslinje så modifierar jag uttrycket något.

SELECT st_boundary(st_union(geom))
FROM "Detaljplan"."test_tabell";

Skärmklipp från 2018-09-08 13:22:45.png

I många fall är det just kantlinjen jag är ute efter och detta ger mig just det resultatet.

Om jag i stället vill ha kantlinjerna runt alla områden indelade efter ”typ_a” så blir uttrycket lite längre.

SELECT st_boundary(st_union(geom))
FROM "Detaljplan"."test_tabell"
GROUP BY "typ_a";

Skärmklipp från 2018-09-08 13:26:28.png

Detta gör samma sak som tidigare, men grupperar först alla objekt efter värdet i kolumnen ”typ_a”. Här introduceras två problem.

  1. Kantlinjen runt hela polygonen är med, vilket kanske inte är önskvärt.
  2. Det är dubbla linjer mellan områdena, vilket gör en del stilar svåra att använda på ett snyggt sätt.

Vi tar problem 1 först. Det kan fixas med lite ”algebra”. Yttergränsen är ju redan beräknad, så det borde gå att ”ta bort” den från det nya resultatet. Det går att göra med st_difference.

SELECT st_difference(typ_a.geom, hela.geom)
FROM
(SELECT st_boundary(st_union(geom)) AS geom
FROM "Detaljplan"."test_tabell") AS hela,
(SELECT st_boundary(st_union(geom)) AS geom
FROM "Detaljplan"."test_tabell"
GROUP BY "typ_a") AS typ_a;

Det nya SQL kommandot är i grunden ett st_difference som tar två geometrier som argument. Dessa är här angivna som två alias som sedan definieras i FROM med varsin parentes där de ”döps” med AS.

Varje SELECT i FROM satsen måste dessutom ges en alias för att st_difference skall kunna hämta resultatet. Det görs genom AS geom för båda.

St_difference tar nu resultatet med alias typ_a och presenterar det som skiljer sig från resultatet med alias hela.

Men vi har fortfarande problem 2 att fixa, eftersom det är många sammanfallande linjer i resultatet (vilket inte direkt syns i bilden ovan).

Jag har i uttrycken nedan dessutom valt att lägga till en alias för den utvalda geometrin med AS geometri för att göra min egen hantering lite enklare, men det är inte nödvändigt så länge man är noga med att kontrollera att man använder rätt geometrikolumn när man lägger till sina lager från DB Manager.

SELECT st_unaryunion(st_collect(st_difference(typ_a.geom, hela.geom))) as geometri
FROM
(SELECT st_boundary(st_union(geom)) AS geom
FROM "Detaljplan"."test_tabell") AS hela,
(SELECT st_boundary(st_union(geom)) AS geom
FROM "Detaljplan"."test_tabell"
GROUP BY "typ_a") AS typ_a;

Skärmklipp från 2018-09-08 13:57:17.png

I uttrycket är två saker tillagda. St_collect som tar alla geometrier och returnerar en ”samling” med geometrier, som sedan st_unaryunion sammanfogar till ett enda linjeobjekt utan överlappande geometrier.

Detta introducerar dock ett nytt problem, nämligen att det är ett multipart-objekt, så stilsättningen ser lite konstig ut på sina ställen med långa linjer eller oregelbundna intervall. Detta kan i och för sig lösas med en geometrigenerator i QGIS (line_merge()), men det går även att göra ytterligare ett tillägg till SQL-kommandot.

SELECT st_linemerge(st_unaryunion(st_collect(st_difference(typ_a.geom, hela.geom)))) as geometri
FROM
(SELECT st_boundary(st_union(geom)) AS geom
FROM "Detaljplan"."test_tabell") AS hela,
(SELECT st_boundary(st_union(geom)) AS geom
FROM "Detaljplan"."test_tabell"
GROUP BY "typ_a") AS typ_a;

Skärmklipp från 2018-09-08 14:03:31.png

Den enda skillnaden mot tidigare är att det dessutom finns ett st_linemerge() kommando som omsluter urvalet.

Om nu detta var allt jag ville få fram så borde jag vara klar. Men det är nu det börjar krångla till sig på riktigt.

Jag kanske dessutom vill ha ett lager med gränser även för typ_b, där dessa inte sammanfaller med gränser i ytterkant eller mellan typ_a områden.

Jag vill ha:

  • Ytterlinjer för sammanfogade typ_b områden,
  • som inte sammanfaller med ytterlinjer för typ_a områden,
  • eller yttergräns.

Den sista punkten kan jag bortse ifrån här då jag kan använda hela typ_a gränserna och inte bara de som ligger mellan typ_a områden. Sedan borde det gå att skapa ett nytt uttryck liknande det ovan, om än lite krångligare.

Ett lager motsvarande typ_a fast för typ_b är inte så svårt att skapa, men för att belysa problemet ytterligare så redigerar jag värdena något i databasen.

Skärmklipp från 2018-09-08 14:39:18.png

Resultatet ovan är från en kopia från tidigare kommandon för typ_a men där referenser ändrats till att grupperas efter värdet i fält typ_b. Ytterlinjerna är borta, men jag skulle vilja att även linjerna mellan typ_a områden skulle försvinna… Vilket visat sig vara lite krångligare att få till.

SELECT st_linemerge(st_unaryunion(st_collect(st_difference(typ_b.geom, typ_a.geom)))) as geometri
FROM
(SELECT st_linemerge(st_unaryunion(st_collect(aa.geom))) AS geom
FROM
(SELECT st_boundary(st_union(geom)) AS geom
FROM "Detaljplan"."test_tabell"
GROUP BY "typ_a") AS aa) AS typ_a
,
(SELECT st_linemerge(st_unaryunion(st_collect(bb.geom))) AS geom
FROM
(SELECT st_boundary(st_union(geom)) AS geom
FROM "Detaljplan"."test_tabell"
GROUP BY "typ_b") AS bb) AS typ_b;

Skärmklipp från 2018-09-08 15:04:33.png

Jag tvivlar på att det måste bli så krångligt som mitt SQL-kommando är ovan, men jag var tvungen att förenkla de båda gränserna för typ_a och typ_b så att de var singelgeometrier och endast ett objekt per geometri. Detta gjorde jag genom att använda de st kommandon som jag använt tidigare och bädda in dessa i varandra i två steg. Först därefter lyckades jag isolera endast de unika gränserna för typ_b.

Nu kan jag däremot skapa lager för enskilda polygontyper i QGIS och med SQL kommandon hämta enskilda unika gränser som inte överlappar varandra, utan att ha skapat ytterligare tabeller i databasen.

Skärmklipp från 2018-09-08 15:11:43.png

Ändrar jag något i ursprungstabellen, geometri eller fältvärde, så behöver kartan bara uppdateras så kommer alla lager att anpassas till de nya förutsättningarna.

Sedan kan jag kombinera detta med andra SQL lager, eller inbyggda funktioner i QGIS för att ytterligare bygga lager som ändå bara har ett enda polygonlager som källa, så länge geometrierna sammanfaller med det jag är ute efter.

Jag kanske bara vill ha gränser typ_b i ytor av typ_a men en annan markering av ytor med värdet ”3,1”. Då kan jag slänga in en WHERE i SQL kommandot för typ_b gränser och ett enkelt virtuellt lager för ”3,1” ytorna är lätt att fixa.

Skärmbild från 2018-09-08 15-23-37.png

Om man inte vill ha gränser så kan man bara utelämna st_boundary i kommandona. Ja det kan bli lite krångligare än så eftersom man byter geometri, men kommandot totalt lär bara bli enklare. Väldigt mycket av dessa polygoner går dessutom ganska enkelt att skapa med virtuella lager i stället för med SQL.

select typ_a, st_union(geometry) from test_tabell group by typ_a;

Skärmklipp från 2018-09-08 15:55:06.png

Om du som här bygger upp ett projekt med flera lager som hämtar information från samma källa, så kan det hjälpa att ställa in databeroenden.

Skärmbild från 2018-09-09 13-22-18.png

Ett SQL lager kanske skall uppdateras automatiskt om det ursprungliga ”rådatalagret” redigeras. För varje lager kan man i lagerinställningarna ställa in så att QGIS håller koll på dessa lager och automatiskt uppdaterar lagret om det sker en förändring.

Sent tillägg

Om du har data för flera områden i en och samma tabell, men vill kunna välja ut endast de data som finns exempelvis i ett specifikt område, då kan du använda fler ”st” kommandon.

Data är här lagrat i tabellen ”test_tabell” och för att testa så skapas ett polygonlager till vid namn ”region”.

Skärmklipp från 2018-09-12 18:26:40.png

Det går att skapa virtuella lager där data från test_tabell (B) som har en geografisk relation till polygonen i regiontabellen (A).

SELECT B.* FROM A, B WHERE st_contains(A, B)
SELECT B.* FROM A, B WHERE st_intersects(A, B)
SELECT B.* FROM A, B WHERE st_overlaps(A, B)

Resultatet från st_contains() är alla polygoner som helt ryms innanför regionen. För st_intersects() är det alla objekt som med någon del rör vid regionen, och för st_overlaps() är det objekt som delvis täcks av regionen.

Det går även att baka in dessa kommandon i befintliga SQL kommandon direkt mot PostGIS.

Skärmklipp från 2018-09-12 18:40:43

Ovanstående är ett enkelt exempel från artikeln, medan övriga uttryck blir något mer omfattande. Det går dock att modifiera alla SQL-lager så att de anpassas efter formen på ett ”regionlager”.

Skärmklipp från 2018-09-12 19:18:03.png

Detta tillägg är främst med anledning av en fråga om just geografiska urval i PostGIS jag fick från Norge nyligen, och det passade bra att lägga in detta i denna sedan tidigare planerade artikel.

Avslutning

Detta inlägg har hjälpt mig att få lite bättre ordning på några av alla st kommandon som finns för databaser. Det är mycket kvar innan det är helt solklart, men varje gång jag sätter mig ner och samtidigt skriver ner lite så blir det tydligare (och enklare att återskapa).

Jag hoppas att du som läser detta också tycker att det blev lite tydligare och att du kan ha användning av dessa kommandon.

Nyhet från Geosupportsystem , orginal inlägg

Lämna ett svar