<div dir="ltr">Hey Sandro thanks for the answer :-)<div><br></div><div><b>Before going into details I tested your suggestion (snap line to point) and it solved the problem.</b></div><div><b>Suppressing some digits from coordinates solves also the problem (6860649 -></b><b>649</b><b> )</b></div>
<div><div><br></div><div>Thanks for the details about ST_Snap, effectively the behavior is understandable. (I still think we need a more general snap function, although it would be very difficult to do for non-trivial cases like line to line with no shared vertex but almost shared path).</div>
<div><br></div><div>I don't understand your answer about precision. I'm not saying postgis processing should be exact (as with the CGAL kernel), which is very difficult and would need a rewrite of everything.<br>
<br>
</div><div>So imprecision in computing is OK (PostGIS is not CGAL !).</div><div><br></div><div>But it is terribly wrong to have spatial relation function misbehaving : PostGis is all about spatial relationship !</div><div>
<br></div><div>I feel PostGIS is made for real GIS data : Currently, using the French projection for mainland France (1000km*1000km) and data with 1 meter of precision, PostGIS answer a point is not on a line because of an error of 10 nano meters !</div>
<div><br></div>
<div>It just makes no sense.</div><div><br></div><div>So I totally agree with your conclusion on <a href="http://trac.osgeo.org/postgis/wiki/ToleranceDiscussion" target="_blank">ToleranceDiscussion</a> about spatial relation operator with a precision given by user, and a defined precision for each data.</div>
<div>As a matter of fact the only rightful solution wouldn't be CGAL-like exact computation , but <b>fuzzy estimation</b> (Again, totally in my dream ;-) , and not ready, even in theory).</div><div><br></div><div><br>
</div>
<div><br></div><div>Unfortunately Snapping to grid is not a solution in many cases (including this one).</div><div><br></div><div><br></div><div>Here I sum up the possible workaround I might try (details after):</div><div>
<ul><li>Use snap with precision p when it is possible</li><li>Use pivot to reduce the number of digits in coordinates</li><li>Use ST_Distance with a cap instead of St_Intersects </li><li>Shift position by p in all direction and do spatial relation test</li>
<li>buffering geom by a radius of precision</li></ul><div>And other possible solutions</div><div><ul><li>Use SFCGAL with exact kernel</li><li>artificially augment number of digits in coordinates to reach precision limit of st_intersect (ad x 9 as a prefix of coordinates)</li>
</ul></div><br>AS a temporary workaround I propose a very ineffective and ugly solution (and again just an approximation):</div><div><br></div></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div>
<div>When trying to find relation between <i><b>A(x,y)</b></i> and <i><b>B</b></i> with precision <b style="font-style:italic">p </b>(15 digits in total in my test), do the computing 4 times:</div></div><div><div><br></div>
</div><div><div><i><b>Result = A_choice_function( Relation( A(x-p,y-p) , B ) Relation( A(x-p,y+p) , B ) OR Relation( A(x+p,y-p) , B ) OR Relation( A(x+p,y+p) , B ) )</b></i></div></div><div><div><br></div></div><div><div>
The choice function depending on the kind of behavior we are looking for (in my case, it would be <b>OR</b> , as I prefer to have false positive than missing positive )</div></div></blockquote><div>
<div><br></div><div><br></div><div>Also I will use ST_Distance with a cap on digit instead of ST_Intersects (like keep only 3 digits after coma, if all 0, intersects, else not intersects).<br>It sure will be slower...</div>
<div><br></div><div>Another working workaround : buffering all geometries with a radius of precision. This works also. Again very ineffective, and false for most spatial relation.</div><div><br></div><div><br></div><div>Testing the geometry with your (hidden) function gives an answer very coherent with the problem: </div>
<div><div><br></div></div></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">
<div><div><div>select topology._st_mintolerance('MULTILINESTRING((651006 6860741.8,650880.2 6860649.4,650872.7 6860644.9))');</div></div></div><div><div>2.46986704800001e-08</div></div><div><br></div><div><br></div>
</blockquote><div><div>Cheers,</div></div><div>Rémi-C</div><div class="gmail_extra"><br><br><div class="gmail_quote">2013/10/9 Sandro Santilli <span dir="ltr"><<a href="mailto:strk@keybit.net" target="_blank">strk@keybit.net</a>></span><br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div>On Wed, Oct 09, 2013 at 04:15:33PM +0200, Rémi Cura wrote:<br>
> Hey fellow postgis users,<br>
> I ran into some precision issue.<br>
> (I use latest postgis/postgres/geos release on a ubuntu 12.0.4 :<br>
> postgis_full_version<br>
> POSTGIS="2.1.0 r11822" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6<br>
> March 2012" GDAL="GDAL 1.10.0, released 2013/04/24" LIBXML="2.8.0" TOPOLOGY<br>
> (topology procs from "2.0.3 r11128" need upgrade) RASTER ).<br>
><br>
</div>> *I would appreciate expert answers on this as precision issued introduce**<br>
> (silently) **very incoherent results. *<br>
<div>><br>
><br>
> During a geometry processing, I project a point on a line using the<br>
> ST_CLosestPoint() function.<br>
><br>
> The issue are<br>
> _the projected point is at a residual distance of 2 x 10^-9 meters of the<br>
> line<br>
<br>
</div>Expected, not all the points in a line can be represented.<br>
<div><br>
> _the projected point is not related to the line (ST_Relate FF*FF****)<br>
<br>
</div>For the same reason as the above.<br>
<div><br>
> _snapping the point to the line does nothing<br>
<br>
</div>Snapping only snaps to vertices, not to segments.<br>
If the snapper tried to snap to segments it would still fail<br>
due to the above (can't create points in arbitrary places).<br>
<div><br>
> _snapping the line to the point moves the line (and they relate :<br>
> 0F*FF****).<br>
<br>
</div>Expected, and currently being the best approach I can think of.<br>
<div><br>
> So it appears to me that :<br>
><br>
</div>> - if spatial testing is more precise than 10^-9, ST_ClosestPoint()<br>
<div>> should return something more precise than this<br>
<br>
</div>Cannot. Like there's no way to write the result of 10/3 with a finite<br>
number of digits.<br>
<br>
> - (in the above example, a point projected on a line is disjoint from<br>
> this line)<br>
<br>
It's supposedly as close as it can get.<br>
<br>
> - ST_Snap() does nothing when snapping point on line with a vastly<br>
<div>> sufficient tolerance, but snapping the line to the point moves the line.<br>
> This is at best incoherent.<br>
<br>
</div>It's documented as the behavior: <a href="http://postgis.org/docs/ST_Snap.html" target="_blank">http://postgis.org/docs/ST_Snap.html</a><br>
"Snaps the vertices and segments of a geometry another Geometry's vertices"<br>
<div><br>
> It seems to be not related to size of coordinates (650880,X being the limit<br>
> for a float)<br>
<br>
</div>The limit of a float depends on the magnitude of the number.<br>
I tried to provide an (internal, don't tell anyone) function to compute<br>
the worst (larger) tolerance for a geometry:<br>
<br>
<br>
strk=# select topology._st_mintolerance('POINT(0 0)');<br>
_st_mintolerance<br>
------------------<br>
3.6e-15<br>
(1 row)<br>
<br>
strk=# select topology._st_mintolerance('POINT(10 0)');<br>
_st_mintolerance<br>
------------------<br>
3.6e-14<br>
(1 row)<br>
<br>
So my suggestion is:<br>
<br>
1) Find the min tolerance for your working extent<br>
(local projection help keeping it low)<br>
2) Snap all vertices to a grid of the tolerated size<br>
3) When you project a point on a line, snap the projected<br>
point on the grid and snap the line to the point<br>
<br>
Automating all of these would require encoding a "precision model"<br>
in each and every geometry. Wouldn't be bad. Whant to help ?<br>
Some ideas are here: <a href="http://trac.osgeo.org/postgis/wiki/ToleranceDiscussion" target="_blank">http://trac.osgeo.org/postgis/wiki/ToleranceDiscussion</a><br>
<br>
<br>
--strk;<br>
_______________________________________________<br>
postgis-devel mailing list<br>
<a href="mailto:postgis-devel@lists.osgeo.org" target="_blank">postgis-devel@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-devel" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-devel</a><br>
</blockquote></div><br></div></div>