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

# # Overlays
# 
# Spatial overlays allow you to compare two GeoDataFrames containing polygon or multipolygon geometries 
# and create a new GeoDataFrame with the new geometries representing the spatial combination *and*
# merged properties. This allows you to answer questions like
# 
# > What are the demographics of the census tracts within 1000 ft of the highway?
# 
# The basic idea is demonstrated by the graphic below but keep in mind that overlays operate at the dataframe level, 
# not on individual geometries, and the properties from both are retained
# 
# ![illustration](http://docs.qgis.org/testing/en/_images/overlay_operations.png)

# Now we can load up two GeoDataFrames containing (multi)polygon geometries...

# In[ ]:


get_ipython().run_line_magic('matplotlib', 'inline')
from shapely.geometry import Point
from geopandas import GeoDataFrame, read_file
from geopandas.tools import overlay
from geodatasets import get_path

# NYC Boros
zippath = get_path("nybb")
polydf = read_file(zippath)

# Generate some circles
b = [int(x) for x in polydf.total_bounds]
N = 10
polydf2 = GeoDataFrame(
    [
        {"geometry": Point(x, y).buffer(10000), "value1": x + y, "value2": x - y}
        for x, y in zip(
            range(b[0], b[2], int((b[2] - b[0]) / N)),
            range(b[1], b[3], int((b[3] - b[1]) / N)),
        )
    ]
)


# The first dataframe contains multipolygons of the NYC boros

# In[ ]:


polydf.plot()


# And the second GeoDataFrame is a sequentially generated set of circles in the same geographic space. We'll plot these with a [different color palette](https://matplotlib.org/examples/color/colormaps_reference.html).

# In[ ]:


polydf2.plot(cmap="tab20b")


# The `geopandas.tools.overlay` function takes three arguments:
# 
# * df1
# * df2
# * how
# 
# Where `how` can be one of:
# 
#     ['intersection',
#     'union',
#     'identity',
#     'symmetric_difference',
#     'difference']
# 
# So let's identify the areas (and attributes) where both dataframes intersect using the `overlay` method. 

# In[ ]:


newdf = polydf.overlay(polydf2, how="intersection")
newdf.plot(cmap="tab20b")


# And take a look at the attributes; we see that the attributes from both of the original GeoDataFrames are retained. 

# In[ ]:


polydf.head()


# In[ ]:


polydf2.head()


# In[ ]:


newdf.head()


# Now let's look at the other `how` operations:

# In[ ]:


newdf = polydf.overlay(polydf2, how="union")
newdf.plot(cmap="tab20b")


# In[ ]:


newdf = polydf.overlay(polydf2, how="identity")
newdf.plot(cmap="tab20b")


# In[ ]:


newdf = polydf.overlay(polydf2, how="symmetric_difference")
newdf.plot(cmap="tab20b")


# In[ ]:


newdf = polydf.overlay(polydf2, how="difference")
newdf.plot(cmap="tab20b")