[GIS] Merge any and all adjacent polygons

mergepostgisunion

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.

enter image description here

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.

Best Answer

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

Related Question