¿Recorte de ráster mediante varios conjuntos de datos o polígonos?

Question

Me gustaría recortar un MDE utilizando una cuadrícula de polígonos. Probablemente es más fácil de usar múltiples polígonos en un archivo de forma, pero no he logrado esto, así que estoy tratando de usar un bucle for para que pueda bucle a través de cada conjunto de datos en un gdb (cada uno contiene sólo 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, sólo se queda ahí, esperando algo más... ¿pero qué? Puedo conseguir que funcione para un clip, pero no con el bucle.

Estoy seguro de que debería estar haciendo algo más para la salida, para nombrar a cada nuevo raster por clase de característica o algo así ... pero de nuevo, no sé cómo. Por favor, hágamelo saber si debo añadir más información.

Answer 1

Una cosa que noto es que su tercer parámetro es una salida codificada (C:/data/lidar). La forma en que está escrito ahora será un bucle a través de cada una de sus características y sobrescribir la salida cada vez, pero ya que no puede haber permitido la sobreescritura automática de archivos, esto podría ser potencialmente el hang-up. 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 la intención de colocar las salidas en la carpeta C:/data llamado lidar ... tenga en cuenta que el tercer parámetro en el clip es la ruta completa de la trama de salida, no una carpeta que se colocará en. Si no especifica una extensión en su nombre de ruta de salida y los coloca en una carpeta estándar, será una cuadrícula, por lo que ahora mismo su programa está intentando crear un nuevo conjunto de datos de cuadrícula llamado 'lidar' en la carpeta C:/data.

Answer 2

Aquí tienes dos opciones:

1. Utilice la herramienta integrada de ArcGIS denominada Split Raster (Datos de datos) .

   Divide un conjunto de datos ráster en piezas más pequeñas, por mosaicos o características de un polígono.

2. Pruebe la herramienta Raster Split Tool, disponible en el USGS (véase adjunto código fuente adjunto).

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)

Answer 3

Para futuros buscadores: He aquí una versión modificada de la Herramienta de división de rásteres del 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 (jpoling@anra.org)
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')