Separación de polígonos en función de la intersección mediante PostGIS

Question

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

• Para un polígono 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 un nuevo polígono tal que la intersección se convierta en un nuevo polígono. Así, 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.

¿Alguna idea?

Answer 1

Como has dicho que obtienes un grupo de polígonos de intersección para cada polígono que te interesa, es posible que quieras crear lo que se conoce como "superposición de polígonos".

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

Creo que la solución de Adam creará un polígono "AB" que cubre tanto el área de "AB!C" como de "ABC", así como un polígono "AC" que cubre "AC!B" y "ABC", y un polígono "BC" que es "BC!A" y "ABC". Así, los polígonos de salida "AB", "AC" y "BC" se solaparían con el á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.

La creación de una superposición de polígonos en PostGIS es en realidad bastante sencilla.

Hay básicamente tres pasos.

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

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

El segundo paso consiste en "nodificar" las líneas (producir un nodo en cada intersección). Algunas bibliotecas como STC tienen clases "Noder" que puedes usar para hacer esto, pero en PostGIS el ST_Union lo hace por ti:

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

El tercer paso consiste en crear todos los posibles polígonos no superpuestos que puedan surgir de todas esas líneas, lo que se hace con el ST_Polygonize función:

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

Puedes guardar la salida de cada uno de esos pasos en una tabla temporal, o puedes combinarlos todos en una sola sentencia:

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 es (normalmente) más conveniente tener una tabla donde cada fila es uno de los polígonos que componen la superposición de polígonos.

Answer 2

Si he entendido bien, usted quiere tomar (A es la geometría de la izquierda, B es la derecha):

Imagen de AB http://img838.imageshack.us/img838/3996/intersectab1.png

Y extracto:

AAB

Imagen de AAB http://img830.imageshack.us/img830/273/intersectab2.png

AB

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

y BAB

Imagen de BAB http://img839.imageshack.us/img839/5458/intersectab4.png

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

En primer lugar, vamos a crear una vista de todas las geometrías de intersección. Suponiendo 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 sólo lectura) que almacena pares de geomas que se intersectan, donde cada par aparece sólo una vez debido al t1.id<t2.id condición.

Ahora vamos a reunir sus intersecciones - AAB , AB y BAB utilizando el SQL UNION en 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 && se utiliza como filtro antes del operador intersects operador, para mejorar el rendimiento.
2. He optado por crear un VIEW en lugar de una consulta gigantesca; esto es sólo por comodidad.
3. Si quieres decir AB es la unión, no la intersección, de A y B - Utilice ST_Union en lugar de st_intersection en el UNION consulta en los lugares adecuados.
4. El es un signo unicode para la diferencia de Set; elimínelo del código si confunde a su base de datos.
5. Fotos por cortesía de La bonita categoría de teoría de conjuntos de Wikimedia Commons .

Answer 3

Lo que está describiendo es la forma en que el Operador de la Unión funciona en ArcGIS, pero es un poco diferente a Union o Intersection en el mundo de GEOS. El manual de Shapely tiene ejemplos de cómo funcionan los conjuntos en GEOS . Sin embargo, la wiki de PostGIS tiene un buen ejemplo utilizando el trabajo de línea que debería servirle.

Como alternativa, podrías calcular tres cosas:

1. ST_Intersección(geom_a, geom_b)
2. ST_Diferencia(geom_a, geom_b)
3. ST_Diferencia(geom_b, geom_a)

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

Answer 4

Algo así como:

INSERT INTO new_table VALUES((select id, the_geom from old_table where st_intersects(the_geom,(select the_geom from old_table where id='123')) = true

EDIT: necesitas la intersección real del polígono.

INSERT INTO new_table values((select id, ST_Intersection(the_geom,(select the_geom from old where id = 123))

a ver si eso funciona.