5

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.

0
8

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)
5

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

4

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.

4

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.

input

Using the following code:

import geopandas as gpd

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

layer = gpd.read_file(_layer)
grid = gpd.read_file(_grid)

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

test

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.