ST_Within and ST_DWithin disagree

[strk]

=# with inp as ( select '010100000070841C1ED7F9204187A97F8DB8AF0141'::geometry as p, '010200000002000000B6813B20D7F92041F5807988B8AF014100000000D6F9204133333333BBAF0141'::geometry as l 
) select ST_AsText(l) as l, ST_AsText(p) as p, 
  ST_Within(p, l) as within, 
  ST_DWithin(p, l, 0) as dwithin 
from inp;
-[ RECORD 1 ]----------------------------------------------------
l       | LINESTRING(556267.562954 144887.066638,556267 144887.4)
p       | POINT(556267.55881132 144887.069091153)
within  | f
dwithin | t

I suspect ST_Within is right, but the point was obtained by ClosestPoint? so is the best we can do to _find_ a point within a line. Having a too-precise ST_Within doesn't really help us.

Maybe we can keep overly-precise ST_Within and use ST_DWithin for topology, but worth an opinion by Martin Davis...

See #1613 for the consequences of this.

[strk]

Shell I add: ST_Distance(p,l) returns 0, which doesn't help

[strk]

I suspect ST_Distance itself is kind of bogus. lw_dist2d_pt_seg in particular seems to be ending with projecting the point on the segment and then computing the distance ?

[strk]

Yup, I think ST_Distance is doing the projection and in doing so merges the query point with the segment thus ending up with a zero distance. This was likely done for ST_ClosestPoints.

[nicklas]

Yes, you are right that it is doing the reprojection. That was the choice when designing ST_ClosestPoint and ST_ShortestLine.

It is mentioned in the comments in the code I think.

I see the consequences of this.

Both #1502 and the related ticket you fixed lately is related to this.

[strk]

Nicklas: so are you planning to fix this ? Meanwhile I committed some crazy computation of a "smallest distance in float space for a given max ordinate value" which was effective to fix #1613, but doesn't feel safe

[nicklas]

One way to quantify how much the two algorithms disagree is to see how much the point have to be buffered to intersect with the line. The answer seems in this case to be 0.0000000001 units which is far less than the precision of the input.

select st_intersects('010200000002000000B6813B20D7F92041F5807988B8AF014100000000D6F9204133333333BBAF0141'::geometry,st_buffer('010100000070841C1ED7F9204187A97F8DB8AF0141'::geometry, 0.0000000001));

This returns true

I have twisted my head when I was looking at #1502 to understand what is right and what is wrong. I don't know, but to get consistent answers beyond the precision of the input I guess is only possible if the approach is identical (so the error is identical).

Does that make sense?

My conclusion is that we might have to start thinking about how to handle tolerance and precision in generic terms.

From what I think SQL-Server is outputting their answers with far less precision than double precision and in ArcGIS world you tell the system what precision you want to use.

[nicklas]

strk

I would love to fix it, but as I wrote above I think it maybe is a bigger discussion.

The quick fix (which I would hate) is to not project the point when calling ST_Dwithin or ST_Distance.

The reason I wouldn't like that is because the consistence between the results of st_shortestline and st_closestpoint and st_distance and st_dwithin also has a value.

But I see no other fix since the projected point cannot have better precision than it gets. So if st_shortestline really is to give the points from where st_distance uses (which was the idea when writing the function) this is what we get :-(

[strk]

+1 for "global tolerance concept" We'd need this in GEOS too, btw.

I do have that for topology already, but not using it much: #785

[nicklas]

But what to do in 2.0?

Leave it as it is?

I can open a wikipage to discuss what functions that need tolerance and how.

Strk, your comment in #1613 about how the floating point grid looks on different distance from origo in each projection is what it is all about I think and makes it quite complicated.

[strk]

Please go ahead with the wiki page. I'm sure ST_Distance doesn't want a tolerance and yes I do see that ST_ClosestPoint maybe be at a non-zero Distance by doing that. I think we need to take for granted that _adding_ a point on a line can't be on the original line, except rare occasions. Our space is not infinite. I think it is important not to try pretending otherwise.

Yeah the expanding float grid is an interesting concept. May be worth an internal function to compute it. Does sound like a reprojection issue doesn't it ? :)

[strk]

No, not for 2.0, it's not really a bug but a design decision

[strk]

Got hurt again by this with #1789.

[pramsey]

We aren't going to be changing tolerance models this time out...

[strk]

For the record, wikipage is here: http://trac.osgeo.org/postgis/wiki/ToleranceDiscussion