<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<meta name="Generator" content="Microsoft Word 14 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
@font-face
{font-family:Consolas;
panose-1:2 11 6 9 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
span.E-postmall17
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri","sans-serif";
mso-fareast-language:EN-US;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:70.85pt 70.85pt 70.85pt 70.85pt;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="SV" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">Thanks Martin,<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">you confirmed what I suspected. A tolerance would have been nice to use.<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">/Paul<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif"">Från:</span></b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif""> postgis-users [mailto:postgis-users-bounces@lists.osgeo.org]
<b>För </b>Martin Davis<br>
<b>Skickat:</b> den 29 augusti 2019 18:06<br>
<b>Till:</b> PostGIS Users Discussion<br>
<b>Ämne:</b> Re: [postgis-users] st_witin<o:p></o:p></span></p>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">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. <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">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].<o:p></o:p></p>
</div>
<div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">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.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">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).<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">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.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">[1] <a href="https://locationtech.github.io/jts/jts-faq.html#D7">https://locationtech.github.io/jts/jts-faq.html#D7</a><o:p></o:p></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<p class="MsoNormal">On Thu, Aug 29, 2019 at 4:58 AM <<a href="mailto:paul.malm@lfv.se">paul.malm@lfv.se</a>> wrote:<o:p></o:p></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-right:0cm">
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">Hi,</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">I have a layer with lines that I have buffered to a polygon layer, those polygons (not multipolygons) are unioned and containes holes.</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">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.</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">I’ve tried like this:</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF">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)</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">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.</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">I run makevalid on the tables, to be sure</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF;background:#E8F2FE">UPDATE polygonlayer SET the_geom = ST_Makevalid(the_geom) WHERE st_isvalid(the_geom)=false</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF;background:#E8F2FE">UPDATE
</span><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF">linelayer
<span style="background:#E8F2FE">SET the_geom = ST_Makevalid(the_geom) WHERE st_isvalid(the_geom)=false</span></span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US"> </span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">I create a primary key</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF">ALTER TABLE linelayer ADD COLUMN \"pkkey\" serial NOT NULL PRIMARY KEY</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US"> </span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">I reindex the tables, to be sure:</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF">REINDEX TABLE linelayer</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF">REINDEX TABLE polygonlayer</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas"> </span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">I change to LineStrings just to be sure not having several linestings in a MultiLineString
</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF">CREATE TABLE dumpedlines AS SELECT *, (ST_Dump(the_geom)).geom AS the_geom2 FROM lineLayer</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF">ALTER TABLE dumpedlines DROP COLUMN IF EXISTS the_geom</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF">ALTER TABLE dumpedlines RENAME COLUMN the_geom2 TO the_geom</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF">ALTER TABLE dumpedlines ALTER COLUMN the_geom TYPE geometry(LineString, 32631)</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:10.0pt;font-family:Consolas;color:#BB005E"> </span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">Then I try to delete the lines outside the buffered polygons:</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US" style="font-size:10.0pt;font-family:Consolas;color:#2A00FF">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))</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US"> </span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">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 _<i>ST</i>_within but with the same result.</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US"> </span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">Any ideas?</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">Thanks in advance,</span><o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="EN-US">Paul</span><o:p></o:p></p>
</div>
</div>
<p class="MsoNormal">_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/postgis-users" target="_blank">https://lists.osgeo.org/mailman/listinfo/postgis-users</a><o:p></o:p></p>
</blockquote>
</div>
</div>
</div>
</body>
</html>