Polygon Intersection with Line | Python Shapely

Question

I have been trying to use shapely to find the intersection of a line and a polygon, but I'm having issues with some floating point numbers.

Example code:

polygon = [(4.0, -2.0), (5.0, -2.0), (4.0, -3.0), (3.0, -3.0), (4.0, -2.0)]
shapely_poly = shapely.geometry.Polygon(polygon)

line = [(4.0, -2.0000000000000004), (2.0, -1.1102230246251565e-15)]
shapely_line = shapely.geometry.LineString(line)

intersection_line = list(shapely_poly.intersection(shapely_line).coords)
print intersection_line

What I would expect is a list of two vertices.

Point 1: point that would be inside the polygon, or (4.0, -2.0000000000000004) in this case.

Point 2: point that is the intersection of [(4.0, -2.0000000000000004), (2.0, -1.1102230246251565e-15)] and [(3.0, -3.0), (4.0, -2.0)].

However, the result I receive is:

[(4.0, -2.0000000000000004)]

I have also checked whether there is an intersection at all with the edge that I'm looking at:

>>> edge = shapely.geometry.LineString([(3.0, -3.0), (4.0, -2.0)])
>>> edge.intersects(shapely_line)
False

If I replace (4.0, -2.0000000000000004) with (4.0, -2.000000000000000) then the edge intersection will evaluate to True.

Does anyone have any ideas for what is going on or what I am missing? Thanks!

enter image description here

Edit:

I have tested using shapely version 1.12 and geos of 3.3.1, 3.3.5, 3.3.6, 3.3.7.

In case anyone is curious as to how I updated the geos version on Windows:

Downloaded the geos-[version].tar.bz2 from the GEOS website. Extracted the files and ran CMake on it, using the Visual Studio 10 Win64 generator. Opened the .sln file and built it then moved the generated geos_c.dll and pasted it over where geos_c.dll had been installed by shapely in the Python directory.

Thanks unutbu for the picture, it looks better than the one I was about to draw! — lisaah
when I do intersection = shapely_poly.intersection(shapely_line) and then print intersection I get LINESTRING (4.0000000000000000 -2.0000000000000004, 4.0000000000000000 -2.0000000000000000) — flup
I found two properties that you can print to see the respective versions: shapely.__version__ and shapely.geos.geos_version_string. The latter ought to be of interest. — flup
it looks like a fun platform dependent robustness issue with GEOS and JTS (1.12, same answer on Win64) — Mike T

flup flup · Accepted Answer · 2013-01-22T00:11:30

Shapely is built on top of a C wrapper around the C++ GEOS library. Somewhere deep inside this C++ library sit the Precision classes which handle roundoff errors. I think we may conclude that your version of Shapely, and the geos libraries, handle this case differently. Unfortunately the code that accesses the precision model is not available in the C api, and therefore neither in Shapely. See http://lists.gispython.org/pipermail/community/2011-February/002898.html

Changing to a higher version of geos might solve your issue, though. It works fine on my machine with shapely 1.2.16 and libgeos 3.3.5-CAPI-1.7.5.

Polygon Intersection with Line | Python Shapely

1 Answers