1

Find the area with the most intersection

This question relates to previous question asked - How to find the intersection areas of overlapping buffer zones in single shapefile?. I tried to find the region which intersects the most. I managed to break my original shapefile into polygon parts and now trying to find the part which intersects the the most with the original to output area of interest in the center. Anyone willing to help me understand where the problem lies or has an alternate way of finding the region of interest?

# creation of the new shapefile with the intersection
for i in threshold:
    with fiona.open("threshold_buffer_test_"+str(i)+".shp","r") as layer1:
        with fiona.open("threshold_buffer_test_"+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']['uid'] = '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['id'])
                    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['id']):
                            feat1 = layer1[fid]
                            geom1 = shape(feat1['geometry'])
                            if geom1.intersects(geom2):
                                print 'hello world'
                                props = feat2['properties']
                                props['uid'] = feat1['properties']['uid']
                                layer3.write({
                                    'properties':props,
                                    'geometry':mapping(geom1.intersection(geom2))
                                    })

enter image description here

Original buffer shapefile

enter image description here

Polygonized shapefile (parts)

Area of interest in center

Area of interest in center (intersects the most with the original polygons)

1 Answer 1

1

I still think that I not understand what you are trying to do with your script (and his complexity)

The original layer

enter image description here

I use the simple script presented in How to find the intersection areas of overlapping buffer zones in single shapefile? (new) (without rtree, many for loops and intersections) to create the resulting shapefile.

enter image description here

By definition, all the resulting polygons intersect the 3 original polygons and vice versa (and none of them intersects the most with the original polygons)

# the original shapefile: 3 blue polygons
ori = [shape(pol['geometry']) for pol in fiona.open('original.shp')] # blue
# the resulting shapefile: 7 polygons
res = [shape(pol['geometry']) for pol in fiona.open('result.shp')] # blue
for i in range(len(res)):
   for j in range(len(ori)):
        print i,j,ori[i].intersects(res[j]

0 0 True
0 1 True
0 2 True
1 0 True
1 1 True
1 2 True
2 0 True
2 1 True
.... (all true)

And the resulting polygon with maximum area is 1 in mi case

Therefore, in this case, yellow polygon (5) has nothing in particular and it is not the region which intersects the most. It is effectively in the centre. This is what you want to find, the central position ?

You can control all the predicates with the Java Topology Suite (JTS) (here original polygon 3 and resulting polygon 6)

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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