CODEWITHSUNDEEP
pythoncountuniongeopandasoverlapping
jan_b
jan_b

Reputation: 123

Count overlapping features using Geopandas

Note: I've asked the same question here, but didn't get any feedback so far.

Is there a way to easily count overlapping polygons using Geopandas, the same way as the ArcGIS Pro Count Overlapping Features works?

enter image description here

So far my approach was to do union overlay and then dissolve with aggfunc='count' but for some reason the results I get are not correct.

I have 3 overlapping polygons in a single geodataframe:

enter image description here

enter image description here

Then I do the overlay:

union = gpd.overlay(gdf, gdf, how='union')

As a result I get only 9 polygons, although I should get 10 (this is what union in QGIS or ArcGIS would return):

enter image description here

Is there anything wrong with my approach? What is the best way to count overlapping polygons in a single geodataframe?

Full code is below. It returns 9 polygons. Based on my understanding on union/intersect operations, it should result in 10 polygons. The intersection of 3 polygons is counted only twice, not three times... The union operation in QGIS for the same set of polygons results in 10 polygons.

import pandas as pd
import matplotlib as plt
import geopandas as gpd
from shapely import wkt
data = {'name':  ['polygon A', 'polygon B', 'polygon C'],
        'id': [1, 2, 3],
         'geom': ['MULTIPOLYGON (((36.00000 11.00000, 36.00000 12.00000, 37.00000 12.00000, 37.00000 11.00000, 36.00000 11.00000)))', 'MULTIPOLYGON (((36.50000 11.50000, 37.50000 11.50000, 37.50000 11.00000, 36.50000 11.00000, 36.50000 11.50000)))', 'MULTIPOLYGON (((36.61799 10.80580, 36.61570 11.19321, 36.86327 11.29637, 37.34925 10.91813, 37.00540 10.71182, 36.61799 10.80580)))']
        }

df = pd.DataFrame (data, columns = ['name','id','geom'])
df['geom'] = df['geom'].apply(wkt.loads)
gdf = gpd.GeoDataFrame(df, geometry='geom')

gdf.plot(alpha=0.5, column='id')
union = gpd.overlay(gdf, gdf, how='union')
len(union)

Expected result:

enter image description here

Upvotes: 2

Views: 2147

Answers (1)

Matthew Borish
Matthew Borish

Reputation: 3096

This should work if you want the counts and output of the ArcGIS Pro Count Overlapping Features tool. If you want to include the original polygon, you can comment out/ remove this line: other_poly_list.remove(idx).

intersection_polygons_list = []

for idx, row in gdf.iterrows():

    main_poly_gdf = gdf.iloc[idx:idx+1, :]

    print('\n' + 'main poly:', main_poly_gdf['id'].tolist()[:])

    other_poly_list = gdf.index.tolist()

    other_poly_list.remove(idx)

    other_poly_gdf = gdf[gdf.index.isin(other_poly_list)]

    print('other polys:', other_poly_gdf['id'].values.tolist()[:])

    intersection_polygons = gpd.overlay(main_poly_gdf, other_poly_gdf, how='intersection')

    intersection_polygons_list.append(intersection_polygons)


intersections_gdf = pd.concat(intersection_polygons_list)

enter image description here

intersections_gdf_gb = intersections_gdf.groupby('name_1').agg({'name_2':'size'})

enter image description here

Upvotes: 1

Related Questions