La realización de polígono drilldown/superposición de capas PostGIS?

Question

Me enfrento a un desafío con PostGIS que yo no puedo envolver mi cabeza alrededor. Sé que puedo resolver esto mediante un lenguaje de programación (y que es mi copia de seguridad-plan), pero realmente no se como resolver esto en PostGIS. He intentado buscar, pero no pudo encontrar ninguna respuesta que coincide con mi problema, esto podría ser porque estoy seguro de mis términos de búsqueda, así que por favor perdonen y me apunte en la dirección correcta de allí, de hecho, es una respuesta.

Mi problema es el siguiente:

• Tengo una tabla con una mezcla de Polígonos/Multipolígonos
• Cada uno (multi), el polígono tiene un atributo que se ubica (prioridad)
• Cada polígono tiene también un valor que me gustaría saber
• Tengo un área de búsqueda (polígono)
• Para mi área de consulta, quiero encontrar el área cubierta por cada polígono valor

Ejemplo:

Decir que tengo los tres polígonos que se representan en rojo, verde, añil y aquí:

Y que la más pequeña, rectángulo azul es mi consulta polígono

Además, los atributos son

geom   | rank | value
-------|------|----  
red    |  3   | 0.1
green  |  1   | 0.2
indigo |  2   | 0.2

Lo que quiero es seleccionar estas geometrías, de modo que el valor más alto (verde) ocupa todo el área de lo que puede (que yo.e la intersección entre mi consulta geom y que geom), a continuación, el siguiente más alto (índigo) llena la intersección entre la consulta geom y la geom MENOS los ya cubiertos) etc.

Algo como esto:

He encontrado esta pregunta: Utilizando ST_Difference para eliminar la superposición de funciones? pero no parece hacer lo que quiero.

Puedo yo entender cómo calcular áreas y tal, así que una consulta que me da las tres geometrías como se muestra en la segunda imagen se ve bien!

Información adicional: - Esta no es una tabla grande (~2000 filas) - no puede ser cero o varias superposiciones (no tres) - no puede ser de cualquiera de los polígonos en mi área de consulta (o sólo en algunas de sus partes) - me estoy quedando postgis 2.3 en postgres 9.6.6

Mi solución alternativa es hacer una consulta como esta:

SELECT 
ST_Intersection(geom, querygeom) as intersection, rank, value
FROM mytable
WHERE ST_Intersects(geom, querygeom)
ORDER by rank asc

Y, a continuación, de forma iterativa "cortar" las partes de las geometrías en el código. Pero, como he dicho, realmente me gustaría hacer esto en PostGIS

Answer 1

Creo que esto funciona.

Es una función de ventanas, obteniendo la diferencia entre la intersección de cada una de las geometrías de la intersección con el cuadro de consulta y la unión de la anterior geometrías.

El unen es necesario ya que la unión de la anterior geometrías para el primer geometría es nula, lo que da un resultado nulo,en lugar de lo que se desea.

WITH a(geom, rank, val) AS
(
    VALUES
    ('POLYGON((1 1, 1 5, 5 5, 5 1, 1 1))'::geometry,3,0.1),
    ('POLYGON((2 3, 2 8, 4 8, 5 3,2 3))'::geometry,1,0.2),
    ('POLYGON((4 4, 4 6, 6 6, 6 4,4 4))'::geometry,2,0.2)
)
,q AS
(
    SELECT 'POLYGON((3 3, 3 4.5, 12 4.5,12 3, 3 3))'::geometry geom
) 
SELECT 
  ST_Difference(
    ST_Intersection(a.geom, q.geom), 
    COALESCE(ST_Union(a.geom) 
           OVER (ORDER BY rank ROWS BETWEEN UNBOUNDED PRECEDING and 1 PRECEDING),
       'POLYGON EMPTY'::geometry)
  ) geom 
FROM a,q
WHERE ST_Intersects(a.geom, q.geom);

No estoy seguro de cómo funciona, aunque. Pero desde que ambos ST_Union y ST_Intersection están marcados inmutable no puede ser que lo malo.

Answer 2

Un poco de un enfoque diferente a este. Hay un problema que yo no sé cómo va a escalar el rendimiento sabio, pero en una tabla indizada se debe aceptar. Se realiza sobre el mismo como Nicklas de la consulta (un poco más lento?), pero la medición en una muestra tan pequeña está cargada.

Su aspecto es mucho más feo que Nicklas de la consulta, pero sí evitar la recursividad en la consulta.

