Covers predicate incorrect for Line/Point due to geometric robustness

[cuteDen-ECNU]

Consider the following code:

from shapely.geometry import LineString, Point

a = LineString([(1, 0), (0, 2)])

b = Point(0.9, 0.2)

covers = a.covers(b)
# expected - {True}; actual - {False}

a1 = LineString([(10, 0), (0, 20)])

b1 = Point(9, 2)

covers = a1.covers(b1)
# True

Line a should cover point b because the distance between point b and line a is 0.
But Shapely doesn't consider line a covers point b.
The issue cannot be reproduced with integers only, indicating a potential floating-point issue. After all coordinate points are magnified ten times, the results are in line with expectations.

[dr-jts]

This is indeed a bug. Ironically it's due to the high-precision evaluation of the line-point orientation predicate.

Not sure what a suitable fix is at the moment. Likely it will involve some sort of precision rounding, but that needs to be done carefully to ensure there are no other impacts.

See also https://trac.osgeo.org/postgis/ticket/5563

[cuteDen-ECNU]

Thank you for your prompt response. It's greatly appreciated.
Could you help me see whether this potential bug has the same core reason as it?
Once again, thank you for your expertise and assistance in addressing this bug. :·)

[dr-jts]

Could you help me see whether this potential bug has the same core reason as it?

The PostGIS ST_3DIntersects issue has the same basic cause: the inaccuracy of finite-precision floating point arithmetic. But it's much less surprising in that case, since there's a lot of arithmetic being carried out, so the likelihood of 2 lines NOT intersecting in 3D is very high.

The bottom line is that you can't rely on topological relationships being preserved through floating-point operations. This is the essence of the well-known issue of geometrical robustness failures. The solution is usually either to use exact arithmetic (often impractical) or use tolerance values (not a panacea, but very often works).

[cuteDen-ECNU]

Thank you for your prompt response.
The explanation of geometrical robustness failures is very clear.
The 'tolerance values' recommendation is very inspiring to me. 😊

[dr-jts]

The 'tolerance values' recommendation is very inspiring to me. 😊

I'm working on a new implementation of spatial predicate evaluation that should allow for a tolerance value. So stay tuned!

[cuteDen-ECNU]

The tolerance value for predicate evaluation sounds good!
We are excited about the new ideas.

[strk]

Could this be the same problem reported in https://lists.osgeo.org/pipermail/geos-devel/2024-March/011012.html ?

[dr-jts]

Could this be the same problem reported in https://lists.osgeo.org/pipermail/geos-devel/2024-March/011012.html ?

Yes, it's the same issue. In this case, the intersection of the two lines is computed as being equal to the endpoint of one of the lines, due to roundoff error. So the relate result for two lines is inconsistent with the result for the boundary points. Unfortunately the new relate algorithm doesn't fix this.

[Scheggetta]

This exact thing is happening also with GEOS C++ API v3.12.1 using within instead of covers, but I guess it is related to the same issue.

Is there currently any fast and easy workaround? I thought about these two solutions:

• isWithinDistance method of Geometry class
• using within on a buffered LineString

I think the first workaround is more efficient, but could you give me your opinion?

[dr-jts]

Is there currently any fast and easy workaround? I thought about these two solutions:

• isWithinDistance method of Geometry class
• using within on a buffered LineString

I think the first workaround is more efficient, but could you give me your opinion?

Yes, using isWithinDistance is more efficient when it can be used instead of buffer for proximity queries.

That's a good work-around for this issue. Perhaps in the future RelateNG will support a distance tolerance, which will have the same effect.

[dr-jts]

One possible for for this is to use FP computation in the Orientation predicate, and use the result if the point and line are collinear; otherwise use the DD computation for more accuracy. This changes a whole bunch of unit test results, however (which may not indicate a problem, just different results).