Merge any and all adjacent polygons

Question

I would like to do adjacency tests on a parcel (polygons) layer and merge them if they fit certain criteria (could be size). Per the picture below, I would like to merge polygons 1,2,3 and 4, but not 5.

I have two problems:

1. ST_TOUCHES returns TRUE if just the corners touch and not a line segment. I think I need ST_RELATE to check for shared line segments.
2. Ideally, I would like to merge ALL adjacent polygons into one, but I am not sure how to scale beyond two--as in, merge 1,2,3 and 4 (and possibly more on actual data) in one round.

The structure I have now is based on a self join on ST_TOUCHES.

Toy data

CREATE TABLE testpoly AS 
SELECT 
1 AS id, ST_PolyFromText('POLYGON ((0 0, 10 0, 10 20, 00 20, 0 0 ))') AS geom UNION SELECT
2 AS id, ST_PolyFromText('POLYGON ((10 0, 20 0, 20 20, 10 20, 10 0 ))') AS geom UNION SELECT
3 AS id, ST_PolyFromText('POLYGON ((10 -20, 20 -20, 20 0, 10 0, 10 -20 ))') AS geom UNION SELECT
4 AS id, ST_PolyFromText('POLYGON ((20 -20, 30 -20, 30 0, 20 0, 20 -20 ))') AS geom  UNION SELECT 
5 AS id, ST_PolyFromText('POLYGON ((30 0, 40 0, 40 20, 30 20, 30 0 ))') AS geom ;

Selection

SELECT 
    gid, adj_gid,
    st_AStext(st_union(l2.g1,l2.g2)) AS geo_combo
from (
    --level 2
    SELECT
      t1.id AS gid,
      t1.geom AS g1,
      t2.id AS adj_gid,
      t2.geom AS g2
     from
      testpoly  t1,
      testpoly  t2
     where
      ST_Touches( t1.geom, t2.geom ) 
      AND t1.geom && t2.geom 
) 
l2

Here is the output:

+-----+---------+-------------------------------------------------------------------------------+
| gid | adj_gid | geo_combo                                                                     |
+-----+---------+-------------------------------------------------------------------------------+
| 1   | 2       | POLYGON((10 0,0 0,0 20,10 20,20 20,20 0,10 0))                                |
+-----+---------+-------------------------------------------------------------------------------+
| 1   | 3       | MULTIPOLYGON(((10 0,0 0,0 20,10 20,10 0)),((10 0,20 0,20 -20,10 -20,10 0)))   |
+-----+---------+-------------------------------------------------------------------------------+
| 2   | 1       | POLYGON((10 20,20 20,20 0,10 0,0 0,0 20,10 20))                               |
+-----+---------+-------------------------------------------------------------------------------+
| 2   | 3       | POLYGON((10 0,10 20,20 20,20 0,20 -20,10 -20,10 0))                           |
+-----+---------+-------------------------------------------------------------------------------+
| 2   | 4       | MULTIPOLYGON(((20 0,10 0,10 20,20 20,20 0)),((20 0,30 0,30 -20,20 -20,20 0))) |
+-----+---------+-------------------------------------------------------------------------------+
| 3   | 1       | MULTIPOLYGON(((10 0,20 0,20 -20,10 -20,10 0)),((10 0,0 0,0 20,10 20,10 0)))   |
+-----+---------+-------------------------------------------------------------------------------+
| 3   | 2       | POLYGON((20 0,20 -20,10 -20,10 0,10 20,20 20,20 0))                           |
+-----+---------+-------------------------------------------------------------------------------+
| 3   | 4       | POLYGON((20 -20,10 -20,10 0,20 0,30 0,30 -20,20 -20))                         |
+-----+---------+-------------------------------------------------------------------------------+
| 4   | 2       | MULTIPOLYGON(((20 0,30 0,30 -20,20 -20,20 0)),((20 0,10 0,10 20,20 20,20 0))) |
+-----+---------+-------------------------------------------------------------------------------+
| 4   | 3       | POLYGON((20 0,30 0,30 -20,20 -20,10 -20,10 0,20 0))                           |
+-----+---------+-------------------------------------------------------------------------------+
| 4   | 5       | MULTIPOLYGON(((30 0,30 -20,20 -20,20 0,30 0)),((30 0,30 20,40 20,40 0,30 0))) |
+-----+---------+-------------------------------------------------------------------------------+
| 5   | 4       | MULTIPOLYGON(((30 0,30 20,40 20,40 0,30 0)),((30 0,30 -20,20 -20,20 0,30 0))) |
+-----+---------+-------------------------------------------------------------------------------+

Note that polygon id=3 shares a point with id=1 and thus is returned as a positive result. If I change the WHERE clause to ST_Touches( t1.geom, t2.geom ) AND t1.geom && t2.geom AND ST_Relate(t1.geom, t2.geom ,'T*T***T**'); I get no records at all.

1. So first, how do I specify ST_Relate to make sure only parcels sharing a line segment are considered.

2. And then, how would I merge polygons 1,2,3,4 in one round, collapsing the results from the above call, all the while recognizing that adjacency 1 to 2 is the same as the reverse?

