How to clip polygon layers with Python, Fiona and Shapely - Tutorial

This tutorial shows the entire procedure to clip a polygon layer to an area of interest in Python with the use of spatial libraries as Fiona and Shapely. The tutorial opens the polygon and clip layer as fiona elements, interpret the geometries as shapely Polygon datatypes, clip the polygons and store results as an output shapefile with the corresponding metadata.

Tutorial

Input data

You can download the input data from this link.

Code

This is the entire Python code used for the tutorial:

import fiona
from shapely.geometry import Polygon, mapping
import matplotlib.pyplot as plt
from descartes import PolygonPatch

#open clip layer
aoi = fiona.open('../shps/AOI.shp')
aoiGeom = Polygon(aoi[0]['geometry']['coordinates'][0])

#open polygon layer and create list of Polygons
polyShp = fiona.open('../shps/Collier_AG_2017.shp')
polyList = []
polyProperties = []
for poly in polyShp:
    polyGeom = Polygon(poly['geometry']['coordinates'][0])
    polyList.append(polyGeom)
    polyProperties.append(poly['properties'])
print(polyList[0])
print(polyProperties[0])

POLYGON ((437765.2307000002 2881459.169199999, 437752.0489999996 2881446.275599999, 437424.8422999997 2881449.459100001, 437425.9841999998 2881481.255100001, 437454.1131999996 2882264.426999999, 437464.8493999997 2882600.848200001, 437483.1067000004 2882600.229699999, 437859.3836000003 2882582.1403, 437822.2763999999 2881571.759500001, 437810.8629000001 2881571.657600001, 437804.5048000002 2881480.533500001, 437804.5022999998 2881480.497300001, 437765.2307000002 2881459.169199999))
OrderedDict([('OBJECTID', 6622), ('Acres', 110.607924009), ('County', 'Collier'), ('FLUCCS', 2140), ('FLUCCSDESC', 'Vegetables - Row crops'), ('FLUCCS_L1', 2000), ('L1_DESC', 'AGRICULTURE'), ('FLUCCS_L2', 2100), ('L2_DESC', 'Cropland and Pastureland'), ('FLUCCS_L3', 2140), ('L3_DESC', 'Vegetables - Row crops'), ('FLUCCS_L4', 0), ('L4_DESC', None), ('Source', 'Balmoral Group'), ('Descript', 'Vegetables'), ('Spring_Veg', None), ('Fall_Veg', None), ('CropSeason', None), ('Irrigated', 'Yes'), ('Irr_Type', 'Drip'), ('WaterSourc', 'GW'), ('DateVerify', '2016-11-15'), ('Date_Fall', None), ('Remarks', None), ('Crop_Type', 'Row Crops - Vegetables'), ('Irr_System', 'Micro'), ('VerifyStat', None), ('WMD', 'SFWMD'), ('WaterSrc2', None)])

#plot sample of polygon and clip layer
fig, ax = plt.subplots(figsize=(12,8))

for poly in polyList:
    polyPatch = PolygonPatch(poly,alpha=0.5,color='orangered')
    ax.add_patch(polyPatch)
clipPatch = PolygonPatch(aoiGeom,alpha=0.8,color='royalblue')
ax.add_patch(clipPatch)
outline = 2000
ax.set_xlim(aoiGeom.bounds[0]-outline,aoiGeom.bounds[2]+outline)
ax.set_ylim(aoiGeom.bounds[1]-outline,aoiGeom.bounds[3]+outline)
plt.show()

#clip polygons to clip layer
clipPolyList = []
clipPolyProperties = []
for index, poly in enumerate(polyList):
    result = aoiGeom.intersection(poly)
    if result.area:
        clipPolyList.append(result)
        clipPolyProperties.append(polyProperties[index])
print(clipPolyList[0])
print(clipPolyProperties[0])

POLYGON ((452015.5206000004 2910386.9307, 451943.4771999996 2910207.359200001, 451922.5979000004 2910095.9868, 451825.9227999998 2910091.8718, 451631.2663000003 2910096.2894, 451625.6064999998 2910206.688899999, 451653.9197000004 2910271.8037, 451695.6659000004 2910370.820900001, 451697.4475999996 2910372.834100001, 451705.8526999997 2910382.373400001, 451740.7489 2910421.917400001, 451764.7916000001 2910445.962200001, 451861.8382999999 2910450.1132, 451862.4753 2910447.4628, 452048.2598999999 2910439.5132, 452015.5206000004 2910386.9307))
OrderedDict([('OBJECTID', 6755), ('Acres', 26.7604326734), ('County', 'Collier'), ('FLUCCS', 2141), ('FLUCCSDESC', 'Vegetables - Tomatoes'), ('FLUCCS_L1', 2000), ('L1_DESC', 'AGRICULTURE'), ('FLUCCS_L2', 2100), ('L2_DESC', 'Cropland and Pastureland'), ('FLUCCS_L3', 2140), ('L3_DESC', 'Vegetables - Row crops'), ('FLUCCS_L4', 2141), ('L4_DESC', 'Tomatoes'), ('Source', 'Balmoral Group; SFWMD'), ('Descript', 'Vegetables - Tomatoes'), ('Spring_Veg', None), ('Fall_Veg', None), ('CropSeason', None), ('Irrigated', 'Yes'), ('Irr_Type', 'Subsurface Seepage'), ('WaterSourc', 'GW'), ('DateVerify', '2016-11-15'), ('Date_Fall', None), ('Remarks', 'Permit number 11-00112-W'), ('Crop_Type', 'Row Crops - Vegetables'), ('Irr_System', 'Flood'), ('VerifyStat', None), ('WMD', 'SFWMD'), ('WaterSrc2', None)])

#plot clipped polygon and clip layer
fig, ax = plt.subplots(figsize=(12,8))

for poly in clipPolyList:
    polyPatch = PolygonPatch(poly,alpha=0.5,color='orangered')
    ax.add_patch(polyPatch)
clipPatch = PolygonPatch(aoiGeom,alpha=0.5,color='royalblue')
ax.add_patch(clipPatch)
outline = 2000
ax.set_xlim(aoiGeom.bounds[0]-outline,aoiGeom.bounds[2]+outline)
ax.set_ylim(aoiGeom.bounds[1]-outline,aoiGeom.bounds[3]+outline)
plt.show()

#export clipped polygons as shapefile
schema = polyShp.schema

outFile = fiona.open('../shps/clippedPolygons.shp',mode = 'w',driver = 'ESRI Shapefile', schema=schema)
for index, poly in enumerate(clipPolyList):
    outFile.write({
        'geometry':mapping(poly),
        'properties':clipPolyProperties[index]
    })
outFile.close()