Count overlapping features using Geopandas

Is there a way to easily count overlapping polygons using Geopandas, the same way as the ArcGIS Pro Count Overlapping Features works? 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:

Then I do the overlay:

union = gpd.overlay(demo_pg, demo_pg, 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): Is there anything wrong with my approach? What is the best way to count overlapping polygons in a single geodataframe?

EDIT: 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'])
gdf = gpd.GeoDataFrame(df, geometry='geom')

demo_pg.plot(alpha=0.5, column='id')
union = gpd.overlay(demo_pg, demo_pg, how='union')
len(union)
• Could you, please, provide your sample code so that others can use it to play around? Feb 19 '21 at 8:40
• @Stefan - I've included a full code in the main post Feb 20 '21 at 19:39

It took quite a bit of head scratching, but I finally got there! Here is how you can do it in GeoPandas.

def count_overlapping_features(in_gdf):
# Get the name of the column containing the geometries
geom_col = in_gdf.geometry.name

# Setting up a single piece that will be split later
input_parts = [in_gdf.unary_union.buffer(0)]

# Finding all the "cutting" boundaries. Note: if the input GDF has
# MultiPolygons, it will treat each of the geometry's parts as individual
# pieces.
cutting_boundaries = []
for i, row in in_gdf.iterrows():
this_row_geom = row[geom_col]
this_row_boundary = this_row_geom.boundary
if this_row_boundary.type[:len('multi')].lower() == 'multi':
cutting_boundaries = cutting_boundaries + list(this_row_boundary.geoms)
else:
cutting_boundaries.append(this_row_boundary)

# Split the big input geometry using each and every cutting boundary
for boundary in cutting_boundaries:
splitting_results = []
for j,part in enumerate(input_parts):
new_parts = list(shapely.ops.split(part, boundary).geoms)
splitting_results = splitting_results + new_parts
input_parts = splitting_results

# After generating all of the split pieces, create a new GeoDataFrame
new_gdf = gpd.GeoDataFrame({'id':range(len(splitting_results)),
geom_col:splitting_results,
},
crs=in_gdf.crs,
geometry=geom_col)

# Find the new centroids.
new_gdf['geom_centroid'] = new_gdf.centroid

# Starting the count at zero
new_gdf['count_intersections'] = 0

# For each of the `new_gdf`'s rows, find how many overlapping features
# there are from the input GDF.
for i,row in new_gdf.iterrows():
new_gdf.loc[i,'count_intersections'] = in_gdf.intersects(row['geom_centroid']).astype(int).sum()
pass

# Dropping the column containing the centroids
new_gdf = new_gdf.drop(columns=['geom_centroid'])[['id','count_intersections',geom_col]]

return new_gdf

This function takes in a GeoDataFrame and spits out a GeoDataFrame that contains all of the "unique" pieces of the original input. It also has a column called "count_intersections" that indicates how many of the features of the original input they intersect with.

When I apply the function to your input data, this is what I get:

# Applying the function to the input data
out_df = count_overlapping_features(gdf)
out_df.plot(alpha=0.5, column='count_intersections')

Caveat

This function works just fine with input GeoDataFrames that contain POLYGON features. However, if you have MULTIPOLYGON features, things might get weird. My code can handle a few corner cases, but there are definitely A LOT more that I didn't foresee. So investigate the results VERY CLOSELY if you're using this function on MULTIPOLYGON features.