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
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
No comments:
Post a Comment