[postgis-devel] precision issue with several functions

[Sandro Santilli]

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;