Durchführen eines Polygon-Drilldowns / Overlays in PostGIS?

8

Ich stehe mit PostGIS vor einer Herausforderung, bei der ich meinen Kopf nicht herumwickeln kann. Ich weiß, dass ich dies mit einer Programmiersprache lösen kann (und das ist mein Backup-Plan), aber ich mag es wirklich, dies in PostGIS zu lösen. Ich habe versucht zu suchen, konnte aber keine Antworten finden, die meinem Problem entsprechen. Dies könnte daran liegen, dass ich bei meinen Suchbegriffen unsicher bin. Entschuldigen Sie dies bitte und weisen Sie mich in die richtige Richtung, da tatsächlich eine Antwort vorliegt.

Mein Problem ist folgendes:

  • Ich habe eine Tabelle mit gemischten Polygonen / MultiPolygonen
  • Jedes (Multi) Polygon hat ein Attribut, das es ordnet (Priorität)
  • Jedes Polygon hat auch einen Wert, den ich gerne wissen würde
  • Ich habe einen Suchbereich (Polygon)
  • Für meinen Abfragebereich möchte ich den Bereich finden, der von jedem Polygonwert abgedeckt wird

Beispiel:

Angenommen, ich habe die drei Polygone hier in Rot, Grün und Indigo dargestellt: Beispiel

Und dass das kleinere blaue Rechteck mein Abfragepolygon ist

Darüber hinaus sind die Attribute

geom   | rank | value
-------|------|----  
red    |  3   | 0.1
green  |  1   | 0.2
indigo |  2   | 0.2

Ich möchte diese Geometrien so auswählen, dass der höchste Rang (grün) den gesamten Bereich ausfüllt, den er kann (dh den Schnittpunkt zwischen meinem Abfrage-Geom und diesem Geom), und der nächsthöhere (Indigo) den Schnittpunkt zwischen dem Abfrage-Geom ausfüllt und das geom MINUS das bereits abgedeckte) etc.

Etwas wie das: Geben Sie hier die Bildbeschreibung ein

Ich habe diese Frage gefunden: Verwenden Sie ST_Difference, um überlappende Features zu entfernen? aber es scheint nicht zu tun, was ich will.

Ich kann selbst herausfinden, wie man Bereiche und dergleichen berechnet, daher ist eine Abfrage, die mir die drei Geometrien gibt, wie im zweiten Bild dargestellt, in Ordnung!

Zusätzliche Informationen: - Dies ist keine große Tabelle (~ 2000 Zeilen) - Es kann keine oder mehrere Überlappungen geben (nicht nur drei) - Es gibt möglicherweise keine Polygone in meinem Abfragebereich (oder nur in Teilen davon) - i ' Ich laufe Postgis 2.3 auf Postgres 9.6.6

Meine Fallback-Lösung besteht darin, eine Abfrage wie folgt durchzuführen:

SELECT 
ST_Intersection(geom, querygeom) as intersection, rank, value
FROM mytable
WHERE ST_Intersects(geom, querygeom)
ORDER by rank asc

Und dann iterativ Teile der Geometrien im Code "abhacken". Aber wie gesagt, ich würde das wirklich gerne in PostGIS machen

atlefren
quelle
2
WITH RECURSIVE ...
Ich kann Ihnen momentan
1
Oh und überprüfen Sie dies
Geozelot
Vielen Dank! Ich werde es morgen versuchen, wenn sich sonst niemand gezwungen fühlt, eine Komplettlösung anzubieten.
Atlefren

Antworten:

7

Ich denke das funktioniert.

Es handelt sich um eine Fensterfunktion, die den Unterschied zwischen dem Schnittpunkt jeder Geometrie mit dem Abfragefeld und der Vereinigung der vorhergehenden Geometrien ermittelt.

Die Koaleszenz wird benötigt, da die Vereinigung der vorhergehenden Geometrien für die erste Geometrie null ist, was anstelle des gewünschten Ergebnisses ein Nullergebnis ergibt.

WITH a(geom, rank, val) AS
(
    VALUES
    ('POLYGON((1 1, 1 5, 5 5, 5 1, 1 1))'::geometry,3,0.1),
    ('POLYGON((2 3, 2 8, 4 8, 5 3,2 3))'::geometry,1,0.2),
    ('POLYGON((4 4, 4 6, 6 6, 6 4,4 4))'::geometry,2,0.2)
)
,q AS
(
    SELECT 'POLYGON((3 3, 3 4.5, 12 4.5,12 3, 3 3))'::geometry geom
) 
SELECT 
  ST_Difference(
    ST_Intersection(a.geom, q.geom), 
    COALESCE(ST_Union(a.geom) 
           OVER (ORDER BY rank ROWS BETWEEN UNBOUNDED PRECEDING and 1 PRECEDING),
       'POLYGON EMPTY'::geometry)
  ) geom 
FROM a,q
WHERE ST_Intersects(a.geom, q.geom);

Ich bin mir allerdings nicht sicher, wie es funktioniert. Da jedoch sowohl ST_Union als auch ST_Intersection als unveränderlich markiert sind, ist dies möglicherweise nicht so schlimm.

Nicklas Avén
quelle
Das hat wie ein Zauber gewirkt! Sie müssen Ihre Abfrage nur in eine andere Abfrage
einschließen
5

