# Calculating polygon area within other polygon using GeoPandas

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

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.

``````polygons = gpd.read_file("Polygons.shp")
`````` Intersection of the two GeoDataFrames:

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

Result Areas of the intersection polygons

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

I don't use geopandas but in any GIS system what you are describing is a Union or Intersect style operation. This you can apply at the dataset level and I suspect would be far more efficient than looping over individual geometries as you are doing now. It appears that geopandas offers up such geoprocessing using overlay().

You're looking for a "Set Overlay" operation. Here's the way to do it in `geopandas`:

``````poly_intersections = Polygons.overlay(Small_polygon, how='intersection')
``````

Here's some more documentation on how to work with it: https://geopandas.org/en/stable/docs/user_guide/set_operations.html

Then you can get the area of the new cut/split polygons by doing this:

``````poly_intersections['area'] = poly_intersections.area
``````

Note: You might still need to find a conversion factor from whatever unit your original CRS was into hectares.

I think you make things a bit complicated, because there is a very useful module called `geopandas.overlay()`.

Let's assume there are two polygon layers 'grid' (comparable with your 'Small_polygon') and 'layer' (comparable with your 'Polygons'), see the image below. Using the following code:

``````import geopandas as gpd

_layer = "C:/Documents/Python Scripts/geopandas/layer.shp"
_grid = "C:/Documents/Python Scripts/geopandas/grid.shp"

_overlay = gpd.overlay(grid, layer, how='intersection')
_overlay['area'] = _overlay['geometry'].area/10**6

print(_overlay)
``````

it is possible to calculate smaller polygon area within large polygons. The output of `print(_overlay)` will look like:

``````        id  fid  ...                                           geometry         area
0     20.0  3.0  ...  POLYGON ((265049.059 5139668.500, 232793.533 5...  1070.267044
1     21.0  3.0  ...  POLYGON ((265049.059 5139668.500, 265049.059 5...  5212.682811
2     22.0  3.0  ...  POLYGON ((265049.059 5039668.500, 265049.059 4...  5912.757196
3     23.0  3.0  ...  POLYGON ((265049.059 4939668.500, 265049.059 4...  5775.394559
4     24.0  3.0  ...  POLYGON ((265049.059 4839668.500, 265049.059 4...  4081.788869
..     ...  ...  ...                                                ...          ...
121  165.0  2.0  ...  POLYGON ((1357832.462 5039668.500, 1265049.059...  2994.101400
122  174.0  2.0  ...  POLYGON ((1365049.059 5339668.500, 1365049.059...  2734.008176
123  175.0  2.0  ...  POLYGON ((1415809.234 5339668.500, 1417035.393...  5137.325491
124  176.0  2.0  ...  POLYGON ((1417035.393 5239668.500, 1418261.551...  5259.941295
125  177.0  2.0  ...  POLYGON ((1418261.551 5139668.500, 1418966.293...  4089.745958

[126 rows x 5 columns]
``````

How to check the way it works?

Let's consider the following area, where there is an overlap of four features at the same time. Three features from 'layer' intersect a feature with `"id" = 74` from 'grid' layer.

This feature (`"id" = 74`) has an area (Cartesian) of 10000.0 km².

When I filter the `_overlay` for features where `"id" = 74`, I will get the following output:

``````_overlay_filter = _overlay.loc[_overlay["id"] == 74]
print(_overlay_filter)

id  fid  ...                                           geometry         area
26  74.0  3.0  ...  POLYGON ((573320.195 5039668.500, 613697.864 4...  4378.454186
56  74.0  1.0  ...  POLYGON ((665049.059 5039668.500, 665049.059 5...  3205.802539
83  74.0  2.0  ...  POLYGON ((665049.059 4939668.500, 656836.249 4...  2415.743275

[3 rows x 5 columns]
``````

And to calculate the full sum of the overlapping features I can use this code:

``````_overlay_filter_sum = _overlay_filter.groupby(['id'])['area'].sum()
print(_overlay_filter_sum)

id
74.0    10000.0
Name: area, dtype: float64
``````

As you can see I got the initial 10000.0 km². However, you could expect the effect of rounding.