Trouver des polygones d'amont

8

Il s'agit d'une question complémentaire à cette question .

J'ai un réseau fluvial (multiligne) et quelques polygones de drainage (voir photo ci-dessous). Mon objectif est de ne sélectionner que les polygones d'amont (vert).

entrez la description de l'image ici

Avec la solution de John, je peux facilement extraire les points de départ de la rivière (étoiles). Cependant, je peux avoir des situations (polygone rouge) où j'ai des points de départ dans un polygone, mais le polygone n'est pas un polygone d'amont, car il est survolé par la rivière. Je veux seulement les polygones d'amont.

J'ai essayé de les sélectionner en comptant le nombre d'intersection entre les polygones et les rivières (justification: un polygone d'amont ne devrait avoir qu'une seule intersection avec la rivière)

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

, où poylg sont les poylgons, start_points de johns answer et stream est mon réseau fluvial.

Cependant, cela prend une éternité et je ne l'ai pas exécuté:

"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)"

Ma question est donc la suivante: comment puis-je interroger efficacement les polygones d'amont?

Mise à jour: j'ai ajouté des exemples de données à ma boîte de dépôt . Les données proviennent du sud-ouest de l'Allemagne. Il s'agit de deux fichiers de formes - un avec des flux et un avec des polygones.

EDi
la source
Donc, pour être clair, vous voulez les polygones qui ne contiennent que des points de départ, pas les points de départ eux-mêmes. Et les points de départ sont définis comme dans votre question précédente (à laquelle j'ai répondu, et pour autant que je sache), correctement?
John Powell
Jupp, seuls les polygones qui contiennent des points de départ ET ne sont pas passés par une rivière / ne sont que des départs de la rivière. Le polygone rouge ci-dessus contient des points de départ, mais n'est PAS un polygone d'amont car la rivière le traverse / ne démarre pas à l'intérieur du polygone ...
EDi
Donc, vous voulez que l'ensemble de polygonsceux-ci ne contienne que les points qui sont des sources fluviales (de la question précédente) et exclure tout point de rencontre de deux fleuves. Désolé, pour toutes les questions, je veux juste être sûr.
John Powell
Non, par exemple, dans le polygone vert inférieur, deux rivières se rencontrent également. Je veux exclure ceux polygonsqui ont une rivière qui passe (la rivière entre et sort du polygone) et garder ceux avec des départs (et les rivières ne quittent que ce polygone).
EDi
1
Je ne connais pas de PostGIS, donc je ne peux pas m'empêcher avec du code direct, cependant, dans ArcGIS, je suivrais ces lignes: (1) intersecter les lignes et les polygones dans un fichier de points. (2) supprimer (spatialement) des points identiques. (3) ajoutez un champ numérique au paramètre de point avec la valeur 1 pour chaque point. (4) joindre spatialement le polygone aux points et en utilisant la somme du champ numérique pour indiquer le type de drainage. Une somme de 1 signifie qu'il s'agit d'un promontoire. Supérieur à 1 signifie qu'il y a plus d'une entrée ou sortie.
Mikkel Lydholm Rasmussen

Réponses:

4

Je pense que le schéma général (partiellement testé jusqu'à présent) est le suivant:

  1. Trouvez les points représentant les sources de flux, comme dans cette réponse .

  2. Intersection avec la table des polygones pour obtenir un nombre de sommets source par polygone.

  3. Utilisez ST_DumpPoints conjointement avec le regroupement par géométrie pour obtenir un décompte de chaque point. L'idée étant de compter le nombre de rivières qui se rencontrent à un point donné.

Un exemple d'une telle requête:

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; 

qui renvoie:

count  |  wkt      
-------+--------------
 2     | POINT(0 0)
 1     | POINT(10 10)
 1     | POINT(20 20)
  1. Exécutez une intersection de 3contre la table des polygones, pour obtenir un compte (somme des sommets) des jonctions de rivière par polygone.

  2. Joignez les polygones à partir 2de 4, en rejetant ceux où le nombre (somme de sommets) de points à une jonction est supérieur à la somme des sources de la rivière, obtenu en sommant les sources par polygone des étapes 1 et 2. Si cette condition est vérifiée, elle signifie qu'au moins une des rivières se rencontrant à une jonction, provient de l'extérieur du polygone en question.

Ceux-ci peuvent tous être combinés ensemble dans une séquence étendue de CTE, à moins que certaines tables aient été créées à partir des étapes impliquant des points (et indexées).

Je n'ai aucune idée de ce que sera l'exécution de cela sur un ensemble de données complet, n'ayant testé qu'une partie de cela sur un sous-ensemble, mais avec un index spatial sur la table des polygones, il y aura une aide - il n'est évidemment pas possible de appliquer un index aux points qui émergent de ST_DumpPoints, donc une analyse complète y sera requise, bien qu'ils devraient être en mémoire d'ici là.

Ce n'est pas affiché comme une réponse complète , mais comme un travail en cours et une chance de trouver des failles dans la logique. Requêtes de travail à venir.

EDIT 1

C'est la requête que j'ai trouvée, qui semble fonctionner sur un petit sous-ensemble de vos données, mais qui s'exécute pendant des heures sur l'ensemble de données complet.

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;

MODIFIER 2 . Bien que cela semble produire des réponses correctes sur un petit sous-ensemble, le temps d'exécution sur l'ensemble de données complet est horrible , probablement parce que la requête finale effectue n ^ 2 comparaisons et n'utilise pas d'index spatial. La solution probable serait de décomposer la requête et de créer des tables à partir des points initiaux et des points dans les requêtes polygonales, qui peuvent ensuite être indexées spatialement avant l'étape finale.

John Powell
la source
La requête est en cours d'exécution sur mon bureau. Je n'ai aucune idée du temps que cela prendra ou si ce sera correct, bien que cela semble raisonnable à partir d'un petit échantillon de vos données. Avez-vous une idée du nombre de polygones qui répondent à vos critères?
John Powell
Je vais exécuter la requête sur un serveur. Je pense que seule une petite partie des polygones répondra aux critères de sélection ...
EDi
C'est ce que j'ai trouvé sur un sous-ensemble. Je posterai ma requête une fois qu'elle sera terminée
John Powell
Simplification demain.
John Powell
Désolé, très occupé aujourd'hui. Je pense que la réponse est d'exécuter la requête source et les requêtes de jonctions fluviales en premier, d'intersecter la table des polygones pour obtenir les nombres par polygone, de les enregistrer en tant que tables, puis de les indexer. Exécutez ensuite la dernière étape, où les géométries sont égales, et comparez le nombre de points des deux tables. J'espère que cela utilisera alors un index plutôt que de faire des comparaisons n² comme présentes. Je reviendrai plus tard.
John Powell
0

En pseudo code, cela devrait fonctionner:

  • tout sélectionner parmi les polygones
  • (FULL OUTER?) Joindre avec des points sur un polygone intersecte des points
  • (FULL OUTER?) Joindre les lignes où le polygone coupe les lignes
  • étaient line.riverid n'est pas égal à point.riverid
  • grouper par polygonide
  • count (pointid)> 0

Je ne sais pas vraiment comment créer la requête, et je ne peux pas la tester sans base de données sur laquelle tester. C'est une question assez folle, je pense. Mais ça devrait marcher!

Alex Leith
la source