How to deal with rounding errors in Shapely

Question

I have a case which is based on projecting a point on a line and then separate this line on it. My use case is slightly more complicated, but my problem can be reproduced with the following code:

from shapely import *
line1 = LineString([(1,1.2), (2,2), (3, 2.), (4,1.2)])
pt = Point(2.5, 1.2)
pr = line1.interpolate(line1.project(pt))

By construction, "pr" should be on line1 and their intersection too:

line1.contains(pr)
line1.intersects(LineString([pt, pr]))

prints two times "True". But changing the input coordinates slightly brakes the workflow:

from shapely import *
line1 = LineString([(1,1.2), (2,2), (3, 2.3), (4,1.2)])
pt = Point(2.5, 1.2)
pr = line1.interpolate(line1.project(pt))
line1.contains(pr)
line1.intersects(LineString([pt, pr]))

prints "False".

I understand the floating precision problem behind this, but does that mean that I can never test for points being on lines? When I construct a line based on a list of points, can I be sure that at least all the "construction" points will be on the line?

Answer 1

Fundamentally, a precision model is needed, and there are various plans to implement this into GEOS at some time (don't hold your breath, as this has been under discussion for several years).

Otherwise, the options are distance-based tests (recommended) or more expensive buffer-based techniques by a small adjustment (see machine epsilon):

from shapely.geometry import LineString, Point

line1 = LineString([(1, 1.2), (2, 2), (3, 2.3), (4, 1.2)])
pt = Point(2.5, 1.2)
pr = line1.interpolate(line1.project(pt))

# Distance based
print(line1.distance(pr) == 0.0)  # True

# Buffer based
EPS = 1.2e-16
print(line1.buffer(EPS).contains(pr))  # True
print(line1.buffer(EPS).intersects(LineString([pt, pr])))  # True

You can also chain cheaper and expensive tests using an or operator, for example:

print(line1.contains(pr) or line1.buffer(EPS).contains(pr))

which only runs the second and more expensive test if the first one returns False.