3

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 A

Segment group B

enter image description here

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:

Segment group C

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.

Expected group

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

3

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)

  • Exists a decomposition operation of "nearly" overlapping MULTILINESTRINGS into "node matching" arc's in the PostGIS/ GEOS/JTS biotop. Things like the CLEAN operation in the old ArcInfo versions? – huckfinn Mar 2 '16 at 1:30
  • @huckfinn it looks like ArcInfo's CLEAN is similar, and even supports a tolerance parameter. – Mike T Mar 2 '16 at 1:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.