[postgis-users] st_witin

[Martin Davis]

My guess is this is caused by numerical precision issues.  Due to numerical
rounding, it is the case that the intersection  of line L with polygon P
does NOT necessarily lie within P.  I.e.  ST_Within( ST_Intersection( L,
P), P )  may or may NOT be true.

In general, due to numerical rounding the overlay operations are not fully
consistent with the spatial predicates.  There's some more information
about this in the JTS FAQ [1].

The reason for this is that in general it is impossible to accurately
represent the intersection point of two line segments defined by floating
point numbers using floating point numbers of the same precision.  So the
result of ST_Intersection may contain a line segment with an endpoint that
falls outside the polygon.  This situation is evaluated precisely by
ST_Within, which thus returns false.

For your case, one fix is to create a new table for the initial
intersections.  This should not contain any lines which did not intersect
the polygon.  (You may also wish to filter out intersection results which
are empty or of very short length, since these can theoretically occur).

In the future it may be possible that PostGIS provides spatial predicates
which accept a distance tolerance, which should allow this issue to be
handled more directly.

[1] https://locationtech.github.io/jts/jts-faq.html#D7

On Thu, Aug 29, 2019 at 4:58 AM <paul.malm at lfv.se> wrote:

> Hi,
>
> I have a layer with lines that I have buffered to a polygon layer, those
> polygons (not multipolygons) are unioned and containes holes.
>
> I would like to create a line layer (from another line layer) with the
> line parts that are within buffered area in the polygon layer.
>
> I’ve tried like this:
>
> update linelayer b set the_geom = ST_MULTI(ST_Intersection(b.the_geom,
> p.the_geom)) FROM polygonlayer p WHERE ST_Intersects(b.the_geom, p.the_geom)
>
> This leaves me with the line parts inside the buffered area and all lines
> that had no intersection with the buffered polygons. That’s ok, but now I
> have to erase the lines that had no intersection with the polygons.
>
> I run makevalid on the tables, to be sure
>
> UPDATE polygonlayer SET the_geom = ST_Makevalid(the_geom) WHERE
> st_isvalid(the_geom)=false
>
> UPDATE linelayer SET the_geom = ST_Makevalid(the_geom) WHERE
> st_isvalid(the_geom)=false
>
>
>
> I create a primary key
>
> ALTER TABLE linelayer ADD COLUMN \"pkkey\" serial NOT NULL PRIMARY KEY
>
>
>
> I reindex the tables, to be sure:
>
> REINDEX TABLE linelayer
>
> REINDEX TABLE polygonlayer
>
>
>
> I change to LineStrings just to be sure not having several linestings in a
> MultiLineString
>
> CREATE TABLE dumpedlines AS SELECT *, (ST_Dump(the_geom)).geom AS
> the_geom2 FROM lineLayer
>
> ALTER TABLE dumpedlines DROP COLUMN IF EXISTS the_geom
>
> ALTER TABLE dumpedlines RENAME COLUMN the_geom2 TO the_geom
>
> ALTER TABLE dumpedlines ALTER COLUMN the_geom TYPE geometry(LineString,
> 32631)
>
>
>
> Then I try to delete the lines outside the buffered polygons:
>
> Delete from dumpedlines b WHERE b.pkkey NOT IN (SELECT b.pkkey FROM
> dumpedlines b, polygon layer c WHERE ST_within(b.the_geom, c.the_geom))
>
>
>
> The result is not correct according to me, all the lines outside the
> buffered polygons are erased, but also SOME lines inside the buffered
> polygons (that are already cut by the polygons). I’ve also tried with _
> *ST*_within but with the same result.
>
>
>
> Any ideas?
>
> Thanks in advance,
>
> Paul
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20190829/7c601242/attachment.html>