Monday, 3 February 2020

shapefile - How to find the area with most intersection for buffer regions?



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)





No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...