13

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?

Fabzi
  • 626
  • 5
  • 15
  • Are you able to choose a more granular unit, lets say, millimeters instead of meters? – Paulo Scardine Jan 19 '15 at 16:40
  • @PauloScardine : thanks. Yes, I can easily give away precision if I gain stability. Multipliying my values by 10 makes the trick. But will it work in **all** cases? Shapely continues to work with floats internally. – Fabzi Jan 20 '15 at 08:14

2 Answers2

14

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.

Georgy
  • 12,464
  • 7
  • 65
  • 73
Mike T
  • 41,085
  • 18
  • 152
  • 203
  • Hi @MikeT, thanks for your answer, that's what I was affraid of. I really wish we had a precision model. Or will the trick of working with integers instead of floats (see Paulo's comment above) ensure stability? My use case is more complicated and based on unions: `shapely.ops.polygonize(LinearRing(line1).union(LineString([pt, pr])))` returns different polygons if the two lines intersect or not... _(I wish I could upvote your answer but I havn't enough reputation for that ;)_ – Fabzi Jan 20 '15 at 08:29
  • Can these distance tests be used to see when two linestrings touch - in the sense that the last point of the one is the first one of the other? I have encountered linestrings where the points differ by those characteristic 0.0000000001 (due to float inaccuracies) – Mr_and_Mrs_D Dec 01 '21 at 17:01
3

As of shapely 2.0 (see release notes), a precision model is now available: the set_precision function and the grid_size parameter for overlay functions.

Using your previous example and setting the same arbitrary precision to the geometries.

from shapely import set_precision
from shapely.geometry import LineString, Point

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

print(line1.intersects(LineString([pt, pr])))  # True

Note: I would have expected the set_precision for the pr object to be redundant since all of the involved geometries have been created with an arbitrary precision, but that's not the case...

contains issue

From my experiments, the contains predicates is never true for any precision and any combination of precision. I'm not qualified enough in GEOS/shapely to give a proper explanation.

Still, intuitively I think it has something to do with the types of geometries involved: using the usual buffer technique explained by @MikeT, we actually change the line/point geometry to a polygon which has consequences on the dimensions of the boundaries and interiors of the involved geometries.

jmon12
  • 1,090
  • 8
  • 17