Performing polygon drilldown/overlay in PostGIS?

Question

I face a challenge with PostGIS that I cannot seem to wrap my head around. I know I can solve this using a programming language (and that is my backup-plan), but I really do like to solve this in PostGIS. I've tried searching, but could not find any answers that matches my problem, this might be because I'm uncertain on my search terms, so please excuse that and point me in the right direction of there indeed is an answer.

My problem is this:

• I have a table with mixed Polygons/MultiPolygons
• Each (multi)polygon has an attribute that ranks it (priority)
• Each polygon also has a value I would like to know
• I have a search area (polygon)
• For my query area I want to find the area covered by each polygon value

Example:

Say I have the three polygons depicted in red, green, and indigo here:

And that the smaller, blue rectangle is my query polygon

In addition, the attributes are

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

What I want is to select these geometries, so that the highest ranked (green) fills all the area it can (i.e the intersection between my query geom and that geom), then the next highest (indigo) fills the intersection between the query geom and the geom MINUS the already covered) etc.

Something like this:

I've found this question: Using ST_Difference to remove overlapping features? but it does not seem to do what I want.

I can myself figure out how to compute areas and such, so a query that gives me the three geometries as depicted in the second image is fine!

Additional info: - This is not a large table (~2000 rows) - there may be zero or multiple overlaps (not just three) - there may not be any polygons in my query area (or just in parts of it) - i'm running postgis 2.3 on postgres 9.6.6

My fallback solution is to do a query like this:

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

And then iteratively "chop off" parts of the geometries in code. But, as I said, I would really like to do this in PostGIS

Answer 1

I think this works.

It is a windowing function, getting the difference between the intersection of each geometries intersection with the query box and the union of the preceding geometries.

The coalesce is needed since the union of the preceding geometries for the first geometry is null which gives a null result,instead of what is desired.

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

I am not sure how it performs though. But since both ST_Union and ST_Intersection are marked immutable it might not be that bad.

Answer 2

A bit of a different approach to this. There is a caveat that I don't know how it will scale performance wise, but on an indexed table it should be ok. It performs about the same as Nicklas's query (a tad slower?), but the measurement on such a small sample is fraught.

It looks a lot uglier than Nicklas's query, but it does avoid recursion in the query.

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

Since I babbled about WITH RECURSIVE I´ll add a quick and dirty answer using it.

This performs about as good as @NicklasAvén's solution on three Polygons, couldn´t test when upscaled.
As both solutions stand, this one has one little benefit over the other: if, for example, the Polygon with rank = 2 is contained by that of rank = 1, the ...WHERE GeometryType = 'POLYGON' filters that out while otherwise there will be a GEOMETRYCOLLECTION EMPTY (I changed the geometry of the respective Polygon in my solution accordingly to give an example; this is also true for other cases when no intersection with the difference is found). This is easily included in the other solutions, though, and might not even be of concern.

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'