¿Ráster de recorte por múltiples conjuntos de datos o polígonos?

Me gustaría recortar un DEM usando una cuadrícula de polígonos. Probablemente sea más fácil usar múltiples polígonos en un archivo de forma, pero no lo he logrado, así que estoy tratando de usar un bucle for para poder recorrer cada conjunto de datos en un gdb (cada uno contiene solo un polígono).

Aquí está mi código (haciéndolo en la ventana de Python).

#creating a workspace and a list of feature classes
arcpy.env.workspace = "C:/data/lidar/lidar.gdb"
fcs = arcpy.ListFeatureClasses()

#looping through each feature class and creating a raster based on the extent of   
#feature class

for fc in fcs:
    arcpy.Clip_management("perth", "#", "C:/data/lidar", fc, "", "ClippingGeometry")

Sin embargo, mi código no se ejecuta, simplemente se queda allí, esperando algo más ... pero ¿qué? Puedo hacer que funcione para un clip, pero no con el bucle.

Estoy seguro de que debería hacer algo más para la salida, para nombrar cada nuevo ráster por clase de entidad o algo así ... pero de nuevo, no sé cómo. Avíseme si debo agregar más información.

Una cosa que noté es que su tercer parámetro es una salida codificada (C: / data / lidar). La forma en que está escrito ahora recorrerá cada una de sus características y sobrescribirá la salida cada vez, pero dado que es posible que no haya permitido la sobrescritura automática de archivos, esto podría ser el problema. Intente crear un nombre de salida único para cada iteración:

#creating a workspace and a list of feature classes
arcpy.env.workspace = "C:/data/lidar/lidar.gdb"
fcs = arcpy.ListFeatureClasses()

#looping through each feature class and creating a raster based on the extent of   
#feature class

i=0
for fc in fcs:
    outputPath = "C:/data/lidar" + str(i)
    i+=1
    arcpy.Clip_management("perth", "#", outputPath, fc, "", "ClippingGeometry")

Además, no estoy seguro de que tuviera la intención de colocar las salidas en la carpeta C: / data llamada lidar ... tenga en cuenta que el tercer parámetro en el clip es la ruta completa de su ráster de salida, no una carpeta que se colocará in. Si no especifica una extensión en el nombre de la ruta de salida y la coloca en una carpeta estándar, será una cuadrícula, por lo que ahora su programa está intentando crear un nuevo conjunto de datos de cuadrícula llamado 'lidar' en C: / carpeta de datos.

para futuros buscadores: Aquí hay una versión modificada del script de la herramienta de división de ráster USGS que no requiere nada por encima del nivel de licencia de ArcGIS Basic (ArcView):

"""
Raster Split Tool 6/16/2011
ArcGIS 10 Script Tool
Python 2.6.5

Contact:
Douglas A. Olsen
Geographer
Upper Midwest Environmental Sciences Center
U.S. Geological Survey
2630 Fanta Reed Road
La Crosse, Wisconsin 54603
Phone: 608.781.6333

#####
modified 2015-04-15 by Jeremiah Poling ([email protected])
Now using only tools available at the ArcGIS Basic license level 
#####

"""
import arcpy
from arcpy import env
import os

# Get Split Shapfile Name
splitShape = arcpy.GetParameterAsText(0)
# Get Field Name
splitField = arcpy.GetParameterAsText(1)
# Get Output Directory
outDirectory = arcpy.GetParameterAsText(2)
# Get Raster to Split
splitRaster = arcpy.GetParameterAsText(3)
# Get type of output tif, img, or GRID
rasterType = arcpy.GetParameterAsText(4)
# Set Workspace environment
env.workspace = outDirectory

# Loop through the rows in the clipping shapefile
cursor = arcpy.SearchCursor(splitShape)
for row in cursor:
    resultName = row.getValue(splitField)

    # Create feature layer of current clipping polygon
    whereClause = '"' + splitField + '" = ' + "'" + resultName + "'"
    arcpy.MakeFeatureLayer_management(splitShape, 'currentMask', whereClause)

    # Replace spaces with underscores
    resultName = resultName.replace(' ','_')

    # Remove .shp suffix
    if resultName[-4:] == '.shp':
        resultName = resultName[:-4]

    arcpy.AddMessage("Processing: " + resultName)

    if rasterType == 'img':
        resultName = resultName + "." + rasterType
    elif rasterType == 'tif':
        resultName = resultName + "." + rasterType
    else:
        print ('No extension')

    # Save the clipped raster
    arcpy.Clip_management(
        in_raster = splitRaster,
        rectangle = "#",
        out_raster = resultName,
        in_template_dataset = 'currentMask',
        nodata_value="255",
        clipping_geometry="ClippingGeometry",
        maintain_clipping_extent="NO_MAINTAIN_EXTENT"
        )

    if arcpy.Exists('currentMask'):
        arcpy.Delete_management('currentMask')

Algunas ideas

1. Pruebe la herramienta de división de ráster , disponible como herramienta de secuencia de comandos de USGS (consulte el código fuente adjunto).
2. Para el recorte / mosaico de ráster simple, use la herramienta integrada ArcGIS 10.1 llamada Split Raster. Puede dividir los rásteres por la cantidad de mosaicos o por el tamaño de los mosaicos.

Código fuente de la herramienta USGS Raster Split:

"""
Raster Split Tool 6/16/2011
ArcGIS 10 Script Tool
Python 2.6.5

Contact:
Douglas A. Olsen
Geographer
Upper Midwest Environmental Sciences Center
U.S. Geological Survey
2630 Fanta Reed Road
La Crosse, Wisconsin 54603
Phone: 608.781.6333
"""
import arcpy
from arcpy import env
import os
from arcpy.sa import *

# Get Shapefile Name
inShape = arcpy.GetParameterAsText(0)
# Get Split Shapfile Name
splitShape = arcpy.GetParameterAsText(1)
# Get Field Name
splitField = arcpy.GetParameterAsText(2)
# Get Output Directory
outDirectory = arcpy.GetParameterAsText(3)
# Get Raster to Split
splitRaster = arcpy.GetParameterAsText(4)
# Get type of output tif, img, or GRID
rasterType = arcpy.GetParameterAsText(5)
# Set Workspace environment
env.workspace = outDirectory

# Run Split command on Input Shapefile
arcpy.Split_analysis(inShape, splitShape, splitField, outDirectory)

# Loop through a list of feature classes in the workspace
for fc in arcpy.ListFeatureClasses():

    # Execute ExtractByMask
    rfc = ExtractByMask(splitRaster, fc)

    # Clean up shp file from work directory
    arcpy.Delete_management(fc, "")    

    arcpy.AddMessage("Processing: " + fc)

    # Replace spaces with underscores
    fc = fc.replace(' ','_')

    # Remove .shp suffix
    fc = fc[:-4]

    # Trim to 13 chars
    fc = fc[0:13]

    if rasterType == 'img':
        fc = fc + "." + rasterType
    elif rasterType == 'tif':
        fc = fc + "." + rasterType
    else:
        print ('No extension')

    # Save the output
    rfc.save(fc)