[postgis-devel] precision issue with several functions

[Sandro Santilli]

On Thu, Oct 10, 2013 at 10:39:22AM +0200, Rémi Cura wrote:

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

Actually 2.1.0 adds CGAL backend to PostGIS, you might want to 
give it a try :)

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

In this case we saw they aren't misbehaving, just not behaving as you
thought they did...

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

Relation matrix is like that, but you can use ST_DWithin to test for
a point being within a given distance from a line.

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

That's a collegial research, not just my conclusions.

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

Agreed

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

Why ?

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

...

>    - Use ST_Distance with a cap instead of St_Intersects

ST_DWithin is what you're looking for.

>    - Shift position by p in all direction and do spatial relation test

You'd very hardly find the point _on_ the line.
It might even be impossible to put it on the line, unless it's on an
existing vertex.

> 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

So does an ST_DWithin() with 2.5e-8 tolerance give the desired result ?

--strk;

> 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;