Separe los polígonos basados ​​en la intersección usando PostGIS

36

Tengo una tabla de polígonos PostGIS donde algunos se cruzan entre sí. Esto es lo que estoy tratando de hacer:

  • Para un polígono determinado seleccionado por id, dame todos los polígonos que se cruzan. Básicamente,select the_geom from the_table where ST_Intersects(the_geom, (select the_geom from the_table where source_id = '123'))
  • A partir de estos polígonos, necesito crear nuevos polígonos para que la intersección se convierta en un nuevo polígono. Entonces, si el polígono A se cruza con el polígono B, obtendré 3 nuevos polígonos: A menos AB, AB y B menos AB.

¿Algunas ideas?

atogle
fuente
1
Hmmm, creo que veo a qué te refieres, pero un simple gráfico podría hacer maravillas para ayudarme (y a otros) a ver exactamente lo que quieres.
Jason
Agregué algo en mi respuesta.
Adam Matan

Respuestas:

29

Como usted dijo que obtiene un grupo de polígonos de intersección para cada polígono que le interesa, puede crear lo que se conoce como "superposición de polígonos".

Esto no es exactamente lo que está haciendo la solución de Adam. Para ver la diferencia, eche un vistazo a esta imagen de una intersección ABC:

Intersección ABC

Creo que la solución de Adam creará un polígono "AB" que cubre el área de "AB! C" y "ABC", así como un polígono "AC" que cubre "AC! B" y "ABC", y un " BC "polígono que es" BC! A "y" ABC ". Por lo tanto, los polígonos de salida "AB", "AC" y "BC" se superpondrían al área "ABC".

Una superposición de polígonos produce polígonos no superpuestos, por lo que AB! C sería un polígono y ABC sería un polígono.

Crear una superposición de polígonos en PostGIS es bastante sencillo.

Básicamente hay tres pasos.

El paso 1 es extraer la línea [Tenga en cuenta que estoy usando el anillo exterior del polígono, se vuelve un poco más complicado si desea manejar correctamente los agujeros]:

SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines

El paso 2 es "nodear" la línea de trabajo (producir un nodo en cada intersección). Algunas bibliotecas como JTS tienen clases "Noder" que puede usar para hacer esto, pero en PostGIS la función ST_Union lo hace por usted:

SELECT ST_Union(the_geom) AS the_geom FROM (...your lines...) AS noded_lines

El paso 3 es crear todos los polígonos no superpuestos posibles que pueden provenir de todas esas líneas, realizado por la función ST_Polygonize :

SELECT ST_Polygonize(the_geom) AS the_geom FROM (...your noded lines...)

Puede guardar la salida de cada uno de esos pasos en una tabla temporal, o puede combinarlos en una sola declaración:

