How to find geometries of all intersections in a polygon and calculate their areas and perimeters using shapely and fiona?

Question

I am able to open the shapefile and convert it into a shapely polygon but I am struggling with the logic required to isolate intersecting polygons and calculate their areas and perimeter.
from shapely.geometry import Polygon
import fiona

# open layer and create shapely polygon
polyShp = fiona.open('../data/areas.shp')
polyList = []
polyProperties = []

for poly in polyShp:
    polyGeom = Polygon(poly['geometry']['coordinates'][0])
    polyList.append(polyGeom)
    polyProperties.append(poly['properties'])

print(polyList[0])
print(polyProperties[0])

xunilk · Answer

You can try following script. It uses 'combinations' itertools method for avoiding unnecessary repetitions in feature intersections.
from shapely.geometry import Polygon, shape
import fiona
import itertools

# open layer and create shapely polygon
polyShp = fiona.open('/home/zeito/pyqgis_data/polygons2.shp')

geom_pol = [ shape(feat["geometry"]) for feat in polyShp ]

comb = range(len(geom_pol))

for i, j in itertools.combinations(comb, 2):
    if geom_pol[i].intersects(geom_pol[j]):
        inter = geom_pol[i].intersection(geom_pol[j])
        print("area[{},{}] = {} m2".format(i, j, inter.area))
        print("perimeter[{},{}] = {} meters".format(i, j, inter.length))

Running above script in Python Console produces following result:
area[0,1] = 538644.0179381497 m2
perimeter[0,1] = 3408.7918178154746 meters
area[1,2] = 1846868.1525187194 m2
perimeter[1,2] = 5971.712074534013 meters

Polygons look as follow in Map View of QGIS 3. Polygons layer uses transparency and id labeling for a better visualization of intersections areas and perimeters.

How to find geometries of all intersections in a polygon and calculate their areas and perimeters using shapely and fiona?

One Answer

Add your own answers!

Ask a Question