3

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

2
  • 2
    Interesting requirement. Can you post a picture of a typical set of polygons?
    – dr_jts
    Apr 6, 2020 at 22:32
  • 1
    I think one approach would be to form a coverage out of the linework from the polygons, and then extract the area created by the coverage polygons which have at least N parents.
    – dr_jts
    Apr 6, 2020 at 22:33

3 Answers 3

3

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?

1
1

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.

0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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