Polygon Intersection with Line | Python Shapely

python line polygon intersection shapely
lisaah picture lisaah · Jan 21, 2013 · Viewed 14.5k times · Source

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.

Answer

flup picture flup · Jan 22, 2013

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.

Related questions

Reading specific lines only

I'm using a for loop to read a file, but I only want to read specific lines, say line #26 and #30. Is there any built-in feature to achieve this? Thanks

Read More
How can I format a list to print each element on a separate line in python?

How can I format a list to print each element on a separate line? For example I have: mylist = ['10', '12', '14'] and I wish to format the list so it prints like this: 10 12 14 so the \n, …

Read More
Search and get a line in Python

Is there a way to search, from a string, a line containing another string and retrieve the entire line? For example: string = qwertyuiop asdfghjkl zxcvbnm token qwerty asdfghjklñ retrieve_line("token") = "token qwerty"

Read More