[postgis-devel] precision issue with several functions

Thu Oct 10 01:39:22 PDT 2013

Hey Sandro thanks for the answer :-)

*Before going into details I tested your suggestion (snap line to point)
and it solved the problem.*
*Suppressing some digits from coordinates solves also the problem (6860649
->**649** )*

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).

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.

So imprecision in computing is OK (PostGIS is not CGAL !).

But it is terribly wrong to have spatial relation function misbehaving :
PostGis is all about spatial relationship !

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 !

It just makes no sense.

So I totally agree with your conclusion on
ToleranceDiscussion<http://trac.osgeo.org/postgis/wiki/ToleranceDiscussion>
about
spatial relation operator with a precision given by user, and a defined
precision for each data.
As a matter of fact the only rightful solution wouldn't be CGAL-like exact
computation , but *fuzzy estimation* (Again, totally in my dream ;-) , and
not ready, even in theory).

Unfortunately Snapping to grid is not a solution in many cases (including
this one).

Here I sum up the possible workaround I might try (details after):

   - Use snap with precision p when it is possible
   - Use pivot to reduce the number of digits in coordinates
   - Use ST_Distance with a cap instead of St_Intersects
   - Shift position by p in all direction and do spatial relation test
   - buffering geom by a radius of precision

And other possible solutions

   - Use SFCGAL with exact kernel
   - artificially augment number of digits in coordinates to reach
   precision limit of st_intersect (ad x 9 as a  prefix of coordinates)

AS a temporary workaround I propose a very ineffective and ugly solution
(and again just an approximation):

When trying to find relation between *A(x,y)* and *B* with precision *p *(15
digits in total in my test), do the computing 4 times:

*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
)  )*

The choice function depending on the kind of behavior we are looking for
(in my case, it would be *OR* , as I prefer to have false positive than
missing positive )

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).
It sure will be slower...

Another working workaround : buffering all geometries with a radius of
precision. This works also. Again very ineffective, and false for most
spatial relation.

Testing the geometry with your (hidden) function gives an answer very
coherent with the problem:

select topology._st_mintolerance('MULTILINESTRING((651006
6860741.8,650880.2 6860649.4,650872.7 6860644.9))');
2.46986704800001e-08

Cheers,
Rémi-C

2013/10/9 Sandro Santilli <strk at keybit.net>

> On Wed, Oct 09, 2013 at 04:15:33PM +0200, Rémi Cura wrote:
> > Hey fellow postgis users,
> > I ran into some precision issue.
> > (I use latest postgis/postgres/geos release on a ubuntu 12.0.4 :
> > postgis_full_version
> > POSTGIS="2.1.0 r11822" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6
> > March 2012" GDAL="GDAL 1.10.0, released 2013/04/24" LIBXML="2.8.0"
> TOPOLOGY
> > (topology procs from "2.0.3 r11128" need upgrade) RASTER ).
> >
> > *I would appreciate expert answers on this as precision issued
> introduce**
> >  (silently) **very incoherent results. *
> >
> >
> > During a geometry processing, I project a point on a line using the
> > ST_CLosestPoint() function.
> >
> > The issue are
> > _the projected point is at a residual distance of 2 x 10^-9 meters of the
> > line
>
> Expected, not all the points in a line can be represented.
>
> > _the projected point is not related to the line (ST_Relate FF*FF****)
>
> For the same reason as the above.
>
> > _snapping the point to the line does nothing
>
> Snapping only snaps to vertices, not to segments.
> If the snapper tried to snap to segments it would still fail
> due to the above (can't create points in arbitrary places).
>
> > _snapping the line to the point moves the line (and they relate :
> > 0F*FF****).
>
> Expected, and currently being the best approach I can think of.
>
> > So it appears to me that :
> >
> >    - if spatial testing is more precise than 10^-9, ST_ClosestPoint()
> >    should return something more precise than this
>
> Cannot. Like there's no way to write the result of 10/3 with a finite
> number of digits.
>
> >    - (in the above example, a point projected on a  line is disjoint from
> >       this line)
>
> It's supposedly as close as it can get.
>
> >    - ST_Snap() does nothing when snapping point on line with a vastly
> >    sufficient tolerance, but snapping the line to the point moves the
> line.
> >    This is at best incoherent.
>
> It's documented as the behavior: http://postgis.org/docs/ST_Snap.html
> "Snaps the vertices and segments of a geometry another Geometry's vertices"
>
> > It seems to be not related to size of coordinates (650880,X being the
> limit
> > for a float)
>
> The limit of a float depends on the magnitude of the number.
> I tried to provide an (internal, don't tell anyone) function to compute
> the worst (larger) tolerance for a geometry:
>
>
>  strk=# select topology._st_mintolerance('POINT(0 0)');
>   _st_mintolerance
>  ------------------
>            3.6e-15
>  (1 row)
>
>  strk=# select topology._st_mintolerance('POINT(10 0)');
>   _st_mintolerance
>  ------------------
>            3.6e-14
>  (1 row)
>
> So my suggestion is:
>
> 1) Find the min tolerance for your working extent
>    (local projection help keeping it low)
> 2) Snap all vertices to a grid of the tolerated size
> 3) When you project a point on a line, snap the projected
>    point on the grid and snap the line to the point
>
> Automating all of these would require encoding a "precision model"
> in each and every geometry. Wouldn't be bad. Whant to help ?
> Some ideas are here:
> http://trac.osgeo.org/postgis/wiki/ToleranceDiscussion
>
>
> --strk;
> _______________________________________________
> postgis-devel mailing list
> postgis-devel at lists.osgeo.org
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-devel/attachments/20131010/a93c1350/attachment.html>