Cómo cortar un shapefile de polígonos con Python, Fiona y Shapely - Tutorial

Este tutorial muestra el procedimiento completo para recortar una capa de polígonos a un área de interés en Python con el uso de bibliotecas espaciales como Fiona y Shapely. El tutorial abre el polígono y la capa de recorte como elementos de Fiona, interpreta las geometrías como polígonos de Shapely, recorta los polígonos y almacena los resultados como un archivo shapefile con los metadatos correspondientes.

Tutorial

Datos de entrada

Puede descargar los datos de entrada de este enlace.

Código

Este es el código completo en Python:

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()