[postgis-users] How to calculate the symmetric difference of linestrings and polygons (i.e. linestrings which do not intersect polygons)

[Stefan Keller]

I think I got a more or less nice solution:

-- Some test data:
DROP TABLE IF EXISTS lines;
CREATE TABLE lines (geom geometry);
INSERT INTO lines(geom)
  SELECT       'LINESTRING(1 1, 2 3, 3 4, 5 4, 6 4)'::geometry -- L1
  UNION SELECT 'LINESTRING(4 4, 5 4, 7 4, 10 2)'::geometry -- L2
  UNION SELECT 'LINESTRING(2 1, 12 2)'::geometry -- L3
  UNION SELECT 'LINESTRING(5 1, 7 1)'::geometry -- L4

DROP TABLE IF EXISTS polygons;
CREATE TABLE polygons (geom geometry);
INSERT INTO polygons(geom)
  SELECT 'POLYGON((0 0, 2 0, 2 3, 0 3, 0 0))'::geometry -- PO1
  UNION
  SELECT 'POLYGON((8 1, 11 1, 11 6, 5 6, 8 1))'::geometry  -- PO2

Note that L1 and L2 intersect each other, L1 intersects PO1+PO2, L2
intersects PO2, L3 touches PO1 and goes through PO2, and L4 is outside
PO1 + PO2.

-- First get all intersection lines
-- (decomposing multilinestring result from ST_Difference):
SELECT ST_AsText( (ST_Dump( ST_Difference(ln.geom, po.multigeom) )).geom )
FROM lines ln, -- (collecting polygon to multipolygon)
  (SELECT ST_Collect(po.geom) as multigeom FROM polygons po) po
WHERE ST_Intersects(ln.geom, po.multigeom)
UNION
-- Second, add those line(s) which don't intersect with any (multi-)polygon:
SELECT ST_AsText(ln.geom) as geom
FROM lines ln,
  (SELECT ST_Collect(po.geom) as multigeom FROM polygons po) po
WHERE NOT ST_Intersects(ln.geom, po.multigeom)

Its without support of multipolygon input and my questions regarding
ST_Difference and ST_SymDifference (at postgis-dev) behavior remain.

Yours, S.

2011/1/31 Stefan Keller <sfkeller at gmail.com>:
> Just FYI:
>
> I defined some test queries over at postgis-dev mailinglist with the
> thread "Behavior of LineString-LineString ST_Intersection and
> LineString-Polygon ST_Difference".
>
> Yours, S.
>
> 2011/1/30 Stefan Keller <sfkeller at gmail.com>:
>> Hi,
>>
>> Given streets linestrings and forest polygons, think of all streets
>> segments which are *not* inside forest areas.
>> How can I calculate all segments from linestrings which do not
>> intersect polygons?
>>
>> Any ideas?
>>
>> - S.
>>
>> Attempt 1: ST_SymDifference sounds good! But even this equivalent
>> produces nonsense:
>> # SELECT ST_AsText(
>> ST_CollectionExtract(ST_Difference(ST_Union(ls.way,po.way),
>> ST_Intersection(ls.way,po.way)),2) )
>> FROM linestring ls, polygon po
>>
>> Attempt 2:
>> # SELECT ST_AsText(ST_Intersection(ls.way, po.way)) FROM linestring
>> ls, polygon po
>> gives me all streets *inside* forests. But I'd like get the
>> "remainder" in between...
>>
>> Attempt 3:
>> # SELECT ST_AsText(ST_Difference(ls.way, po.way)) FROM linestring ls, polygon po
>> gives me two displaced lines for each street (because the resultse
>> contains the difference between one line and one polygon), Thesse two
>> could be intersected - but ST_Intersection of two overlapping lines
>> (with same nodes except end nodes) does'nt seem to work.
>>
>> Attempt 4:
>> In a very desperate idea, I calculcated the extent of both tables in
>> order generate a huge boundary polygon with has as many wholes as
>> there are polygons. Then I could do an ST_Intersect of this big swiss
>> cheese polygon with the lines - but no luck.
>>
>> Attempt 5:
>> I finally tried to dump all points from Attempt 2 (intersection) to
>> get hold of the boarder points where lines and polygons cross.
>> But when tried to split a line at a point and to sort out all lines
>> outside polygons things got again complicated..
>>
>