6

While working with shapely I faced a strange issue. There're 2 points p1 and p2, where the first one belongs to polygon and the second is not. When I tried to find intersection between the LineString containing these 2 points as endpoints with boundary lines of the Polygon, I received message that no intersection was found. I wonder, how is it possible?

from shapely.geometry import Polygon as SPolygon, Point, LineString

p1 = Point(5.414213562373095, 2.585786437626905)
p2 = Point(15.17279752753168, -7.172797527531679)

l = LineString([p1, p2])

l1 = LineString([(2, 2), (2, 6)])
l2 = LineString([(2, 6), (6, 6)])
l3 = LineString([(6, 6), (6, 2)])
l4 = LineString([(6, 2), (2, 2)])


sp = SPolygon([(2, 2), (2, 6), (6, 6), (6, 2)])

print "Polygon contains p1:", sp.contains(p1)
print "Polygon contains p2:", sp.contains(p2)

for i, line in enumerate((l1, l2, l3, l4)):
    res = l.intersects(line)
    print "Line {0} intersects l1: {1}".format(i, res) 

And here are output:

Polygon contains p1: True
Polygon contains p2: False
Line 0 intersects l1: False
Line 1 intersects l1: False
Line 2 intersects l1: False
Line 3 intersects l1: False
Georgy
  • 12,464
  • 7
  • 65
  • 73
SliceOfTuna
  • 81
  • 1
  • 4
  • 1
    If I am not mistaken, your line `l` passes right through a corner of the polygon. My guess would be that numeric inaccuracies cause it to miss both line segments meeting in that corner, even though it should hit either (or both!) of them. – Ture Pålsson Apr 23 '19 at 07:14
  • 1
    That's weird. It should work for at least one line. `l.distance(l3)` gives `0.0` and `l.distance(l4)` gives `6.435738891173273e-17`. Here is a close-up: https://i.stack.imgur.com/Foa7j.png Consider opening an issue on GitHub: https://github.com/Toblerity/Shapely/issues – Georgy May 31 '19 at 12:01

1 Answers1

3

I changed l.intersects(line) to l.intersection(line) and I did get an intersection at

LINESTRING (6 6, 6 2) at Point (6, 2) 

Not sure why .instersects() is acting differently.

Then I rounded p1 and p2

p1 = Point(round(5.414213562373095, 2), round(2.585786437626905, 2))
p2 = Point(round(15.17279752753168, 2), round(-7.172797527531679, 2))

And I got 2 intersections at

LINESTRING (6 6, 6 2) at POINT (6 2)
LINESTRING (6 2, 2 2) at POINT (6 2)

This fix also worked with .intersects() (2 True's)

Shapely can be a bit touchy with floating point precision, and I usually can fix an issue by rounding. Though this may not be acceptable for you.

Matt M
  • 691
  • 2
  • 6
  • 17