Creating lines to visualize differences between polygon features in PostGIS

Here are a few new tricks, using:

  • EXCEPT to remove geometries from either table that are the same, so we can focus only on geometries that are unique to each table (A_only and B_only).
  • ST_Snap to get exact noding for overlay operators.
  • Use the ST_SymDifference overlay operator to find the symmetric difference between the two geometry sets to show the differences. Update: ST_Difference shows the same result for this example. You can try either function to see what they get.

This should get what you expect:

-- CREATE OR REPLACE VIEW polygon_SymDifference AS
SELECT row_number() OVER () rn, *
FROM (
  SELECT (ST_Dump(ST_SymDifference(ST_Snap(A, B, tol), ST_Snap(B, A, tol)))).*
  FROM (
    SELECT ST_Union(DISTINCT A_only.geom) A, ST_Union(DISTINCT B_only.geom) B, 1e-5 tol
    FROM (
      SELECT ST_Boundary(geom) geom FROM polygon_a
      EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_b
    ) A_only,
    (
      SELECT ST_Boundary(geom) geom FROM polygon_b
      EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_a
    ) B_only
  ) s
) s;

 rn |                                        geom
----+-------------------------------------------------------------------------------------
  1 | LINESTRING(206.234028204842 -92.0360704110685,219.846021625456 -92.5340701703592)
  2 | LINESTRING(18.556700448873 -36.4496098325257,44.44438533894 -40.5104231486146)
  3 | LINESTRING(-131.974995802602 -38.6145334122719,-114.067738329597 -39.0215165366584)
(3 rows)

three lines


To unpack this answer a bit more, the first step with ST_Boundary gets the boundary of each polygon, rather than just the exterior. For instance, if there were holes, these would be traced by the boundary.

The EXCEPT clause is used to remove geometries from A that are part of B, and rows from B that are part of A. This reduces the number of rows that are part of A only, and part of B only. E.g., to get A_only:

SELECT ST_Boundary(geom) geom FROM polygon_a
EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_b

Here are the 6 rows of A_only, and 3 rows of B_only: A_only B_only

Next, ST_Union(DISTINCT A_only.geom) is used to combine the linework into a single geometry, typically a MultiLineString.

ST_Snap is used to snap nodes from one geometry to another. For instance ST_Snap(A, B, tol) will take the A geometry, and add more nodes from the B geometry, or move them to the B geometry, if they are within tol distance. There are probably several ways to use these functions, but the idea is to get coordinates from each geometry that are exact to each other. So the two geometries after snapping look like this:

A snapped B snapped

And to show differences, you can choose to use either ST_SymDifference or ST_Difference. They both show the same result for this example.


I think it is a little bit tricky, because of the different node sets of your both polygons (green polygon A, red different segments of polyon B). Comparing the segments of both polygons give a clue which segments of polygon B will be modified.

Nodes polygon A

poly a

Nodes of the "different" segments polygon B

seg diff

Unfortunally this show only the difference in segment structure, but I hope it is a starting point and it works like this:

After a download and unzip process I've imported the data set using PostgrSQL 9.46, PostGIS 2.1 under Debian Linux Jessie with the commands.

$ createdb gis-se
$ psql gis-se < /usr/share/postgis-2.1/postgis.sql
$ psql gis-se < /usr/share/postgis-2.1/spatial_ref_sys.sql
$ shp2pgsql -S polygon_a | psql gis-se
$ shp2pgsql -S polygon_b | psql gis-se

Assuming that the segments of polygon A are not in B and vice vera, I try to build the difference between the segments of both polygon sets, neglecting the segment membership to the polygons in each group (A or B). For didactic reasons I formulate the SQL stuff in several views.

Corresponding to this GIS-SE post, I decompose the both polygons into segment tables segments_a and segments_b

-- Segments of the polygon A
CREATE VIEW segments_a AS SELECT sp, ep
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   (SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
    FROM
    -- extract the individual linestrings
     (SELECT (ST_Dump(ST_Boundary(geom))).geom
      FROM polygon_a
     ) AS linestrings
    -- be sure that nothing is scrambled
    ORDER BY sp, ep
) AS segments;

The segment table polygon A:

SELECT 
  st_astext(sp) AS sp, 
  st_astext(ep) AS ep 
FROM segments_a 
LIMIT 3;
                    sp                     |                 ep
-------------------------------------------+--------------------------------------------
POINT(-292.268907321861 95.0342877387557)  | POINT(-287.118411917425 99.4165242769195)
POINT(-287.118411917425 99.4165242769195)  | POINT(-264.62129248575 93.2470010145007)
POINT(-277.459563916327 -44.5629543976138) | POINT(-292.268907321861 95.03428773875

The same procedure was applied to polygon B.

-- Segments of the polygon B
CREATE VIEW segments_b AS SELECT sp, ep
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   (SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
    FROM
    -- extract the individual linestrings
     (SELECT (ST_Dump(ST_Boundary(geom))).geom
      FROM polygon_b
     ) AS linestrings
    -- be sure that nothing is scrambled
    ORDER BY sp, ep
) AS segments;

The segment table polygon B

SELECT
  st_astext(sp) AS sp, 
  st_astext(ep) AS ep 
FROM segments_b 
LIMIT 3;
                    sp                     |                    ep
-------------------------------------------+-------------------------------------------
POINT(-292.268907321861 95.0342877387557)  | POINT(-287.118411917425 99.4165242769195)
POINT(-287.118411917425 99.4165242769195)  | POINT(-264.62129248575 93.2470010145007)
POINT(-277.459563916327 -44.5629543976138) | POINT(-292.268907321861 95.0342877387557)
...                        

I can build a difference table view named segments_diff_{a,b}. The difference is given by the non occurence of sorted start or end points in the segment set A and B.

CREATE VIEW segments_diff_a AS
SELECT st_makeline(b.sp, b.ep) as geom
FROM segments_b as b
LEFT JOIN segments_a as a ON (a.sp=b.sp and a.ep = b.ep)
-- filter segments without corresponding stuff in polygon A
WHERE a.sp IS NULL;

segs diff b

And the complemental stuff:

CREATE VIEW segments_diff_b AS
SELECT st_makeline(a.sp, a.ep) as geom
FROM segments_a as a
LEFT JOIN segments_b as b ON (a.sp=b.sp and a.ep = b.ep)
-- filter segments without corresponding stuff in polygon B
WHERE b.sp IS NULL;

segs diff a

Conclusion: To get a proper result for the small little segments you marked with the red arrow, both polygons must have same the node strucure and a intersection step on a node level (inserting vertexes of polygon A in B) is required. The Intersection could be done by:

CREATE VIEW segments_bi AS 
SELECT distinct sp, ep
FROM (
 SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
 FROM (
   SELECT st_difference(b.seg, a.seg) as geom FROM 
      segments_diff_a as a, segments_diff_b as b 
      WHERE st_intersects(a.seg, b.seg)
    ) as cut
  ) as segments
  WHERE sp IS NOT NULL AND ep IS NOT NULL
ORDER BY sp, ep;

But with strange results...

cut version