[postgis-devel] precision issue with several functions

Thu Oct 10 02:11:01 PDT 2013

I'm well aware of SFCGAL possibilities !

It's like a christmas gift for my phd,
and I'll definitively use it when I can.
For the moment I'm prototyping so I try to avoid it.

Snapping to grid is not a solution exactly for the same reason than
snapping a point to a line doesn't work.
Say the projected point is exactly on the grid, the extremities of the line
are exactly on the grid, but the line does not go exactly trough this
point.
Snapping to grid won't change a anything, and I'll still have the same
precision issue.

Unfortunately I'm having this precision issue while using ST_Split, the
exact issue being that splitting line by point doesn't split if point is
not sufficiently on the line. Here ST_DWithin is of no help, but snapping
line to point works.

For other intersection test I'll use ST_DWithin, or directly BBox test when
possible (bbox are already less precise if I remember right)

Cheers,

Rémi-C

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

> 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;
> _______________________________________________
> 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/abfe1a27/attachment.html>