[postgis-devel] Behavior of LineString-LineString ST_Intersection and LineString-Polygon ST_Difference

[Stefan Keller]

Hi,

The goal is to get all linestring features which are completely
outside of polygons including those segments up to their intersection
with polygons. Think of "all streets that don't cross forests" (see my
other post in postgis-users).

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

To me the following query would be close to the solution - if all
these irregularities reported below would'nt exist:

SELECT DISTINCT ST_Astext(ST_Difference(ln1.geom, ln2.geom))
FROM
  lines ln1,
  (
   SELECT ST_Intersection(ln.geom, po.geom) as geom
   FROM lines ln, polygons po
   WHERE (ST_Intersects(ln.geom, po.geom) OR ST_Contains(po.geom, ln.geom))
   AND NOT ST_Touches(ln.geom, po.geom)
  ) ln2
"LINESTRING(2 3,3 4,5 4,6 4)"     -- (*) still is not clipped by PO2.
Could be merged later on with (**)
"MULTILINESTRING((2 1,9.58018867924528
1.75801886792453),(9.58018867924528 1.75801886792453,12 2))"
"LINESTRING(4 4,5 4,6.2 4)"       -- (**) ok, at least a segment.
Could be merged later on with (*)
"LINESTRING(5 1,7 1)"             -- ok
"LINESTRING(2 1,12 2)"            -- not clipped by PO2, processed
into MULTILINESTRING!
"LINESTRING(1 1,2 3,3 4,5 4,6 4)" -- should'nt belong here
"LINESTRING(4 4,5 4,7 4,10 2)"    -- should'nt belong here

So, let's begin with some tests about LineString-LineString-intersection:

-- ST_Intersection LineString-LineString:
SELECT ST_Astext(ST_Intersection(ln1.geom, ln2.geom))
FROM lines ln1, lines ln2
WHERE ST_Intersects(ln1.geom, ln2.geom)
"MULTILINESTRING((1 1,2 3),(2 3,3 4),(3 4,5 4),(5 4,6 4))"
"MULTILINESTRING((4 4,5 4),(5 4,6 4))"
"LINESTRING(2 1,12 2)"
"MULTILINESTRING((4 4,5 4),(5 4,6 4))"
"MULTILINESTRING((4 4,5 4),(5 4,7 4),(7 4,10 2))"
"LINESTRING(5 1,7 1)" --???
=> Why those 'pseudo' MULTILINESTRINGs?
=> Why are *all* lines returned - even L4 "(5 1,7 1)" which does *not*
intersect?

Now the intersection and difference between LineString and Polygons:

-- ST_Intersection Line-Polygon:
SELECT ST_Astext(ST_Intersection(ln.geom, po.geom))
FROM lines ln, polygons po
WHERE (ST_Intersects(ln.geom, po.geom) OR ST_Contains(po.geom, ln.geom))
AND NOT ST_Touches(ln.geom, po.geom)
"LINESTRING(1 1,2 3)"
"LINESTRING(7.66037735849057 1.56603773584906,11 1.9)"
"LINESTRING(6.2 4,7 4,10 2)"
=> OK.

-- ST_Difference LineString-Polygon:
SELECT ST_Astext(ST_Difference(ln.geom, po.geom))
FROM lines ln, polygons po
WHERE ST_Intersects(ln.geom, po.geom)
"LINESTRING(2 3,3 4,5 4,6 4)" -- L1-PO1
"LINESTRING(2 1,12 2)" --???
"MULTILINESTRING((2 1,7.66037735849057 1.56603773584906),(11 1.9,12 2))" -L3-PO2
"LINESTRING(4 4,5 4,6.2 4)" -- L2-PO2
=> Ok, except superfluous L3 "(2 1,12 2)": intersects with PO2 but was
already processed!?

-- ST_SymDifference LineString-Polygon
SELECT ST_Astext(ST_SymDifference(ln.geom, po.geom))
FROM lines ln, polygons po
WHERE NOT ST_Intersects(ln.geom, po.geom)
"GEOMETRYCOLLECTION(LINESTRING(1 1,2 3,3 4,5 4,6 4),POLYGON((8 1,5
6,11 6,11 1,8 1)))"
"GEOMETRYCOLLECTION(LINESTRING(4 4,5 4,7 4,10 2),POLYGON((0 0,0 3,2
3,2 0,0 0)))"
"GEOMETRYCOLLECTION(LINESTRING(5 1,7 1),POLYGON((0 0,0 3,2 3,2 0,0 0)))"
"GEOMETRYCOLLECTION(LINESTRING(5 1,7 1),POLYGON((8 1,5 6,11 6,11 1,8 1)))"
=> What the heck is this...?

SymDifference is defined as "Returns a geometry that represents the
portions of A and B that do not intersect." It seems that there are
only existing nodes returned and there is no "real" intersection
calculation behind (look at the fractional intersection points from
ST_Difference)...
Regina backs this up here:
http://www.bostongis.com/blog/index.php?/archives/48-Spatial-Differences.html

So I think the ST_SymDifference LineString-Polygon is wrong. I'd
expect something like all remaining geometries calculated from
ST_Difference(ln.geom, po.geom), which would almost solve my initial
problem statement...

Any comments on this?

Yours, S.