Shapely Polygons that should touch don't

Question

Here is some Shapely code that creates three square polygons, p1, p2, and p3. p2 is positioned immediately to the right of p1, and p3 is positioned immediately underneath.

The problem is that Shapely tells me p1 and p2 don't touch, whereas p1 and p3 do. I can't see what is going wrong here.

from shapely.geometry import Polygon

DELTA = 0.2

def polygonFromPoint(p):
    return Polygon([(p[0]-DELTA*0.5, p[1]-DELTA*0.5),
                    (p[0]-DELTA*0.5, p[1]+DELTA*0.5),
                    (p[0]+DELTA*0.5, p[1]+DELTA*0.5),
                    (p[0]+DELTA*0.5, p[1]-DELTA*0.5)])

p1 = polygonFromPoint([-118.8,35.0])
p2 = polygonFromPoint([-118.6,35.0])
p3 = polygonFromPoint([-118.8,34.8])

print(p1)
print(p2)
print(p3)

print(p1.overlaps(p2), p1.intersects(p2), p1.crosses(p2), p1.contains(p2),
      p1.disjoint(p2), p1.touches(p2))
print(p1.overlaps(p3), p1.intersects(p3), p1.crosses(p3), p1.contains(p3),
      p1.disjoint(p3), p1.touches(p3))

Running this produces the following output:

POLYGON ((-118.9 34.9, -118.9 35.1, -118.7 35.1, -118.7 34.9, -118.9 34.9))
POLYGON ((-118.7 34.9, -118.7 35.1, -118.5 35.1, -118.5 34.9, -118.7 34.9))
POLYGON ((-118.9 34.7, -118.9 34.9, -118.7 34.9, -118.7 34.7, -118.9 34.7))
False False False False True False
False True False False False True

Which shows that Shapely thinks p1 and p2 don't intersect or touch, whereas p1 and p3 intersect and touch.

EDIT: As Gilles-Philippe Paillé and others remarked, this is a precision problem with the polygon coordinates. Using the following tweak fixes the issue in this case:

def polygonFromPoint(p):
    return Polygon( [(round(p[0]-DELTA*0.5,1), round(p[1]-DELTA*0.5,1)),
                        (round(p[0]-DELTA*0.5,1), round(p[1]+DELTA*0.5,1)),
                        (round(p[0]+DELTA*0.5,1), round(p[1]+DELTA*0.5,1)),
                        (round(p[0]+DELTA*0.5,1), round(p[1]-DELTA*0.5,1))] )

Computers can only represent an approximation of most floating point numbers, so you cannot expect comparisons of values calculated from them to be exact. You may be able to work around this limitation by using only integer coordinates. — martineau, Jan 27 '20 at 18:54
See related [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) — martineau, Jan 27 '20 at 21:03
Thanks - your point is well taken! However, I am not making any comparisons in my code, I am asking the Shapely library to do whatever math is required to determine if the polygons are touching etc., based on the precise coordinates given. Should it not take care of handling the precision of its own calculations? Edit: perhaps I'm not giving it as precise coordinates as I think, as Paille shows below. — Julian, Jan 27 '20 at 21:31
Yes, I understand that…and it sounds to me like it's due to a bug in Shapely's implementation of `touches()` — it's a Comp Sci 101 kind of mistake (and your mention of "precise coordinates" makes me suspect you may still not get that there's no such thing on a computer, either). — martineau, Jan 27 '20 at 21:37
No, you are correct. See my revised question, where I round to one decimal place and it works. It's curious that Shapely, when printing the polygon, doesn't show the full precision of the numbers it is using, otherwise it would have been immediately obvious that there was a precision issue :-) — Julian, Jan 27 '20 at 21:44
Even if it was the "full precision" of the number, there's a related problem that it's often not possible to correctly print the decimal representation of the binary floating point number stored internally — so whatever you see printed is also an approximation (of something that may itself be an approximation of something else `;¬)`. — martineau, Jan 27 '20 at 21:49
I'd like to add one more thing after seeing your edit — a better way would be to see if the difference between the two values was something less than a relatively small value, like say `0.00001`. Using `round()` makes me a little nervous because it may involve doing more floating-point math. — martineau, Jan 27 '20 at 22:01
Does this answer your question? [How to deal with rounding errors in Shapely](https://stackoverflow.com/questions/28028910/how-to-deal-with-rounding-errors-in-shapely) — Georgy, Jan 27 '20 at 22:51

Gilles-Philippe Paillé · Accepted Answer · 2020-01-27T19:01:58.687

Even if the string representation of the polygon show that the coordinates are the same, the underlying floating point representation is not exactly what is printed and may include some imprecision. Using your coordinates and the same calculation, I get:

DELTA = 0.5

a = -118.6 - 0.2 * DELTA
b = -118.8 + 0.2 * DELTA
print(a)
print(b)
print(a <= b)

a = 35.0 - 0.2 * DELTA
b = 34.8 + 0.2 * DELTA
print(a)
print(b)
print(a <= b)

which gives the following output:

-118.69999999999999
-118.7
False
34.9
34.9
True

Shapely Polygons that should touch don't

1 Answers1