[GIS] How to simplify borders between polygons in shapely

geoprocessingpythonshapely

I have some fairly complex, and not completely accurate shapes stored as shapely Polygon objects. I'm trying to find the common borders, which should be as easy as finding the LineString returned from poly1.intersection(poly2).

However, due the the slight inaccuracies I'm finding that intersection is returning Polygon and Point objects, or sometimes a GeometryCollection object. Never a simple LineString except in my simple test code where the Polygon shapes are accurate.

Is there any way of simplifying the borders, like a "snap-to" method, or a tolerance parameter that might ease my pain?

Best Answer

Since version 1.5, shapely provides a snap function that can be used in that way.

From the shapely documentation

from shapely.ops import snap
square = Polygon([(1,1), (2, 1), (2, 2), (1, 2), (1, 1)])
line = LineString([(0,0), (0.8, 0.8), (1.8, 0.95), (2.6, 0.5)])
result = snap(line, square, 0.5)

Related Solutions

Shapely – Translating Between Polygons and MultiPolygons Using Shapely

Further to relet's answer on how to get individual polygons, you can then run an intersection on all the polygons to create the holes. If your dataset contains overlapping polygons though you're out of luck.

Explain again what is wrong with existing shapefile readers?

Would it not be easier to export feature IDs and M values from the shapefile and then join them back to the polygons after using an existing shapefile reader?

For multipatches you can use the same technique of assigning polygon IDs to a "patch ID" and then adding this attribute back to the features.

Edit: Whilst you say you don't want to use OGR, just in case you change your mind..

import ogr
# Get the driver
driver = ogr.GetDriverByName('ESRI Shapefile')
# Open a shapefile
shapefileName = "D:/temp/myshapefile.shp"
dataset = driver.Open(shapefileName, 0)

layer = dataset.GetLayer()
for index in xrange(layer.GetFeatureCount()):
    feature = layer.GetFeature(index)
    geometry = feature.GetGeometryRef()
    #geometry for polygon as WKT, inner rings, outer rings etc. 
    print geometry

The geometry should be output as follows:

POLYGON ((79285 57742,78741 54273...),(76087 55694,78511 55088,..))

The first bracket contains the coords of the exterior ring, subsequent brackets the coords of interior rings. If you have Z values points should be in the format 79285 57742 10 (where the last coord is a height).

Otherwise you could use the Shapely Contains and Within functions to assess every polygon with each other and apply a spatial index beforehand - http://pypi.python.org/pypi/Rtree/ to speed up processing.

[GIS] Python shapely difference between two polygons

Once you've calculated the difference between your polygons (let's call this delta), take the difference between the exterior of delta and the exterior of your brown polygon2: result = delta.exterior.difference(polygon2.exterior). Note that the sense of the difference will be significant and it may be tricky to generalize this to all possible polygons. Good luck!

Best Answer

Related Solutions

Shapely – Translating Between Polygons and MultiPolygons Using Shapely

[GIS] Python shapely difference between two polygons

Related Question