WITH 
    a(geom, rank, val) AS
    (
        VALUES
        ('POLYGON((1 1, 1 5, 5 5, 5 1, 1 1))'::geometry,3,0.1),
        ('POLYGON((2 3, 2 8, 4 8, 5 3, 2 3))'::geometry,1,0.2),
        ('POLYGON((4 4, 4 6, 6 6, 6 4, 4 4))'::geometry,2,0.2)
    ),
    polygonized AS (
        -- get the basic building block polygons
        SELECT (ST_Dump(         -- 5. Dump out the polygons
            ST_Polygonize(line)) -- 4. Polygonise the linework
            ).geom AS mypoly
        FROM (
            SELECT 
                ST_Node(                  -- 3. Node lines on intersection
                    ST_Union(             -- 2. Union them for noding
                        ST_Boundary(geom) -- 1. Change to linestrings
                    ) 
                ) 
                AS line
            FROM a
        ) line
    ),
    overlayed AS ( 
        -- Overlay with original polygons and get minimum rank.
        -- Union and dissolve interior boundaries for like ranks.
        SELECT (ST_Dump(ST_UnaryUnion(geom))).geom, rank 
        FROM ( 
            -- union up the polygons by rank, unary union doesn't count as an aggregate function?
            SELECT ST_Union(mypoly) geom, rank 
            FROM ( 
                -- get the minimum rank for each of the polygons
                SELECT p.mypoly, min(rank) rank
                FROM polygonized p 
                    INNER JOIN a ON ST_Contains(a.geom,ST_PointOnSurface(p.mypoly))
                GROUP BY p.mypoly
                ) g
            GROUP BY rank
            ) r
    )
-- get the intersection of the query area with the overlayed polygons
SELECT ST_Intersection(o.geom,'POLYGON((3 3, 3 4.5, 12 4.5,12 3, 3 3))'::Geometry), rank
FROM overlayed o
WHERE ST_Intersects(o.geom,'POLYGON((3 3, 3 4.5, 12 4.5,12 3, 3 3))'::Geometry) and
    -- Intersection can do funky things
    GeometryType(ST_Intersection(o.geom,'POLYGON((3 3, 3 4.5, 12 4.5,12 3, 3 3))'::Geometry)) like '%POLYGON';

Answer 3

Desde que me balbuceó sobre WITH RECURSIVE Enfermo añadir un rápido y sucio respuesta usando.

Esta se realiza sobre tan buena como la de @NicklasAvén la solución en tres Polígonos, no podía prueba cuando inflados.
Como ambas soluciones stand, esta tiene un poco de ventaja sobre los otros: si, por ejemplo, el Polígono con rango = 2 es contenido por el de rank = 1, ...WHERE GeometryType = 'POLYGON' filtros que fuera, mientras que de lo contrario habrá un GEOMETRYCOLLECTION EMPTY (he cambiado la geometría de los respectivos Polígonos en mi solución, en consecuencia, para dar un ejemplo; esto también es cierto para otros casos, cuando no hay intersección con la diferencia se encuentra). Esto es fácilmente incluidos en el resto de soluciones, sin embargo, y podría incluso no ser motivo de preocupación.

WITH RECURSIVE
    a(geom, rank, val) AS (
        VALUES
           ('POLYGON((1 1, 1 5, 5 5, 5 1, 1 1))'::geometry,3,0.1),
           ('POLYGON((2 3, 2 8, 4 8, 5 3,2 3))'::geometry,1,0.2),
           ('POLYGON((2.1 3.1, 2.1 7.9, 3.9 7.9, 4.9 3.1,2.1 3.1))'::geometry,2,0.2)
    ),
    q AS (
        SELECT 'POLYGON((3 3, 3 4.5, 12 4.5,12 3, 3 3))'::geometry geom
    ),
    src AS (
        SELECT ROW_NUMBER() OVER(ORDER BY rank) AS rn,
               ST_Intersection(q.geom, a.geom) AS geom,
               rank,
               val
        FROM a
        JOIN q
           ON ST_Intersects(a.geom, q.geom)
    ),
    res AS (
        SELECT s.geom AS its,
               ST_Difference(q.geom, s.geom) AS dif,
               s.rank,
               s.val,
               2 AS icr
        FROM src AS s,
             q
        WHERE s.rn = 1
        UNION ALL
        SELECT ST_Intersection(s.geom, r.dif) AS its,
               ST_Difference(r.dif, s.geom) AS dif,
               s.rank,
               s.val,
               icr + 1 AS icr
        FROM src AS s,
             res AS r
        WHERE s.rank = r.icr
    )

SELECT its AS geom,
       rank,
       val
FROM res
WHERE GeometryType(its) = 'POLYGON'