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