Get resulting polygon that is overlapped by others

Question

What I basically want is:

I want to have the polygon which is numbered as 13. As you can see there are many overlapping polygons but the most overlapping one is numbered 13. The image is from this answer.

Also there is one another question on stackexchange that is related to my question.

In my case I have these polygons:

And I need the red part.

What I have tried so far is this:

WITH pois AS (
    SELECT unnest(ARRAY['1','2','3','4','5','6','7']) AS id, unnest(ARRAY[41.720302,41.732075701,41.727253557,41.7503713,41.723427674,41.786769,41.77869826276586]) AS lat, unnest(ARRAY[38.488551,38.525743929,38.5174403600001,38.397713078,38.4954947360001,38.605068,38.4268728196621]) AS lon
    ORDER BY id ASC
), buffers AS (
    SELECT id, ST_Transform(ST_Buffer(ST_Transform(ST_SetSRID(ST_Point(lon, lat), 4326), 5275), 3000), 4326) buff, ST_SetSRID(ST_Point(lon, lat), 4326) pt
    FROM pois
), nodes AS (
    SELECT MAX(id) as id, (ST_Dump(ST_Node(ST_Collect(ST_ExteriorRing(buff))))).geom AS geom
    FROM buffers
)
SELECT ,
ST_Difference(
    geom, (
        SELECT ST_Union(buff)
        FROM buffers
        WHERE
            buffers.id = ANY(z.ids) AND
            z.id != buffers.id
    )
)
FROM (
    SELECT count() AS count, array_agg(c.id) ids, p.id, ST_GeomFromText(MAX(ST_AsEWKT(p.buff))) geom, p.pt
    FROM buffers p 
    JOIN buffers c 
    ON ST_Contains(c.buff, ST_PointOnSurface(p.buff)) 
    GROUP BY (p.id, p.pt)
    ORDER BY count ASC
) AS z
WHERE count > 1

EDIT

So turns out, in python I can iteratively intersect areas, it maybe slower than some database query but looks like it gets work done. It is hard to make this example into a working code without my database and modules but at least it may give some ideas to other people out there and it also is not optimized by any means at this state:

targets = {
    b_id: {
        "point": GEOSPoint(i["lat"], i["lon"], srid=4326),
        "buffer": BufferAnalysis.create_buffer(i["lat"], i["lon"], distance / 1000)
    } for b_id, i in df.iterrows()
}
sections = set()
for b_id, polygon in targets.items():
    intersections = []
    for sub_id, sub_polygon in targets.items():
        if polygon["buffer"] & sub_polygon["buffer"] and id != sub_id:
            intersections.append([sub_id, sub_polygon])
    if len(intersections) > 1:
        for i in range(len(intersections), 1, -1):
            tmp_poly = polygon["buffer"]
            areas = intersections[:i]
            info = "|".join(set([i[0] for i in areas] + [b_id]))
            if info not in sections:
                print(b_id, info)
                sections.add(info)
                for sub_id, sub_polygon in areas:
                    tmp_poly = tmp_poly.intersection(sub_polygon["buffer"])

I added sections and range selection with length of intersections are just to get other intersections like, 4 regions were to intersect but there are other combinations 3 regions in these 4 regions, it is not correct in this example but at least it gives me the most intersecting area. If I come up with a better idea, I will add it here.

Answer 1

There's a reason this very popular problem and task has no generic answer: it needs to be tailored to your data to produce predictable results - and has to be implemented with different algorithmic approaches. In terms of algorithmic theory you are required to set up specific rules to address your specific case.

On top of that you will often need to create temporary auxiliary columns or on-the-fly inline storage (e.g. arrays) that cannot benefit from indexes, apart from adding these to the base tables in advance - which, together with the iterative ST_Intersection calls, makes this a very expensive procedure.

Many of these stacked overlay applications need to operate on tables as well as sets of rows, both in intermediate steps and to generate the results - packing these into procedural language functions tends to be awkward, despite the arguably better (but definitely more intuitive) handling of iteration.