CREATE TABLE my_poly_overlay AS
SELECT geom FROM ST_Dump((
    SELECT ST_Polygonize(the_geom) AS the_geom FROM (
        SELECT ST_Union(the_geom) AS the_geom FROM (
            SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines
        ) AS noded_lines
    )
)

Estoy usando ST_Dump porque la salida de ST_Polygonize es una colección de geometría, y (generalmente) es más conveniente tener una tabla donde cada fila sea uno de los polígonos que componen la superposición de polígonos.

Jeff
fuente
Tenga en cuenta que ST_ExteriorRingcae cualquier agujero. ST_Boundarypreservará los anillos interiores, pero también creará un polígono dentro de ellos.
jpmc26
La unión de los anillos exteriores se bloquea cuando hay demasiados polígonos. Desafortunadamente, esta solución "sencilla" solo funciona para coberturas pequeñas.
Pierre Racine
13

Si entiendo correctamente, quieres tomar (A es la geometría izquierda, B es la derecha):

Imagen de A∪B http://img838.imageshack.us/img838/3996/intersectab1.png

Y extraer:

A ∖ AB

Imagen de A ∖ AB http://img830.imageshack.us/img830/273/intersectab2.png

AB

Imagen de AB http://img828.imageshack.us/img828/7413/intersectab3.png

y B ∖ AB

Imagen de B ∖ AB http://img839.imageshack.us/img839/5458/intersectab4.png

Es decir, tres geometrías diferentes para cada par de intersección.

Primero, creemos una vista de todas las geometrías de intersección. Asumiendo que el nombre de su tabla es polygons_table, usaremos:

CREATE OR REPLACE VIEW p_intersections AS    -- Create a view with the 
SELECT t1.the_geom as t1_geom,               -- intersecting geoms. Each pair
       t2.the_geom as t2_geom                -- appears once (t2.id<t2.id)
    FROM polygons_table t1, polygons_table t2  
         WHERE t1.id<t2.id AND t1.the_geom && t2.the_geom 
                           AND intersects t1.the_geom, t2.the_geom;

Ahora tenemos una vista (prácticamente, una tabla de solo lectura) que almacena pares de gemas que se cruzan, donde cada par aparece solo una vez debido a t1.id<t2.id condición.

Ahora reunamos sus intersecciones A∖AB, ABy B∖AB, usando SQL UNIONen las tres consultas:

--AB:
SELECT ST_intersection(t1.the_geom, t2.the_geom) 
    AS geom 
    FROM p_intersections

UNION 

--AAB:
SELECT ST_Difference(t1.the_geom, t2.the_geom) 
    AS geom 
    FROM p_intersections

UNION

--BAB:
SELECT ST_Difference(t2.the_geom, t1.the_geom) 
    AS geom 
    FROM p_intersections;

Notas:

  1. El &&operador se utiliza como filtro antes que el intersectsoperador, para mejorar el rendimiento.
  2. Elegí crear una VIEWconsulta gigantesca en lugar de una gigantesca; Esto es solo por conveniencia.
  3. Si quiso decir ABes la unión, no la intersección, de Ay B- Use ST_Union en lugar de st_intersection en la UNIONconsulta en los lugares apropiados.
  4. los signo es un signo unicode para establecer diferencia; elimínelo del código si confunde su base de datos.
  5. Fotografías cortesía de la agradable categoría de teoría de conjuntos de Wikimedia Commons .
Adam Matan
fuente
Mi boleto de error en meta: meta.gis.stackexchange.com/questions/79/…
Adam Matan
Buena explicación! Los resultados son los mismos que en la solución scw, pero su camino debería ser más rápido (no calcula / ni almacena / intersecciones adicionales de A y B)
stachu
¡Gracias! Creo que no almaceno ninguna información adicional, ya que solo creo VIEW SQL, no tablas.
Adam Matan
Sí, eso es cierto, pero usted calcula Intersección adicional de A y B, lo cual no es necesario
stachu
55
Las imágenes en esta respuesta ya no funcionan.
Fezter
8

Lo que está describiendo es la forma en que el operador de Union trabaja en ArcGIS, pero es un poco diferente de Union o Intersection en el mundo GEOS. El manual de Shapely tiene ejemplos de cómo funcionan los conjuntos en GEOS . Sin embargo, el wiki de PostGIS tiene un buen ejemplo con el trabajo de línea que debería ser el truco para usted.

Alternativamente, podría calcular tres cosas:

  1. ST_Intersection (geom_a, geom_b)
  2. ST_Difference (geom_a, geom_b)
  3. ST_Difference (geom_b, geom_a)

Esos deberían ser los tres polígonos que mencionaste en tu segundo punto.

scw
fuente
2
El ejemplo wiki de PostGIS es bueno
fmark
¿ST_Intersects no devolvería un valor booleano si se cruzan o no? Creo que ST_Intersection funcionaría.
Jason
Sí, error tipográfico de mi parte - arreglado en el original ahora, ¡gracias Jason!
scw
-2

Algo como:

INSERTAR EN NUEVOS VALORES DE TABLA ((seleccione id, the_geom de old_table donde st_intersects (the_geom, (seleccione the_geom de old_table donde id = '123')) = true

EDITAR: necesita la intersección real del polígono.

INSERTAR EN los valores de new_table ((seleccione id, ST_Intersection (the_geom, (seleccione the_geom de antiguo donde id = 123))

A ver si eso funciona.

George Silva
fuente