Ein etwas anderer Ansatz. Es gibt eine Einschränkung, bei der ich nicht weiß, wie die Leistung skaliert werden soll, aber bei einer indizierten Tabelle sollte dies in Ordnung sein. Es funktioniert ungefähr genauso wie Nicklas 'Abfrage (ein bisschen langsamer?), Aber die Messung an einer so kleinen Stichprobe ist schwierig.

Es sieht viel hässlicher aus als die Abfrage von Nicklas, vermeidet jedoch eine Rekursion in der Abfrage.

WITH 
    a(geom, rank, val) AS
    (
        VALUES
        ('POLYGON((1 1, 1 5, 5 5, 5 1, 1 1))'::geometry,3,0.1),
        ('POLYGON((2 3, 2 8, 4 8, 5 3, 2 3))'::geometry,1,0.2),
        ('POLYGON((4 4, 4 6, 6 6, 6 4, 4 4))'::geometry,2,0.2)
    ),
    polygonized AS (
        -- get the basic building block polygons
        SELECT (ST_Dump(         -- 5. Dump out the polygons
            ST_Polygonize(line)) -- 4. Polygonise the linework
            ).geom AS mypoly
        FROM (
            SELECT 
                ST_Node(                  -- 3. Node lines on intersection
                    ST_Union(             -- 2. Union them for noding
                        ST_Boundary(geom) -- 1. Change to linestrings
                    ) 
                ) 
                AS line
            FROM a
        ) line
    ),
    overlayed AS ( 
        -- Overlay with original polygons and get minimum rank.
        -- Union and dissolve interior boundaries for like ranks.
        SELECT (ST_Dump(ST_UnaryUnion(geom))).geom, rank 
        FROM ( 
            -- union up the polygons by rank, unary union doesn't count as an aggregate function?
            SELECT ST_Union(mypoly) geom, rank 
            FROM ( 
                -- get the minimum rank for each of the polygons
                SELECT p.mypoly, min(rank) rank
                FROM polygonized p 
                    INNER JOIN a ON ST_Contains(a.geom,ST_PointOnSurface(p.mypoly))
                GROUP BY p.mypoly
                ) g
            GROUP BY rank
            ) r
    )
-- get the intersection of the query area with the overlayed polygons
SELECT ST_Intersection(o.geom,'POLYGON((3 3, 3 4.5, 12 4.5,12 3, 3 3))'::Geometry), rank
FROM overlayed o
WHERE ST_Intersects(o.geom,'POLYGON((3 3, 3 4.5, 12 4.5,12 3, 3 3))'::Geometry) and
    -- Intersection can do funky things
    GeometryType(ST_Intersection(o.geom,'POLYGON((3 3, 3 4.5, 12 4.5,12 3, 3 3))'::Geometry)) like '%POLYGON';
MickyT
quelle
1

Da ich darüber WITH RECURSIVEgeredet habe, werde ich damit eine schnelle und schmutzige Antwort hinzufügen.

Dies ist ungefähr so ​​gut wie die Lösung von @ NicklasAvén für drei Polygone, die beim Hochskalieren nicht getestet werden konnte.
Nach dem Stand beider Lösungen hat diese einen kleinen Vorteil gegenüber der anderen: Wenn zum Beispiel das Polygon mit Rang = 2 in dem von Rang = 1 enthalten ist , werden die ...WHERE GeometryType = 'POLYGON'Filter herausgefiltert, während es sonst ein gibt GEOMETRYCOLLECTION EMPTY(ich habe die Geometrie geändert des jeweiligen Polygons in meiner Lösung entsprechend, um ein Beispiel zu geben; dies gilt auch für andere Fälle, in denen kein Schnittpunkt mit der Differenz gefunden wird). Dies ist jedoch leicht in den anderen Lösungen enthalten und möglicherweise nicht einmal von Belang.

WITH RECURSIVE
    a(geom, rank, val) AS (
        VALUES
           ('POLYGON((1 1, 1 5, 5 5, 5 1, 1 1))'::geometry,3,0.1),
           ('POLYGON((2 3, 2 8, 4 8, 5 3,2 3))'::geometry,1,0.2),
           ('POLYGON((2.1 3.1, 2.1 7.9, 3.9 7.9, 4.9 3.1,2.1 3.1))'::geometry,2,0.2)
    ),
    q AS (
        SELECT 'POLYGON((3 3, 3 4.5, 12 4.5,12 3, 3 3))'::geometry geom
    ),
    src AS (
        SELECT ROW_NUMBER() OVER(ORDER BY rank) AS rn,
               ST_Intersection(q.geom, a.geom) AS geom,
               rank,
               val
        FROM a
        JOIN q
           ON ST_Intersects(a.geom, q.geom)
    ),
    res AS (
        SELECT s.geom AS its,
               ST_Difference(q.geom, s.geom) AS dif,
               s.rank,
               s.val,
               2 AS icr
        FROM src AS s,
             q
        WHERE s.rn = 1
        UNION ALL
        SELECT ST_Intersection(s.geom, r.dif) AS its,
               ST_Difference(r.dif, s.geom) AS dif,
               s.rank,
               s.val,
               icr + 1 AS icr
        FROM src AS s,
             res AS r
        WHERE s.rank = r.icr
    )

SELECT its AS geom,
       rank,
       val
FROM res
WHERE GeometryType(its) = 'POLYGON'
Geozelot
quelle