That being said, DO blocks (or the new PROCEDURES) are great if the output is a new table, and functions are effective when the overlay is applied on a per-row basis.

Here, however, we can utilize the built-in SQL level iteration concept of (the misleadingly named) RECURSIVE CTEs. This query will suit your example well, and is sufficiently performant. However, as mentioned above, it (or rather, the underlying concept) is limited to the assumption that all geometries of a cluster (of intersecting geometries) participate in the final common overlap area.

Unfortunately, RECURSIVE CTEs are lengthy and perplexing if not familiar; I at least assume here that the inline data in your example (buffer -> geom) was created into a table called aois only to lessen the lines of code:

WITH RECURSIVE
    clst AS (
        SELECT *,
               MAX(_rnk) OVER(PARTITION BY _cid) AS __mrnk
        FROM   (
            SELECT *,
                   ROW_NUMBER() OVER(PARTITION BY _cid) AS _rnk
            FROM   (
                SELECT *,
                       ST_ClusterDBSCAN(geom, 0, 1) OVER() AS _cid
                FROM   aois
            ) q
        ) q
    ),
    its AS (
        SELECT ARRAY[id] AS ids,
               _rnk,
               __mrnk,
               _cid,
               geom
        FROM   clst
        WHERE  _rnk = 1
        UNION ALL
        SELECT its.ids || a.id,
               a.ida._rnk,
               a.__mrnk,
               a._cid,
               ST_Intersection(a.geom, its.geom) AS geom
        FROM   clst AS a
        JOIN   its
          ON   a._cid = its._cid
         AND   a._rnk = its._rnk+1
    )
SELECT ids,
       geom
FROM   its
WHERE  _rnk = __mrnk
;

Here

• the clst block adds a set of auxiliary columns to the base data:
  • _cid is an id of clusters of intersecting geometries to operate the (spatial) overlap iteration on
  • _rnk is a sequential id of rows within each cluster to enable iteration per cluster without the need for spatial filters
  • __mrnk is the maximum _rnk per cluster, solely to help with the final filtering of the result set (a RECURSIVE CTE will add rows from every iteration to the result set, which we are not interested in)
• the its block is the actual iteration process, which produces an array of participating aois.id for your reference

No attempts to homogenize intermediate GEOMETRYCOLLECTIONs or different geometry types in the output were made!

Answer 2

So, there are many ways to solve this problem and this is one of them based on the vector geometry method :-).

1. Create a "fun" and "dangerous" custom function ST_CarvesPolygons();

2. Run this construct, where I tried to adjust it to your test data and check the result:

WITH pois AS (SELECT unnest(ARRAY['1','2','3','4','5','6','7']) AS id, unnest(ARRAY[41.720302,41.732075701,41.727253557,41.7503713,41.723427674,41.786769,41.77869826276586]) AS lat,
unnest(ARRAY[38.488551,38.525743929,38.5174403600001,38.397713078,38.4954947360001,38.605068,38.4268728196621]) AS lon ORDER BY id ASC),
    buffers AS (SELECT id, ST_Transform(ST_Buffer(ST_Transform(ST_SetSRID(ST_Point(lon, lat), 4326), 5275), 3000), 4326) bufgeom, 
          ST_SetSRID(ST_Point(lon, lat), 4326) geom FROM pois),
    carves AS (SELECT DISTINCT ST_CarvesPolygons(a.bufgeom, b.bufgeom) geom FROM buffers a JOIN buffers b ON true),
    dis_buffers AS (SELECT geom FROM ST_Dump((SELECT ST_Polygonize(geom) AS geom FROM (SELECT ST_Union(geom) AS geom FROM (SELECT ST_ExteriorRing(geom) AS geom FROM carves) AS lines) AS foo)))
    SELECT DISTINCT ST_Difference(a.geom, ST_Union(b.geom)) geom FROM dis_buffers a JOIN carves b ON true GROUP BY a.geom;

Applied to your drawings and data, the proposed design should behave as expected...

Original spatial solutions...

https://www.deepl.com/ru/translator