I try to calculate the symmetric difference for a set of segment groups A and B with precision of 1E-5 (ST_SnapedToGrid(geom,1E-5). I've formulated the an ideal test case for a H based shape segment group A and the segment group B containing only the vertical bars of the H. If I build the symmetric difference, I get the expected result, the horizontal bar of the H.
gis-se=# SELECT
ST_AsText((ST_Dump(
ST_SymDifference(
ST_GeometryFromText('MULTILINESTRING(
(-1 1, -1 0),
(-1 0, 1 0),
( 1 0, 1 1),
(-1 -1, -1 0),
(-1 0, 1 0),
( 1 0, 1 -1))'),
ST_GeometryFromText('MULTILINESTRING(
(-1 1, -1 -1),
( 1 1, 1 -1))')
))).geom) AS DUMP;
dump
----------------------
LINESTRING(-1 0,1 0)
If I use a similar setup like in this real world example, I get an unexpected result.
SELECT
ST_AsText((ST_Dump(
ST_SymDifference(
-- Segment group A the "H"
ST_GeometryFromText('MULTILINESTRING(
(219.84602 -92.53407,223.95468 -131.19481),
(212.94651 -27.6126,219.84602 -92.53407),
(-114.06774 -39.02152,-107.63268 -99.71427),
(-125.67524 -98.00064,-131.975 -38.61453),
(-131.975 -38.61453,-114.06774 -39.02152),
(44.44439 -40.51042,52.34093 -114.89627),
(26.19878 -108.47147,18.5567 -36.44961),
(210.25568 -129.89414,206.23403 -92.03607),
(219.84602 -92.53407,206.23403 -92.03607),
(206.23403 -92.03607,194.92038 14.4658),
(-131.975 -38.61453,-137.6955 15.31113),
(-120.60474 22.63271,-114.06774 -39.02152),
(35.93373 39.66043,44.44439 -40.51042),
(18.5567 -36.44961,10.55734 38.93942)
)'),
-- Segment group B The only the vertical bars of the "H"
ST_GeometryFromText('MULTILINESTRING(
(212.94651 -27.6126,223.95468 -131.19481),
(210.25568 -129.89414,194.92038 14.4658),
(-125.67524 -98.00064,-137.6955 15.31113),
(-120.60474 22.63271,-107.63268 -99.71427),
(35.93373 39.66043,52.34093 -114.89627),
(26.19878 -108.47147,10.55734 38.93942))')
))).geom);
Segment group A
Segment group B
Result is the complete segment group A.
dump
----------------------------------------------------------------------
LINESTRING(219.84602 -92.53407,223.95468 -131.19481)
LINESTRING(212.94651 -27.6126,219.84602 -92.53407)
LINESTRING(-114.06774 -39.02152,-107.63268 -99.71427)
LINESTRING(-125.67524 -98.00064,-131.975 -38.61453)
LINESTRING(-131.975 -38.61453,-131.974995810398 -38.6145300952198)
LINESTRING(-131.974995810398 -38.6145300952198,-114.06774 -39.02152)
LINESTRING(44.44439 -40.51042,52.34093 -114.89627)
LINESTRING(26.19878 -108.47147,18.5567 -36.44961)
LINESTRING(210.25568 -129.89414,206.23403 -92.03607)
LINESTRING(219.84602 -92.53407,206.234032029431 -92.0360700742475)
LINESTRING(206.234032029431 -92.0360700742475,206.23403 -92.03607)
LINESTRING(206.23403 -92.03607,194.92038 14.4658)
LINESTRING(-131.975 -38.61453,-137.6955 15.31113)
LINESTRING(-120.60474 22.63271,-114.06774 -39.02152)
LINESTRING(35.93373 39.66043,44.44439 -40.51042)
LINESTRING(18.5567 -36.44961,10.55734 38.93942)
LINESTRING(212.94651 -27.6126,223.95468 -131.19481)
LINESTRING(210.25568 -129.89414,206.234032029431 -92.0360700742475)
LINESTRING(206.234032029431 -92.0360700742475,194.92038 14.4658)
LINESTRING(-125.67524 -98.00064,-131.974995810398 -38.6145300952198)
LINESTRING(-131.974995810398 -38.6145300952198,-137.6955 15.31113)
LINESTRING(-120.60474 22.63271,-107.63268 -99.71427)
LINESTRING(35.93373 39.66043,52.34093 -114.89627)
LINESTRING(26.19878 -108.47147,10.55734 38.93942)
which looks like this:
I expect to get the LINESTRINGS
LINESTRING(-131.975 -38.61453,-114.06774 -39.02152)
LINESTRING(18.5567 -36.44961,10.55734 38.93942)
LINESTRING(219.84602 -92.53407,206.23403 -92.03607)
as marked with the black circles.
What I'm doing wrong or is it a matter of mathmatical precision? The setup of the segment groups is related to this GIS-SE post.
I use Debian Jessie, QGIS 2.10, PostGIS-2.1 and PostgreSQL 9.4.1
Answer
The issue is that the two geometries don't share common nodes. For instance, A has vertices mid-way that don't exist in B. While they visually overlap, there are tiny round-off errors from the interpolation used in the algorithm to determine differences, and you see unexpected results.
To avoid the round-off errors with overlay operators available with GEOS (PostGIS, Shapely) and JTS, you require exact noding. This can be done in PostGIS with ST_Snap. The pseudo code for the use could be something like:
ST_Overlay(ST_Snap(A, B, tol), ST_Snap(B, A, tol))
which can be adapted to the question to use a tolerance tol of 1e-5, because that's how many decimal places there are with the input datasets:
SELECT ST_AsText((ST_Dump(
ST_SymDifference(
ST_Snap(A, B, tol),
ST_Snap(B, A, tol))
)).geom)
FROM (
SELECT ST_GeometryFromText('MULTILINESTRING(
(219.84602 -92.53407,223.95468 -131.19481),
(212.94651 -27.6126,219.84602 -92.53407),
(-114.06774 -39.02152,-107.63268 -99.71427),
(-125.67524 -98.00064,-131.975 -38.61453),
(-131.975 -38.61453,-114.06774 -39.02152),
(44.44439 -40.51042,52.34093 -114.89627),
(26.19878 -108.47147,18.5567 -36.44961),
(210.25568 -129.89414,206.23403 -92.03607),
(219.84602 -92.53407,206.23403 -92.03607),
(206.23403 -92.03607,194.92038 14.4658),
(-131.975 -38.61453,-137.6955 15.31113),
(-120.60474 22.63271,-114.06774 -39.02152),
(35.93373 39.66043,44.44439 -40.51042),
(18.5567 -36.44961,10.55734 38.93942)
)') A,
ST_GeometryFromText('MULTILINESTRING(
(212.94651 -27.6126,223.95468 -131.19481),
(210.25568 -129.89414,194.92038 14.4658),
(-125.67524 -98.00064,-137.6955 15.31113),
(-120.60474 22.63271,-107.63268 -99.71427),
(35.93373 39.66043,52.34093 -114.89627),
(26.19878 -108.47147,10.55734 38.93942))') B,
1e-5 tol
) f;
st_astext
-----------------------------------------------------
LINESTRING(-131.975 -38.61453,-114.06774 -39.02152)
LINESTRING(219.84602 -92.53407,206.23403 -92.03607)
(2 rows)
(Note that the WKT for A in this question has two horizontal bars, thus two lines are expected here)
No comments:
Post a Comment