Encuentra polígonos de cabecera

8

Esta es una pregunta de seguimiento a esta pregunta .

Tengo una red fluvial (multilínea) y algunos polígonos de drenaje (ver imagen a continuación). Mi objetivo es seleccionar solo los polígonos de cabecera (verde).

ingrese la descripción de la imagen aquí

Con la solución de John puedo extraer fácilmente los puntos de inicio del río (estrellas). Sin embargo, puedo tener situaciones (polígono rojo) en las que tengo puntos de inicio en un polígono, pero el polígono no es un polígono de cabecera, porque el río lo vuela. Solo quiero los polígonos de cabecera.

Traté de seleccionarlos contando el número de intersección entre polígonos y ríos (justificación: un polígono de cabecera debería tener solo 1 intersección con el río)

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

, donde poylg son los polígonos, start_points de johns answer y stream es mi red fluvial.

Sin embargo, esto lleva una eternidad y no lo ejecuté:

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

Entonces mi pregunta es: ¿cómo puedo consultar eficientemente los polígonos de cabecera?

Actualización: agregué algunos datos de muestra a mi Dropbox . Los datos son del suroeste de Alemania. Son dos archivos de forma: uno con secuencias y otro con polígonos.

EDi
fuente
Entonces, para ser claros, desea los polígonos que solo contienen puntos de inicio, no los puntos de inicio en sí. ¿Y los puntos de inicio se definen como en su pregunta anterior (que respondí, y que yo sepa), correctamente?
John Powell el
Jupp, solo los polígonos que contienen puntos de inicio Y no son pasados ​​por un río / son solo comienzos del río. El polígono rojo anterior contiene puntos de inicio, pero NO es un polígono de cabecera ya que el río fluye a través de él / no comienza dentro del polígono ...
EDi
Por lo tanto, desea que el conjunto de polygonsesos contenga solo aquellos puntos que son fuentes de ríos (de la pregunta anterior) y excluya cualquier lugar donde se unan dos ríos. Lo sentimos, para todas las preguntas, solo quiero estar seguro.
John Powell el
No, por ejemplo, en el polígono verde inferior también se encuentran dos ríos. Quiero excluir a aquellos polygonsque tienen un río que pasa (el río entra y sale del polígono) y mantener a aquellos con comienzos (y los ríos solo dejan este polígono).
EDi
1
No conozco ningún PostGIS, por lo que no puedo ayudar con el código directo, sin embargo, en ArcGIS, iría a lo largo de estas líneas: (1) se cruzan entre líneas y polígonos en un archivo de puntos. (2) eliminar (espacialmente) puntos idénticos. (3) agregue un campo numérico al parámetro de punto con el valor de 1 para cada punto. (4) unir espacialmente el polígono en los puntos y usar la suma del campo numérico para indicar el tipo de drenaje. Una suma de 1 significa que es un promontorio. Superior a 1 significa que hay más de 1 entrada o salida.
Mikkel Lydholm Rasmussen

Respuestas:

4

Creo que el esquema general (parcialmente probado hasta ahora) es:

  1. Encuentre los puntos que representan las fuentes de flujo, como en esta respuesta .

  2. Interseca con la tabla de polígonos para obtener un recuento de vértices de origen por polígono.

  3. Use ST_DumpPoints junto con group by geometry para obtener un recuento de cada punto. La idea es contar cuántos ríos se encuentran en un punto dado.

Un ejemplo de tal consulta:

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; 

que devuelve:

count  |  wkt      
-------+--------------
 2     | POINT(0 0)
 1     | POINT(10 10)
 1     | POINT(20 20)
  1. Ejecute intersecciones de 3contra la tabla de polígonos, para obtener un recuento (suma de vértices) de uniones de ríos por polígono.

  2. Une los polígonos desde 2adelante 4, rechazando aquellos donde el recuento (suma de vértices) de puntos en una unión es mayor que la suma de las fuentes del río, obtenido al sumar las fuentes por polígono de los pasos 1 y 2. Si esta condición se cumple, significa que al menos uno de los ríos que se unen en un cruce, se originó fuera del polígono en cuestión.

Todos estos se pueden combinar en una secuencia grande de CTE, a menos que se hayan creado algunas tablas a partir de los pasos que involucran puntos (e indexados).

No tengo idea de cuál será el tiempo de ejecución de esto en un conjunto de datos completo, habiendo probado solo una parte de esto en un subconjunto, pero con un índice espacial en la tabla de polígonos, habrá algo de ayuda, obviamente no es posible aplique un índice a los puntos que emergen de ST_DumpPoints, por lo que se requerirá un escaneo completo allí, aunque para entonces deberían estar en la memoria.

Esto no se publica como una respuesta completa , sino como un trabajo en progreso y una oportunidad para encontrar fallas en la lógica. Consultas de trabajo próximamente.

EDITAR 1

Esta es la consulta que se me ocurrió, que parece funcionar en un pequeño subconjunto de sus datos, pero se ejecuta durante horas en el conjunto de datos completo.

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;

EDITAR 2 . Si bien esto parece producir respuestas correctas en un pequeño subconjunto, el tiempo de ejecución en el conjunto de datos completo es horrible , presumiblemente porque la consulta final está haciendo n ^ 2 comparaciones y no está usando un índice espacial. La solución probable sería dividir la consulta y crear tablas a partir de los puntos iniciales y el punto en consultas de polígonos, que luego se pueden indexar espacialmente antes del paso final.

John Powell
fuente
La consulta se está ejecutando actualmente en mi escritorio. No tengo idea de cuánto tiempo llevará o si será correcto, aunque parecía razonable a partir de una pequeña muestra de sus datos. ¿Tienes idea de cuántos polígonos cumplen con tus criterios?
John Powell
Ejecutaré la consulta en un servidor. Creo que solo una pequeña parte de los polígonos cumplirá con los criterios de selección ...
EDi
Eso es lo que encontré en un subconjunto. Publicaré mi consulta una vez que termine
John Powell
Simplificación mañana.
John Powell
Lo siento, muy ocupado hoy. Creo que la respuesta es ejecutar primero la consulta de origen y las consultas de cruce de ríos, intersectarse con la tabla de polígonos para obtener recuentos por polígono, guardarlos como tablas y luego indexarlos. Luego ejecute el paso final, donde las geometrías son iguales, y compare el recuento de puntos de las dos tablas. Espero que esto use un índice en lugar de hacer n² comparaciones como presente. Volveré a publicar más tarde.
John Powell
0

En pseudocódigo, esto debería funcionar:

  • seleccionar todo de polígonos
  • (¿COMPLETO EXTERIOR?) Unirse con puntos en puntos de intersección de polígonos
  • (¿COMPLETO EXTERIOR?) Unir líneas donde el polígono intersecta líneas
  • donde line.riverid no es igual a point.riverid
  • grupo por polígono
  • cuenta (pointid)> 0

No estoy realmente seguro de cómo construir la consulta, y no puedo probarla sin una base de datos para probar. Es una consulta bastante loca, creo. ¡Pero debería funcionar!

Alex Leith
fuente