Polígonos separados según la intersección usando PostGIS

32

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

  • Para un polígono dado seleccionado por id, dame todos los polígonos que se intersecan. 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 interseca con el polígono B, obtendré 3 polígonos nuevos: A menos AB, AB y B menos AB.

¿Alguna idea?

    
pregunta atogle 22.07.2010 - 23:43

4 respuestas

25

Como dijo que obtiene un grupo de polígonos que se cruzan para cada polígono en el que está interesado, es posible que desee 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:

CreoquelasolucióndeAdamcrearáunpolígono"AB" que cubre tanto el área de "AB! C" y "ABC", como un polígono "AC" que cubre "AC! B" y "ABC", y un polígono "BC" 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, de modo 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.

Hay básicamente 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 quieres manejar correctamente los agujeros]:

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

El paso 2 es "nudar" la línea (producir un nodo en cada intersección). Algunas bibliotecas como JTS tienen clases de "Noder" que puede usar para hacer esto, pero en PostGIS el ST_Union la función 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 posibles polígonos no superpuestos que pueden provenir de todas esas líneas, realizado por ST_Polygonize función:

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

    
respondido por el Jeff 05.08.2010 - 17:29
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 que se interseca.

Primero, creemos una vista de todas las geometrías que se cruzan. Suponiendo que el nombre de su tabla sea 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 geometros que se cruzan, donde cada par aparece solo una vez debido a la condición t1.id<t2.id .

Ahora reunamos sus intersecciones - A∖AB , AB y B∖AB , usando el UNION de SQL en las tres consultas:

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

UNION 

--A∖AB:
SELECT ST_Difference(t1.the_geom, t2.the_geom) 
    AS geom 
    FROM p_intersections

UNION

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

Notas:

  1. El operador && se usa como filtro antes del operador intersects , para mejorar el rendimiento.
  2. He elegido crear un VIEW en lugar de una consulta gigantesca; Esto es solo por conveniencia.
  3. Si quiso decir que AB es la unión, no la intersección, de A y B - Use ST_Union en lugar de st_intersection en la consulta UNION en los lugares apropiados.
  4. El signo es un signo de Unicode para establecer la diferencia; elimínelo del código si confunde su base de datos.
  5. Imágenes cortesía de categoría de teoría de conjuntos de Wikimedia Commons .
respondido por el Adam Matan 01.08.2010 - 13:51
8

Lo que estás describiendo es la forma en que El operador de la Unión funciona en ArcGIS, pero es un poco diferente a la Unión o la Intersección 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 que utiliza el trabajo de línea que debería hacer el truco por usted.

Alternativamente, podrías 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 deben ser los tres polígonos que mencionó en su segundo punto.

    
respondido por el scw 23.07.2010 - 04:36
-2

Algo como:

INSERTA EN VALORES new_table ((selecciona id, the_geom from old_table donde st_intersects (the_geom, (selecciona the_geom desde old_table donde id = '123')) = true

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

INTRODUZCA EN los valores new_table ((seleccione id, ST_Intersection (the_geom, (seleccione the_geom desde old donde id = 123))

ver si eso funciona.

    
respondido por el George Silva 23.07.2010 - 04:41

Lea otras preguntas en las etiquetas