#!/usr/bin/env python
# coding: utf-8

# # Overlay polygons ("cascaded intersection + difference")
# 
# 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)

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')

import geopandas
import shapely.geometry
import shapely.ops

import itertools


# In[2]:


circle1 = shapely.geometry.Point(0, 0).buffer(1)
circle2 = shapely.geometry.Point(1, 1).buffer(1)
circle3 = shapely.geometry.Point(1, 0).buffer(1)


# In[3]:


s = geopandas.GeoSeries([circle1, circle2, circle3])


# In[4]:


s.plot(alpha=.5, cmap='Set1')


# Calculate all pairwise intersections:

# In[5]:


all_intersections = [a.intersection(b) for a, b in list(itertools.combinations(s, 2))]


# Combine those intersections with the original polygons:

# In[6]:


s_all = pd.concat([s, geopandas.GeoSeries(all_intersections)])


# Then, we take the boundaries of all those polygons:

# In[7]:


s_all.boundary.plot(cmap='Set1', alpha=.5)


# From those boundaries, we can polygonize the union of all linestrings:

# In[8]:


polys = list(shapely.ops.polygonize(s_all.boundary.unary_union))


# In[9]:


polys = geopandas.GeoSeries(polys)


# In[10]:


polys.plot(cmap='Set1')


# In[ ]: