Encontrar en las cabeceras de los polígonos

Question

Esta es una pregunta de seguimiento a esta pregunta.

Tengo una red fluvial (multiline) y algunos de drenaje polígonos (ver imagen de abajo). Mi objetivo es seleccionar sólo las cabeceras de los polígonos (verde).

Con Juan solución me puede fácilmente extraer el río los puntos de inicio (estrellas). Sin embargo, puede haber situaciones (rojo polígono) donde he startpoints en un polígono, pero el polígono no es una cabecera poligonal, porque se volado aunque por el río. Yo sólo quiero las cabeceras de los polígonos.

Traté de seleccionar mediante el conteo del número de intersección entre los polígonos y los ríos (justificación: una cabecera de polígono debe tener sólo 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 poylgons, start_points de la universidad de johns respuesta y la secuencia es mi red fluvial.

Sin embargo, esta toma para siempre y yo no se ejecuta:

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

Así que mi pregunta es: ¿Cómo puedo eficiente de consulta en las cabeceras de los polígonos?

Actualización: He añadido algunos datos de ejemplo para mi dropbox. De datos es desde el sur-oeste de Alemania. Se trata de dos archivos shape - uno con arroyos y uno con polígonos.

Answer 1

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

1. Encontrar los puntos que representan los orígenes de secuencias, como en esta respuesta.

2. Se cruzan con los polígonos de la tabla para obtener un recuento de origen de los vértices del polígono.

3. Uso ST_DumpPoints en conjunto con el grupo de geometría para obtener un recuento de cada punto. Con la idea de obtener un recuento de cómo muchos ríos se juntan en un punto dado.

Un ejemplo de este tipo de 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)

4. Ejecutar una cruza de 3 contra el polígono de la tabla, para obtener un recuento (suma de los vértices) de río cruces por polígono.

5. Únete a los polígonos de 2 a 4, rechazando aquellos donde la cuenta (suma de los vértices) de los puntos de unión es mayor que la suma de las fuentes del río, obtenida sumando las fuentes por el polígono de los pasos 1 y 2. Si se cumple esta condición, esto significa que al menos uno de los ríos de la reunión, en un cruce, se originó fuera del polígono en cuestión.

Todos estos pueden ser combinados en un largish secuencia de Cte, a menos que algunas de las tablas fueron creadas a partir de los pasos que involucran puntos (y el índice).

No tengo idea de lo que el tiempo de ejecución de este será en un conjunto de datos completo, después de haber probado solo una parte de este en un subconjunto, pero con un índice espacial en los polígonos de la tabla, habrá un poco de ayuda -- es obvio que no es posible aplicar un índice de los puntos que surgen de ST_DumpPoints, por lo que un análisis completo será necesario, aunque deben estar en la memoria.

Esto no es de ser publicado como una respuesta completa, pero como un trabajo en progreso, y una oportunidad para encontrar fallos en la lógica. Trabajando consultas próximamente.

EDICIÓN 1

Esta es la consulta que se me ocurrió, que aparece a trabajar en un pequeño subconjunto de los datos, pero funciona durante horas en la totalidad del conjunto de datos.

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;

EDICIÓN 2. Mientras que esto parece producir respuestas correctas en un pequeño subconjunto, el tiempo de ejecución en el conjunto de datos completo es horrible, presumiblemente debido a la final de la consulta está haciendo n^2 comparaciones y no el uso de un índice espacial. Probable solución sería romper consulta hacia abajo y crear tablas a partir de los puntos inicial y el punto en el polígono de las consultas, que luego pueden ser espacialmente indexado antes, el paso final.

Answer 2

En pseudo-código, esto debería funcionar:

• seleccione todos los de polígonos
• (FULL OUTER?) unir con puntos del polígono cruza puntos
• (FULL OUTER?) líneas de unión polígono donde se cruza con las líneas de
• fueron de la línea.riverid no es igual que el punto.riverid
• grupo por polygonid
• count (pointid) > 0

No estoy realmente seguro de cómo construir la consulta, y yo no puedo probarlo sin una base de datos para hacer pruebas. Es bastante loco consulta, creo. Pero debería funcionar!