#!/usr/bin/env python # coding: utf-8 # 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 # In[1]: from IPython.core.display import Image Image(url="http://docs.qgis.org/testing/en/_images/overlay_operations.png") # Now we can load up two GeoDataFrames containing (multi)polygon geometries... # In[2]: get_ipython().run_line_magic('matplotlib', 'inline') from shapely.geometry import Point from geopandas import datasets, GeoDataFrame, read_file from geopandas.tools import overlay # NYC Boros zippath = datasets.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[3]: 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[4]: 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` tool. # In[5]: from geopandas.tools import overlay newdf = overlay(polydf, 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[6]: polydf.head() # In[7]: polydf2.head() # In[8]: newdf.head() # Now let's look at the other `how` operations: # In[9]: newdf = overlay(polydf, polydf2, how="union") newdf.plot(cmap='tab20b') # In[10]: newdf = overlay(polydf, polydf2, how="identity") newdf.plot(cmap='tab20b') # In[11]: newdf = overlay(polydf, polydf2, how="symmetric_difference") newdf.plot(cmap='tab20b') # In[12]: newdf = overlay(polydf, polydf2, how="difference") newdf.plot(cmap='tab20b')