Polygon intersection error in Shapely: "shapely.geos.TopologicalError: The operation 'GEOSIntersection_r' produced a null geometry"
I have been trying to debug this problem but unable to do so. I am trying to find the intersection of two Polygon objects. It works most of the time but for the following case, it raises the following exception:
P1 area: 13.125721955
P2 area: 1.0
Traceback (most recent call last):
File "geom2d.py", line 235, in <module>
print p1.intersection(p2)
File "/usr/local/lib/python2.7/dist-packages/shapely/geometry/base.py", line 334, in intersection
return geom_factory(self.impl['intersection'](self, other))
File "/usr/local/lib/python2.7/dist-packages/shapely/topology.py", line 47, in __call__
"The operation '%s' produced a null geometry. Likely cause is invalidity of the geometry %s" % (self.fn.__name__, repr(this)))
shapely.geos.TopologicalError: The operation 'GEOSIntersection_r' produced a null geometry. Likely cause is invalidity of the geometry <shapely.geometry.polygon.Polygon object at 0x8e5ad6c>
The code is as follows.
from shapely.geometry import Point,Polygon,MultiPolygon
poly1 = [(35.0041000000000011, -88.1954999999999956), (34.9917999999999978, -85.6068000000000069), (32.8404000000000025, -85.1756000000000029), (32.2593000000000032, -84.8927000000000049), (32.1535000000000011, -85.0341999999999985), (31.7946999999999989, -85.1358000000000033), (31.5199999999999996, -85.0438000000000045), (31.3384000000000000, -85.0836000000000041), (31.2092999999999989, -85.1069999999999993), (31.0023000000000017, -84.9943999999999988), (30.9953000000000003, -87.6008999999999958), (30.9422999999999995, -87.5926000000000045), (30.8538999999999994, -87.6256000000000057), (30.6744999999999983, -87.4072000000000031), (30.4404000000000003, -87.3687999999999931), (30.1463000000000001, -87.5240000000000009), (30.1545999999999985, -88.3863999999999947), (31.8938999999999986, -88.4742999999999995), (34.8937999999999988, -88.1020999999999930), (34.9478999999999971, -88.1721000000000004), (34.9106999999999985, -88.1461000000000041)]
poly2 = [(34.7998910000000024, -88.2202139999999986), (34.7998910000000024, -87.2202139999999986), (35.7998910000000024, -87.2202139999999986), (35.7998910000000024, -88.2202139999999986)]
p1 = Polygon(poly1)
p2 = Polygon(poly2)
print 'P1 area:',p1.area
print 'P2 area:',p2.area
print p1.intersection(p2)
Since it prints the areas of the two polygons, I assume that the polygons are formed correctly. I also (somehow) printed the first polygon to make sure that it is indeed a simple polygon.
Could anyone please explain why I am getting this exception?
Edit: I printed p1.is_valid and it turns out to be False. There is some explanation here. Search for the string is_valid. It says that
A valid Polygon may not possess any overlapping exterior or interior rings.
Could someone please explain what this means and if there is a possible work-around?
BTW, I also noticed that if I remove the last coordinate from poly1, the whole thing works. Perhaps the whole list of coordinates makes the polygon complex.
Answer
As previously pointed out, p1 is not valid. On plotting it, I noticed a little 'bowtie' at the lower right. I assume you don't need this in your polygon; if not, you can try Shapely's buffer(0) trick (documented in the Shapely Manual) to fix that:
In [382]: p1.is_valid
Out[382]: False
In [383]: p1 = p1.buffer(0)
In [384]: p1.is_valid
Out[384]: True
buffer(0) has the following effect:
Before:
After:
And you can now do this:
print p1.intersection(p2)
POLYGON ((34.9396324323625151 -88.1614025927056559, 34.8937999999999988 -88.1020999999999930, 34.7998910000000024 -88.1137513649788247, 34.7998910000000024 -87.2202139999999986, 34.9994660069532983 -87.2202139999999986, 35.0041000000000011 -88.1954999999999956, 34.9396324323625151 -88.1614025927056559))
Note that I have had some instances (with areas that looked more like "birds' nests" than simple bowties) where this didn't work; check to make sure you get back a single Polygon object and not a MultiPolygon one.