3

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))
                                    })

enter image description here

Overlapping buffer regions

6

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)]

enter image description hereenter image description hereenter image description hereenter image description here

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))

enter image description here

Then to use the polygonize function

enter image description here

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

enter image description here

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

enter image description here

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
  • I managed to integrate your code in my solution. It gives separate polygons (circle parts) for newly created shapefile layer. I'm now trying to find the intersection of my new shapefile (parts) with original shapefile to get the part which overlaps/intersects the most with original (hence the middle/centre which has intersection of 3). I've tried intersect, but it gives shapefile similiar to the newly created shapefile layer. This is due to the fact that it says that far left part intersects with far right polygon (original layer). @gene – Jamesy Apr 4 '16 at 13:56
  • look above ---- – gene Apr 4 '16 at 16:18
  • My problem is a new question where I am trying to find the part which the polygons in my layer intersect the most. So for example: If poly in result intersect/(falls within) a polygon in original layer, then count and store. The poly with most count will be the most intersected part right? That's what I want to achieve. So far it wont give me the most intersected part. It says that the polygons all intersect each other. If created a new question here- [link] (gis.stackexchange.com/questions/187998/…) @gene – Jamesy Apr 5 '16 at 6:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.