Intersect geometries in a geodataframe

Question

I have a geodataframe that contains some polygons (the buffer points). I am trying to find which of these polygons intersect (in the same gpd), but the option gpd.overlay(how=intersection) does not work. Although visibly the two buffers are intersected, as you can see in the image:

After finding those that intersect those elements assign the same centroid to them in order to identify them as the same element in the map (multipart) preserving all the elements.

I hope this reproducible example help to find an answer.

index	geometry
1	POINT (0.61328 41.61887)
2	POINT (0.61335 41.61897)
3	POINT (0.60776 41.61706)

x = gdf.geometry.buffer(0.000015)

With this buffer the points 1 and 2 intersect, while point 3 remains isolated. I need those buffers that intersect to have only one indicator and one coordinate.

Answer 1

• unary_union can be used to get merged geometries of overlapping buffered points. Then expand out multi-polygon to polygons
• add centroid to this geodataframe of merged points
• simple case of sjoin() to bring buffered points and merged polygons back together and have centroid to group overlapping buffered points (can also use index_right)
• have visualised with folium

import geopandas as gpd
import shapely

world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))

# get some points,  let's use Ukraine to show support for a country
# that has been invaded by a dictator with a little man complex
gdf = gpd.GeoDataFrame(
    world.loc[world["name"].eq("Ukraine"), "geometry"]
    .apply(lambda g: [shapely.geometry.Point(xy) for xy in g.exterior.coords])
    .explode(),
    crs=world.crs,
).reset_index(drop=True)

# get some buffered points some of which overlap
pb_test = gpd.GeoDataFrame(
    geometry=gdf.to_crs(gdf.estimate_utm_crs()).buffer(15000).to_crs(gdf.crs),
    crs=gdf.crs,
).reset_index(drop=True)

# merge overlapping polygons
merged = gpd.GeoDataFrame(
    geometry=[g for g in pb_test.unary_union.geoms], crs=pb_test.crs
).reset_index(drop=True)
# add centroid
merged["centroid"] = merged.centroid

# put centroid back on buffered points
pb_test = gpd.sjoin(pb_test, merged)

pb_test.explore(column="index_right", height=430, width=600)

Answer 2

If you are tying to find intersecting polygons within a single gdf, you can use sjoin on itself. Here's a toy example.

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
import matplotlib.pyplot as plt

df = pd.DataFrame()
df['x'] = [0.61328, 0.61335, 0.60776]
df['y'] = [41.61887, 41.61897, 41.61706]

gdf.crs= 4326

gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.x, df.y))

gdf['geometry'] = gdf['geometry'].buffer(0.0005)

buffer_sjoin_gdf = gpd.sjoin(gdf, gdf)
overlaying_buffers_gdf = buffer_sjoin_gdf.loc[buffer_sjoin_gdf.index != buffer_sjoin_gdf.index_right]

Output of overlaying_buffers_gdf:

x_left  y_left  geometry    index_right x_right y_right
1   0.61335 41.61897    POLYGON ((0.61385 41.61897, 0.61385 41.61892, ...   0   0.61328 41.61887
0   0.61328 41.61887    POLYGON ((0.61378 41.61887, 0.61378 41.61882, ...   1   0.61335 41.61897

Let's plot the original gdf:

f, ax=plt.subplots()
gdf.plot(color='red',ax=ax)

for idx, row in gdf.iterrows():
    ax.annotate(str(idx), (df['x'][idx], df['y'][idx]))

Now let's plot the overlaying_buffers_gdf:

f, ax=plt.subplots()
overlaying_buffers_gdf.plot(color='red',ax=ax)

for idx, row in overlaying_buffers_gdf.iterrows():
    ax.annotate(str(idx), (df['x'][idx], df['y'][idx]))

For the last part of your question, you can use an exploded unary union and then update the geometry with the centroid attribute.

poly_geom = overlaying_buffers_gdf.geometry.unary_union
uu_gdf_df = gpd.GeoDataFrame(geometry=[poly_geom])

uu_gdf_df_exp = uu_gdf_df.explode().reset_index(drop=True)
uu_gdf_df_exp['geometry'] = uu_gdf_df_exp['geometry'].centroid