How to create lines to visualize differences between polygon features in PostGIS?

Question

I've a PostGIS table polygon_b with some polygon features. There is also a table polygon_a which contains the same polygons as polygon_b but with minor changes. Now I want to create lines to visualize the differences between the polygon features.

I suppose that ST_ExteriorRing and ST_Difference will do the job but the WHERE clause seems to be quite tricky.

CREATE VIEW line_difference AS SELECT
row_number() over() AS gid,
g.geom::geometry(LineString, yourSRID) AS geom
FROM 
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(ST_ExteriorRing(polygon_a.geom), ST_ExteriorRing(polygon_b.geom))))).geom AS geom
    FROM polygon_a, polygon_b
    WHERE 
    -- ?
    ) AS g;

Can anyone help me?

EDIT 1

As posted by 'tilt' I've tried ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom) but the result is not as expected.

CREATE VIEW line_difference AS SELECT
row_number() over() AS gid,
g.geom::geometry(LineString, your_SRID) AS geom
FROM 
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(ST_ExteriorRing(polygon_a.geom), ST_ExteriorRing(polygon_b.geom))))).geom AS geom
    FROM polygon_a, polygon_b
    WHERE 
    ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom))
     AS g;

EDIT 2

workupload.com/file/J0WBvRBb (example dataset)

---

I've tried to turn the polygons into multilines before using ST_Difference, but the results are still strange.

CREATE VIEW multiline_a AS SELECT
row_number() over() as gid,
ST_Union(ST_ExteriorRIng(polygon_a.geom))::geometry(multilinestring, 4326) AS geom
FROM
polygon_a;

CREATE VIEW multiline_b AS SELECT
row_number() over() as gid,
ST_Union(ST_ExteriorRIng(polygon_b.geom))::geometry(multilinestring, 4326) AS geom
FROM
polygon_b;

CREATE VIEW line_difference AS SELECT
row_number() over() as gid,
g.geom
FROM
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(multiline_a.geom, multiline_b.geom)))).geom::geometry(linestring, 4326) AS geom
    FROM
    multiline_a, multiline_b)
As g;

Answer 1

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)

---

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:

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:

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

Answer 2

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

Nodes of the "different" segments polygon B

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;

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;

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...

Answer 3

Looking at the example, the change implies that features from the new table that have been changed will always be overlapping features from the old table. Therefore you would be done with

ST_Overlaps(geoma, geomb) AND !ST_Touches(geoma,geomb)

The negation on touches is because features also overlap if only their borders share the same vertice locations.