Realizando drilldown / overlay de polígonos no PostGIS?

Enfrento um desafio com o PostGIS que não consigo entender. Eu sei que posso resolver isso usando uma linguagem de programação (e esse é o meu plano de backup), mas eu realmente gosto de resolver isso no PostGIS. Eu tentei pesquisar, mas não consegui encontrar respostas que correspondam ao meu problema, isso pode ser porque eu não tenho certeza sobre meus termos de pesquisa. Por isso, desculpe-me e aponte-me na direção certa; de fato, há uma resposta.

Meu problema é este:

• Eu tenho uma tabela com polígonos / poligonos múltiplos
• Cada polígono (multi) possui um atributo que o classifica (prioridade)
• Cada polígono também tem um valor que eu gostaria de saber
• Eu tenho uma área de pesquisa (polígono)
• Para a minha área de consulta, quero encontrar a área coberta por cada valor de polígono

Exemplo:

Digamos que eu tenha os três polígonos representados em vermelho, verde e índigo aqui:

E que o menor retângulo azul é o meu polígono de consulta

Além disso, os atributos são

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

O que eu quero é selecionar essas geometrias, para que a classificação mais alta (verde) preencha toda a área possível (por exemplo, a interseção entre o meu geom de consulta e esse geom), e a próxima mais alta (índigo) preencha a interseção entre o geom da consulta e o geom menos o já coberto) etc.

Algo assim:

Encontrei esta pergunta: Usando ST_Difference para remover recursos sobrepostos? mas parece não fazer o que eu quero.

Eu mesmo posso descobrir como calcular áreas e coisas assim, então uma consulta que me fornece as três geometrias mostradas na segunda imagem é boa!

Informações adicionais: - Esta não é uma tabela grande (~ 2000 linhas) - pode haver zero ou várias sobreposições (não apenas três) - pode não haver polígonos na minha área de consulta (ou apenas em partes dela) - i ' m executando o postgis 2.3 no postgres 9.6.6

Minha solução de fallback é fazer uma consulta como esta:

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

E, em seguida, iterativamente "corta" partes das geometrias no código. Mas, como eu disse, eu realmente gostaria de fazer isso no PostGIS

Eu acho que isso funciona.

É uma função de janela, obtendo a diferença entre a interseção de cada geometria com a caixa de consulta e a união das geometrias anteriores.

A coalescência é necessária, pois a união das geometrias anteriores para a primeira geometria é nula, o que fornece um resultado nulo, em vez do desejado.

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

Não tenho certeza de como ele funciona. Mas como ST_Union e ST_Intersection são marcados como imutáveis, pode não ser tão ruim assim.

Um pouco de uma abordagem diferente para isso. Há uma ressalva de que eu não sei como ele será dimensionado em termos de desempenho, mas em uma tabela indexada deve estar ok. Ele executa quase o mesmo que a consulta de Nicklas (um pouco mais lenta?), Mas a medição em uma amostra tão pequena é complicada.

Parece muito mais feio que a consulta de Nicklas, mas evita recursão na 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';

Desde que eu falei sobre isso, WITH RECURSIVEadicionarei uma resposta rápida e suja ao usá-lo.

Isso tem um desempenho tão bom quanto a solução da @ NicklasAvén em três polígonos, que não podia testar quando aumentado.
Como ambas as soluções estão, esta tem um pequeno benefício sobre a outra: se, por exemplo, o polígono com rank = 2 estiver contido no rank = 1 , os ...WHERE GeometryType = 'POLYGON'filtros que sairão enquanto houver um GEOMETRYCOLLECTION EMPTY(eu mudei a geometria do respectivo polígono em minha solução para dar um exemplo; isso também é válido para outros casos em que nenhuma interseção com a diferença é encontrada). Isso é facilmente incluído nas outras soluções e pode até não ser motivo de preocupação.

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'