Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Polygon Intersection with Line | Python Shapely

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.

like image 852
lisaah Avatar asked Jan 21 '13 22:01

lisaah


1 Answers

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.

like image 128
flup Avatar answered Sep 22 '22 11:09

flup