GeoPandas – Calculating Polygon Area within Another Polygon

areageopandasoverlapping-featurespolygonpython

I have a polygon GeoDataFrame 'Polygons' with attribute "AREA" (total polygon area in ha) and another GeoDataFrame 'Small_polygon' with smaller polygons within the 'Polygons'.

How do I calculate smaller polygon area within large polygons?

What I have already is Python code where I loop through each polygon in 'Polygons' and find corresponding smaller polygons within the 'Polygons' using the index.

After that I need to find each smaller polygon area within a large polygon.
For example, here is a polygon with Polygons['AREA'] = 100 ha and the result is calculated area (ha) for each small polygon.

I found some solutions in ArcGIS and QGIS, however I would like this to be done using Python.

enter image description here

I do not need to store the result anywhere, I will use them for further calculations within a loop, so this is kind of middle stage of the analyses. The further steps would be using the 'Some_points' layer calculate each smaller polygon average and max diameter, height and so on.

import pandas as pd
import geopandas as gpd

Polygons = gpd.read_file(r'E:\...\Polygons.shp')
Small_polygon = gpd.read_file(r'E:\...\Small_polygon.shp')
some_points = gpd.read_file(r'E:\...\points.shp')

print(Polygons.crs == Small_polygon.crs == some_points.crs)
True

# define the ID for each polygon
Polygons['ID'] = Polygons.index
id = Polygons['ID']
Polygons.drop(labels=['ID'], axis=1, inplace=True)
Polygons.insert(0, 'ID', id)

small_poly_join = gpd.sjoin(Small_polygon, Polygons)[['ID', 'small_id', 'geometry']]
points_join = gpd.sjoin(small_poly_join, some_points)[['ID',  'small_id', 'height', 'diameter']]

for i in Polygons.index:
    df_1 = small_poly_join[small_poly_join['FID'] == i]
    
    # this gives a wrong result
    area = (df_1['geometry'].area/ 10**6) # km2
    area = (df_1['geometry'].area/10000 # ha

    for n in df_1.index:
        df_2 = points_join[points_join['id'] == n]
        max_height =  df_2['height'].max()
        ave_height = df_2['height'].mean()
        
        max_diam = df_2['diameter'].max()
        ave_diam = df_2['diameter'].mean()

I had some problems with the area calculation that I found in other question, for instance, area = df_1['geometry'].area / 10000 gives me a wrong result compared to total polygon area. I also checked the CRS and projection which are equal to all used shapefiles.

Best Answer

Use GeoPandas Overlay

polygons = gpd.read_file("Polygons.shp")
small_polygon = gpd.read_file("Small_polygon.shp")

enter image description here

Intersection of the two GeoDataFrames:

result = gpd.overlay(polygons,small_polygon, how='intersection')

Result

enter image description here

Areas of the intersection polygons

result['area'] =result.apply(lambda row: row.geometry.area,axis=1)
Related Question