How to find the intersection areas of overlapping buffer zones in single shapefile?

Question

Find the intersections of the buffer zones

I found an example where they tried finding the intersection of two different shapefiles -Intersect Shapefiles using Shapely. I've tried to adapt my code to do the same intersection, but I am only working with one shapefile with buffer features inside. The program is able to write the shapefiles, but the intersection aren't occurring, thus no data written to my output shapefile Anyone willing to help me understand where the problem lies or has an alternate way of finding intersection in a single shapefile?

# creation of the new shapefile with the intersection
for i in threshold:
    with fiona.open("threshold_buffer_shape_"+str(i)+".shp","r") as layer1:
        with fiona.open("threshold_buffer_shape_"+str(i)+".shp","r") as layer2:
            # We copy schema and add the  new property for the new resulting shp
            schema=layer2.schema.copy()
            schema['properties']['uGID'] = 'int:10'
            # create an empty spatial index object
            with fiona.open("intersection_region_"+str(i)+".shp","w","ESRI Shapefile",schema) as layer3:
                index = rtree.index.Index()
                # populate the spatial index
                for feat1 in layer1:
                    fid = int(feat1['properties']['GID'])
                    geom1 = shape(feat1['geometry'])
                    index.insert(fid,geom1.bounds)

            # get list of fids where bounding boxes intersect
                for feat2 in layer2:
                    geom2=shape(feat2['geometry'])

                    for fid in list(index.intersection(geom2.bounds)):
                        if fid!= int(feat2['properties']['GID']):
                            feat1 = layer1[fid]
                            geom1 = shape(feat1['geometry'])
                            if geom1.intersects(geom2):
                                props = feat2['properties']
                                props['uGID'] = feat1['properties']['uGID']
                                layer3.write({
                                    'properties':props,
                                    'geometry':mapping(geom1.intersection(geom2))
                                    })

Overlapping buffer regions

Answer 1

The problem here is that you need to iterate over two polygons of the layer of polygons as you do when you iterate through two layers (look at intersecting two shapefiles from Python or command line)

You have many solutions to iterate over every two polygons as explained in Shapely/ Python looping through a number of polygons (an you can even use rtree as you do)

import itertools
result = [poly1.intersection(poly2) for poly1,poly2 in  itertools.combinations(layer, 2) if poly1.intersects(poly2)]

But the problem is that the result depends on the order of the polygons in the layer and the resulting intersections are not intersected themselve.

The solution is to use the properties of unary_union with the LinearRings/LineString of the polygons (it cuts the lines at each intersection,Planar graph))

Then to use the polygonize function

rings = [LineString(list(pol.exterior.coords)) for pol in layer]
from shapely.ops import unary_union, polygonize
union = unary_union(rings)
result = [geom for geom in polygonize(union)]

And you can compute the areas of the resulting polygons.

New

I don't understand your problem. You say that you want to find the intersections in one shapefile (polygons looking at figure):

If I use a shapefile with 5 five polygons

With the proposed solution I don't need a spatial tree (rtree) nor the intersection of geometries (done by unary_union)

with fiona.open("polygons.shp") as layer:
   rings = [LineString(list(shape(pol['geometry']).exterior.coords)) for pol in layer]
   from shapely.ops import unary_union, polygonize
   union = unary_union(rings)
   result = [geom for geom in polygonize(union)]
   schema=layer.schema.copy()
   schema['properties']['area'] = 'float:10.2'
   with fiona.open('result.shp', 'w', 'ESRI Shapefile', schema) as c:
       for index, pol in enumerate(result):
           c.write({'geometry': mapping(pol),'properties': {'id': index, 'area': pol.area}})

Result

Control

with fiona.open('result.shp') as input:
    for pol in input:
      print pol['properties']['id'],pol['properties']['area']

0 1362.75
1 4905.55
2 2162.83
3 2818.66
4 5688.51
5 4536.22
6 98.17
7 1484.88
8 38.13
9 1271.71
10 4885.7
11 99.09
12 1594.85