Finde Quellwasserpolygone

Dies ist eine Folgefrage zu dieser Frage .

Ich habe ein Flussnetz (mehrzeilig) und einige Entwässerungspolygone (siehe Bild unten). Mein Ziel ist es, nur die Quellwasserpolygone (grün) auszuwählen.

Mit Johns Lösung kann ich leicht die Flussstartpunkte (Sterne) extrahieren. Es kann jedoch Situationen geben (rotes Polygon), in denen ich Startpunkte in einem Polygon habe, das Polygon jedoch kein Quellwasserpolygon ist, da es vom Fluss geflogen wird. Ich möchte nur die Quellwasserpolygone.

Ich habe versucht, sie auszuwählen, indem ich die Anzahl der Schnittpunkte zwischen Polygonen und Flüssen gezählt habe (Begründung: Ein Quellwasserpolygon sollte nur einen Schnittpunkt mit dem Fluss haben).

SELECT 
    polyg.*
FROM 
    polyg, start_points, stream
WHERE 
    st_contains(polyg.geom, start_points.geom)
    AND ST_Npoints(ST_Intersection(poly.geom, stream.geom)) = 1

, wo poylg die poylgons sind, start_points von johns antworten und stream ist mein flussnetz .

Dies dauert jedoch ewig und ich habe es nicht ausgeführt:

"Nested Loop  (cost=0.00..20547115.26 rows=641247 width=3075)"
"  Join Filter: _st_contains(ezg.geom, start_points.geom)"
"  ->  Nested Loop  (cost=0.00..20264906.12 rows=327276 width=3075)"
"        Join Filter: (st_npoints(st_intersection(ezg.geom, rivers.geom)) = 1)"
"        ->  Seq Scan on ezg_2500km2_31467 ezg  (cost=0.00..2161.52 rows=1648 width=3075)"
"              Filter: ((st_area(geom) / 1000000::double precision) < 100::double precision)"
"        ->  Materialize  (cost=0.00..6364.77 rows=39718 width=318)"
"              ->  Seq Scan on stream_typ rivers  (cost=0.00..4498.18 rows=39718 width=318)"
"  ->  Index Scan using idx_river_starts on river_starts start_points  (cost=0.00..0.60 rows=1 width=32)"
"        Index Cond: (ezg.geom && geom)"

Meine Frage lautet also: Wie kann ich Quellwasserpolygone effizient abfragen?

Update: Ich habe meiner Dropbox einige Beispieldaten hinzugefügt . Die Daten stammen aus Südwestdeutschland. Es sind zwei Formdateien - eine mit Streams und eine mit Polygonen.

Ich glaube, der allgemeine Umriss (teilweise getestet) lautet:

1. Suchen Sie die Punkte, die Stream-Quellen darstellen, wie in dieser Antwort .

2. Schneiden Sie mit der Polygontabelle, um die Anzahl der Quellscheitelpunkte nach Polygon zu ermitteln.

3. Verwenden Sie ST_DumpPoints in Verbindung mit der Gruppierung nach Geometrie, um die Anzahl der einzelnen Punkte zu ermitteln. Die Idee war, zu zählen, wie viele Flüsse sich an einem bestimmten Punkt treffen.

Ein Beispiel für eine solche Abfrage:

SELECT count(geom), ST_AsText(geom) as wkt
FROM 
   (SELECT (ST_DumpPoints(foo.geom)).geom 
   FROM 
     (SELECT 
        ST_Collect(ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(10,10)),
                   ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(20,20))
        ) as geom
     ) foo 
 ) bar 
 GROUP BY geom; 

was zurückgibt:

count  |  wkt      
-------+--------------
 2     | POINT(0 0)
 1     | POINT(10 10)
 1     | POINT(20 20)

4. Führen Sie einen Schnittpunkt von 3gegen die Polygontabelle aus, um eine Anzahl (Summe der Eckpunkte) der Flussübergänge pro Polygon zu erhalten.

