Monday 27 January 2020

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


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.


enter image description here enter image description here enter image description here


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;

enter image description here


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;

enter image description here



Answer



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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...