Union the polygon boundaries into one multiline. Polygonize it and intersect with original df. This will create duplicate geometries where there are overlaps. Drop duplicates on centroid coordinates and area, keeping last (=topmost).
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.ops import unary_union, polygonize
#Read a shapefile into a dataframe and plot
file = r"C:\Users\bera\Desktop\gistest\Overlapping_squares_1k.shp"
df = gpd.read_file(file)
df["id"] = range(df.shape[0]) #Create an id column
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(10,10))
df.plot(ax=ax[0][0], alpha=0.3, color="red")
ax[0][0].set_title("Polygons with overlaps")
#Union the polygon boundaries into one multiline
multiline = unary_union(df.boundary)
gpd.GeoSeries(data=multiline, crs=df.crs).plot(ax=ax[0][1])
ax[0][1].set_title("The unioned boundaries")
#Polygonize the unioned boundaries
poly_df = gpd.GeoDataFrame(geometry=list(polygonize(multiline)), crs=df.crs)
poly_df.plot(ax=ax[1][0], color="red", alpha=0.3)
ax[1][0].set_title("The polygonized boundaries")
#Intersect the polygonized multiline boundary with the original df
inter = gpd.overlay(df1=df, df2=poly_df, how="intersection",
keep_geom_type=True, make_valid=True)
inter.plot(ax=ax[1][1], alpha=0.3, color="green")
ax[1][1].set_title("Original df and polygonized df inters.")
#Create columns of rounded centroid coordinates and area
# to use for detecting duplicates
inter["cent_x"] = round(inter.centroid.x)
inter["cent_y"] = round(inter.centroid.y)
inter["area_m2"] = inter.area.round()
#Sort ascending so larger ids (polygons on top) come last
inter = inter.sort_values(by="id")
#Drop duplicates on centroid coords and area, keep last/topmost
inter = inter.drop_duplicates(subset=["cent_x", "cent_y", "area_m2"], keep="last")
inter = inter.dissolve(by="id") #Dissolve split polygons back to original shapes
inter.to_file(r"C:\Users\bera\Desktop\gistest\confusion\no_overlaps.gpkg")
