[postgis-users] Retrieve the portion of line external to polygon

[strk]

On Wed, Feb 16, 2011 at 07:08:03PM +0100, strk wrote:
> On Wed, Feb 16, 2011 at 06:47:53PM +0100, Andrea Peri wrote:
> > Look this simple example
> > 
> > the difference between the line and the same line intersected with the
> > polygon is equal to the line :)
> > 
> > select 1, ST_Difference(ST_GeomFromText('LINESTRING(10.9 2, 11
> > 8)'),ST_Intersection(ST_GeomFromText('LINESTRING(10.9 2, 11 8)'),
> > ST_GeomFromText('POLYGON((12 2, 10 10, 20 21, 21 5, 12 2))')))
> > 
> > the trick is that in the intersection point between polygon and line, the
> > vertex added move the line so the difference fail to
> > remove the internal portion.
> 
> Ok, I think I know what's going on.
> The original linestring isn't noded with the polygon boundary.
> When computing the intersection, a node is added.
> Such a node will NOT fall on the original line anymore, due
> to precision constraints.
> 
> See this:
>  CREATE TABLE a AS SELECT
>    'LINESTRING(10.9 2, 11 8)'::geometry as g;
>  CREATE TABLE b AS SELECT
>    'POLYGON((12 2, 10 10, 20 21, 21 5, 12 2))'::geometry as g;
>  CREATE TABLE c AS SELECT
>    ST_Intersection(a.g, b.g) FROM a, b;
>  SELECT ST_Covers(a.g, ST_StartPoint(c.g)) FROM a,c; -- false
>  SELECT ST_Covers(a.g,   ST_EndPoint(c.g)) FROM a,c; -- true
> 
> What you could do is _snap_ the original line to the intersection,
> or to node the input before proceeding.

For a working example:
select st_covers(a.g, st_startpoint(st_snap(a.g,c.g,1e-10))) from a,c;

--strk;

  ()   Free GIS & Flash consultant/developer
  /\   http://strk.keybit.net/services.html