Update

If I add this to the where clause I obviously only get polygons and not multipolygons, thus weeding out false positives for my purposes--corner touches will be ignored.

GeometryType(st_union(t1.geom,t2.geom)) != 'MULTIPOLYGON'

While this is not ideal (I would rather use topology checks with ST_RELATE as a more general solution), it is a way forward. Then remains the matter of de-duping and union'ing these. Possibly, if I could generate a sequence for only polygons touching, I could union on that.

Update II

This one seems to work for selecting polygons sharing lines (but not corners) and is thus a more general solution than the above MULTIPOLYGON test. My where clause now looks like this:

WHERE
              ST_Touches( t1.geom, t2.geom ) 
              AND t1.geom && t2.geom 
              
              -- 'overlap' relation
              AND ST_Relate(t1.geom, t2.geom)='FF2F11212') t2 

Now what remains is still how to do the merge for more than just a pair of polygons, but for an arbitrary number fitting the criteria, in one go.

Answer 1

I couldn't help thinking that your example is actually a raster and although you mentioned that you would like to merge based on "certain criteria (could be size)" I would like to give it a shot with a raster conversion.

For your specific example this would work:

WITH rast AS (
  SELECT 
  ST_UNION(ST_AsRaster(geom,10, 20, '2BUI')) r
  FROM testpoly 
)
,p AS (
    SELECT (ST_DumpAsPolygons(r)).geom FROM rast
)
SELECT t.id,p.* 
FROM p
LEFT JOIN testpoly  t ON ST_Equals(p.geom, t.geom)

What happens is that since your polygons are perfectly aligned cells, they will convert nicely into a raster (10x20 cellsize). The dumpaspolygons helps you here by merging all adjacent cells into one and by comparing with the original polygons you will even be able to get the id back for non-merged poly's.

Having explained this, I am very curious how this would scale and how large your dataset is :D

Answer 2

Here's an example of how to do this in procedural style with multiple passes under the hood.

CREATE TABLE joined_testpoly AS SELECT array[id] ids, geom FROM testpoly; 

You should be able to carry along more columns and apply extra criteria for joining by modifying how the LIMIT 1 select below works:

CREATE OR REPLACE FUNCTION reduce_joined_testpoly()
RETURNS void
AS $$
DECLARE
  joined_row joined_testpoly%ROWTYPE;
BEGIN
  LOOP
     SELECT array_cat(a.ids, b.ids), st_union(a.geom, b.geom)
         INTO joined_row 
     FROM joined_testpoly a INNER JOIN joined_testpoly b
           on a.ids != b.ids
              and ST_Touches(a.geom, b.geom) and a.geom && b.geom 
              and ST_Relate(a.geom, b.geom)='FF2F11212'
         LIMIT 1;
     IF NOT FOUND THEN
           EXIT;
     END IF;
     INSERT INTO joined_testpoly VALUES (joined_row.ids, joined_row.geom);
     DELETE FROM joined_testpoly
         WHERE joined_testpoly.ids <@ joined_row.ids 
           AND joined_testpoly.ids != joined_row.ids;
  END LOOP;
  RETURN;
END;
$$ LANGUAGE plpgsql;

Run the thing:

SELECT reduce_joined_testpoly();

Proper unions, no multipolygons:

SELECT ids, st_geometrytype(geom), st_area(geom), st_numgeometries(geom) 
FROM joined_testpoly;
    ids    | st_geometrytype | st_area | st_numgeometries 
-----------+-----------------+---------+------------------
 {5}       | ST_Polygon      |     200 |                1
 {1,2,3,4} | ST_Polygon      |     800 |                1

Answer 3

Here's another (not working) strategy for reference (that I couldn't get to exclude the single touching point case). It should be faster than my other answer as it only takes 'one pass'.

SELECT st_numgeometries(g), (SELECT st_union(x.geom) FROM st_dump(g) x GROUP BY g)
FROM (
    SELECT unnest(st_clusterintersecting(geom)) g, id < 100 as other_arbitrary_grouping 
    FROM testpoly
    GROUP BY other_arbitrary_grouping) c;

(feel free to amend and post another answer if anyone can get the id=5 geometry in it's own group)

To get back the list of ids etc. you'd have to use st_contains to rejoin on the testpoly table as detailed in the following answer: https://stackoverflow.com/a/37486732/6691 but I couldn't get that to work for polygons for some reason.

Answer 4

Here's a quick stab at it using your original query adjusted a bit:

with gr as (SELECT 
    gid, adj_gid,
    st_AStext(st_union(l2.g1,l2.g2)) AS geo_combo
from (
    --level 2
    SELECT
      t1.id AS gid,
      t1.geom AS g1,
      t2.id AS adj_gid,
      t2.geom AS g2
     from
      testpoly  t1,
      testpoly  t2
     where
      ST_Touches( t1.geom, t2.geom ) 
      AND ST_Relate(t1.geom,t2.geom, '****1****')
      AND t1.geom && t2.geom 
) 
l2) select ST_AsText(st_union(gr.geo_combo)) from gr;

References: https://postgis.net/docs/using_postgis_dbmanagement.html#DE-9IM