Using PostGIS to generate a polygon representing the area shared by some (but not all) of a set of polygons

Question

I have an application where several users submit polygons, and I am looking to create an aggregate polygon that represents a kind of average of their submissions. We use PostGIS for our spatial analysis needs.

Let's say we received polygons from 10 users. If we used ST_Intersection on those polygons, the remaining polygon would only represent the points included in all 10 polygons. If we used ST_Union, the output would repesent the points included in at least 1 polygon.

Can anyone recommend a way to output a polygon that represents the points that are in n polygons, where n is greater than 1 and less than the total number of polygons (10 in this case)?

Answer 1

Just spit-balling here:

• ST_GeneratePoints() on each input geometry with N points.
• Randomize that point set and take a 1/M of them.
• Build voronoi polygons of that set.
• Spatial join the voronoi polygons to the original point set and only retain those polygons with more than P points in them
• Union those polygons.
• Output the result.

What do you think?

Answer 2

Some options, from the top of my head; not tested and not optimized.

• To get the 'average' point set area for points with 1 < cnt < total overlaps:

  SELECT  ST_ConvexHull(ST_Collect(geom)) AS geom
  FROM    (
      SELECT  dmp.geom, COUNT(DISTINCT a.*) AS cnt
      FROM    mask AS a,
              mask AS b,
              LATERAL ST_DumpPoints(b.geom) AS dmp
      WHERE   ST_Intersects(a.geom, dmp.geom)
      GROUP BY
              1
  ) q
  WHERE   cnt > 1 AND cnt < (SELECT COUNT(*) FROM mask)
  ;
  

• To get the precise polygon area from parts with 1 < cnt < total overlaps:

  SELECT  ST_Union(geom) AS geom
  FROM    (
      SELECT  a.geom, COUNT(b.*) AS cnt
      FROM    (
          SELECT  (ST_Dump(ST_Split(ST_Union(a.geom), ST_Union(ST_ExteriorRing(b.geom))))).geom
          FROM    mask AS a
          JOIN    mask AS b
            ON    a.id <> b.id AND ST_Intersects(a.geom, b.geom)
      ) AS a
      JOIN    mask AS b
        ON    ST_Intersects(ST_Centroid(a.geom), b.geom)
      GROUP BY
              a.geom
  ) q
  WHERE cnt > 1 AND cnt < (SELECT COUNT(*) FROM mask)
  ;
  

---

Note that the first query guarantees a single-part polygon, but with a generalized shape over the point set, while the second will stitch together exact polygon parts, but may result in a multi-part polygon.

Answer 3

One of the people who replied (see Paul Ramsey above) posted a great answer on this blog here: https://info.crunchydata.com/blog/polygon-averaging-in-postgis.