What I basically want is:

enter image description here

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:

enter image description here

And I need the red part.

enter image description here

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

    geom, (
        SELECT ST_Union(buff)
        FROM buffers
            buffers.id = ANY(z.ids) AND
            z.id != buffers.id
    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


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

  • Could you create intersections from all permutations of circles and keep the one that overlaps the most with the circles in the first place? Nov 1 '21 at 8:40
  • I mean, by just looking the image i need, my eyes basically can see there is darkest blue region but when it comes to make it into SQL, I just can't write it. It's the same thing when you say that, thanks for the help but not sure if I can make it in SQL. Nov 1 '21 at 10:47

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:

    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,
        FROM   clst
        WHERE  _rnk = 1
        UNION ALL
        SELECT its.ids || a.id,
               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
FROM   its
WHERE  _rnk = __mrnk


  • 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!

  • I was thinking about iterative intersection too but couldn't think of a way to do that in postgres so I went with python this time, I basically went with two for loops, got the intersecting areas, extracted the sum intersection. I think it would be great in postgis if there was some implementation of st_intersection aggregation. Could have been much easier than your code I think. Nov 1 '21 at 13:14
  • Thanks for this great example of work, going through detailing it helped me better understand that but looks like my current python code is in a more flexible state, making it a better choice for me this time. Nov 1 '21 at 13:16
  • Objectively easier is not necessarily better - this query is accurate and a lot faster than other approaches I have seen, including attempts to pack this up in functions etc. Despite that, ST_Intersection (or any other function) does not handle this because of the multitude of cases that will arise (and break certain approaches) when handling cascading overlap operations.
    – geozelot
    Nov 1 '21 at 13:31

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...



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.