Friday 19 August 2016

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


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



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


enter image description here


 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

enter image description here


 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)


enter image description here



 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


enter image description here



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...