Goal: given a set of overlapping polygons, create a new polygon layer with all polygons overlayed
Based on https://gis.stackexchange.com/a/326583/9828 (this is only all intersections, without the non-intersecting parts of the original polygons)
%matplotlib inline
import geopandas
import shapely.geometry
import shapely.ops
import itertools
circle1 = shapely.geometry.Point(0, 0).buffer(1)
circle2 = shapely.geometry.Point(1, 1).buffer(1)
circle3 = shapely.geometry.Point(1, 0).buffer(1)
s = geopandas.GeoSeries([circle1, circle2, circle3])
s.plot(alpha=.5, cmap='Set1')
<matplotlib.axes._subplots.AxesSubplot at 0x7fae967601d0>
Calculate all pairwise intersections:
all_intersections = [a.intersection(b) for a, b in list(itertools.combinations(s, 2))]
Combine those intersections with the original polygons:
s_all = pd.concat([s, geopandas.GeoSeries(all_intersections)])
Then, we take the boundaries of all those polygons:
s_all.boundary.plot(cmap='Set1', alpha=.5)
<matplotlib.axes._subplots.AxesSubplot at 0x7fae9467e7b8>
From those boundaries, we can polygonize the union of all linestrings:
polys = list(shapely.ops.polygonize(s_all.boundary.unary_union))
polys = geopandas.GeoSeries(polys)
polys.plot(cmap='Set1')
<matplotlib.axes._subplots.AxesSubplot at 0x7fae94605860>