Geopandas [Shapely] spatial difference: TopologyException: no outgoing dirEdge found

Question

I am working on the "solver" that sequentially selects best location from the list, basing on poi and overlapping polygons within the certain buffer, then removes those points and location from the loop, subtracts buffer for the selected location from buffers for all remaining locations, and starts again, looking for the next location.

In terms of tools, I use py2.7 and geopandas, which, in its turn, uses shapely -> geos.

It works fine, until at certain moment stops with the following error:

TopologyException: no outgoing dirEdge found at 37.584887352337432 55.775210354236364
TopologyException: found non-noded intersection between LINESTRING (37.8155 55.7494, 37.8155 55.7496) and LINESTRING (37.8156 55.7495, 37.8154 55.7496) at ...

Basing on suggestions, I tried to use .buffer(0) and even with a very small amount - it didn't help, and I am totally lost...

Answer 1

It is a GEOS well-known problem (with Shapely, RGeos-RGdal example or PotGIS example and ...). It simply means there is a line crossing another line and no intermediate coordinate records the intersection or by intersecting line segments which are nearly parallel

What are the solutions ?

1) rounding the numerical precision of the input geometries (limitations of using finite-precision numerical value in geometric algorithms)
2) create a buffer area around the polygons (not working for you)
3) use unary_union before the intersection (if they are MultiPolygons in your data)

Example with GeoPandas

import geopandas as gp
shape = gp.read_file("polys.shp")
shape

 for i in range(len(shape)):
     print i,  shape.iloc[i]['geometry']
 0 MULTIPOLYGON (((0 0, 0 10, 10 10, 10 0, 0 0)), ((0 0, 10 0, 10 -10, 0 -10, 0 0)))
 1 MULTIPOLYGON (((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)), ((0.5 0.5, 10.5 0.5, 10.5 -9.5, 0.5 -9.5, 0.5 0.5)))
 2 POLYGON ((-7.711842576354295 14.98571372222998, -3.183406557396946 13.24400756109254, -6.169188547918274 3.341163959196798, -11.14549186545382 10.40751467009727, -7.711842576354295 14.98571372222998))

Now the intersection

shape.iloc[0]['geometry'].intersection(shape.iloc[1]['geometry'])
# errors

But

from shapely.ops import unary_union
print unary_union(shape.iloc[0]['geometry']).intersection(unary_union(shape.iloc[1]['geometry'])).wkt
POLYGON ((10 0, 10 -9.5, 0.5 -9.5, 0.5 0.5, 0.5 10, 10 10, 10 0))

Inserting a new column in the GeoDataFrame

 new_shape = shape.copy()
 new_shape ['geometry']= shape.geometry.map(lambda x: unary_union(x))
 new_shape

 for i in range(len(new_shape)):
     print i,  new_shape.iloc[i]['geometry'].type
 0 Polygon
 1 Polygon
 2 Polygon
 # and
 print new_shape.iloc[0]['geometry'].intersection(new_shape.iloc[1])['geometry']).wkt
 POLYGON ((10 0, 10 -9.5, 0.5 -9.5, 0.5 0.5, 0.5 10, 10 10, 10 0))

New

How can I change numeric precision ?

Look at Is it possible to round all coordinates in shapely?

With GeoPandas we can do the same

Original GeoDataFrame (geoda)

 import numpy as np
 from shapely.geometry import shape, mapping
 def around(geom,p):
      geojson = mapping(geom)
      geojson['coordinates'] = np.round(np.array(geojson['coordinates']),p)
      return  shape(geojson)
 geoda.geometry= geoda.geometry.apply(lambda x: around(x,2))

Result