5. Verbinden Sie die Polygone von Anfang 2an 4und lehnen Sie diejenigen ab, bei denen die Anzahl (Summe der Eckpunkte) der Punkte an einer Kreuzung größer ist als die Summe der Flussquellen, die durch Summieren der Quellen durch Polygon aus den Schritten 1 und 2 erhalten werden. Wenn diese Bedingung zutrifft, gilt dies bedeutet, dass mindestens einer der Flüsse, die sich an einer Kreuzung treffen, außerhalb des betreffenden Polygons entstanden ist.

Diese können alle in einer größeren Folge von CTEs kombiniert werden, es sei denn, einige Tabellen wurden aus den Schritten mit Punkten erstellt (und indiziert).

Ich habe keine Ahnung, wie die Laufzeit für einen vollständigen Datensatz sein wird, da nur ein Teil davon für eine Teilmenge getestet wurde, aber mit einem räumlichen Index für die Polygontabelle wird es eine gewisse Unterstützung geben - dies ist offensichtlich nicht möglich Wenden Sie einen Index auf die Punkte an, die aus ST_DumpPoints hervorgehen. Daher ist dort ein vollständiger Scan erforderlich, obwohl sie sich bis dahin im Speicher befinden sollten.

Dies wird nicht als vollständige Antwort veröffentlicht , sondern als laufende Arbeit und als Chance, Fehler in der Logik zu finden. Arbeitsanfragen folgen in Kürze.

BEARBEITEN 1

Dies ist die Abfrage, die ich mir ausgedacht habe und die anscheinend für eine kleine Teilmenge Ihrer Daten funktioniert, die jedoch stundenlang für den gesamten Datensatz ausgeführt wird.

CREATE TABLE good_polys as  
   WITH 
     rivers as 
       (SELECT (ST_DUMP(ST_LineMerge(geom))).geom as geom FROM streams),
     start_points as
       (SELECT ST_StartPoint(geom) as geom FROM rivers),
     end_points as 
        (SELECT ST_EndPoint(geom) as geom FROM rivers),
     junctions as 
        (SELECT (ST_DumpPoints(geom)).geom 
        FROM (SELECT geom FROM streams) s),
     source_polygons as 
        (SELECT 
            count(rivers.geom) as source_count, 
            polygons.geom, 
            polygons.gid 
         FROM rivers, polygons
         WHERE st_intersects(polygons.geom, rivers.geom) 
         GROUP BY polygons.geom, polygons.gid),
     junction_polygons as 
        (SELECT 
            count(junctions.geom) as junction_count, 
            polygons.geom, 
            polygons.gid 
         FROM junctions, polygons
         WHERE st_intersects(polygons.geom, junctions.geom) 
         GROUP BY polygons.geom, polygons.gid)
    SELECT 
       jp.gid 
    FROM 
       junction_polygons jp, source_polygons sp 
    WHERE ST_Equals(jp.geom, sp.geom) 
    AND junction_count <= source_count;

BEARBEITEN 2 . Während dies für eine kleine Teilmenge korrekte Antworten zu liefern scheint, ist die Laufzeit für den gesamten Datensatz schrecklich , vermutlich weil die letzte Abfrage n ^ 2 Vergleiche durchführt und keinen räumlichen Index verwendet. Eine wahrscheinliche Lösung wäre, die Abfrage aufzuschlüsseln und Tabellen aus den Anfangspunkten und dem Punkt in Polygonabfragen zu erstellen, die dann vor dem letzten Schritt räumlich indiziert werden können.

Im Pseudocode sollte dies funktionieren:

• Wählen Sie alle aus Polygonen
• (FULL OUTER?) Mit Punkten auf Polygon verbinden schneidet Punkte
• (FULL OUTER?) Verbinden Sie Linien, bei denen das Polygon Linien schneidet
• waren line.riverid ist nicht gleich point.riverid
• Gruppe nach Polygonid
• count (pointid)> 0

Ich bin mir nicht sicher, wie ich die Abfrage erstellen soll, und ich kann sie nicht ohne eine Datenbank testen, auf der ich testen kann. Es ist eine ziemlich verrückte Frage, denke ich. Aber es sollte funktionieren!