Cortando a varredura por vários conjuntos de dados ou polígonos?

8

Gostaria de recortar um DEM usando uma grade de polígonos. Provavelmente é mais fácil usar vários polígonos em um arquivo de forma, mas eu não consegui fazer isso, então estou tentando usar um loop for para poder percorrer cada conjunto de dados em um gdb (cada um contém apenas um polígono).

Aqui está o meu código (fazendo isso na janela 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")

Meu código não é executado no entanto, apenas fica lá, esperando por outra coisa ... mas o que? Posso fazê-lo funcionar para um clipe, mas não com o loop.

Tenho certeza de que deveria estar fazendo outra coisa pela saída, para nomear cada nova varredura por classe de recurso ou algo assim ... mas, novamente, não sei como. Informe-me se devo adicionar mais informações.

Rosie Bell
fonte
2
Talvez teste primeiro se o seu loop está funcionando, comentando o clipe e colocando uma impressão ou AddMessage para ver se ele cospe o nome de cada classe de recurso em OK.
PolyGeo
Você poderia explicar por que precisa fazer esse recorte? A menos que a saída seja necessária para a comunicação com outro software, geralmente existem métodos muito mais eficientes para analisar dados raster polígono por polígono; portanto, talvez a melhor resposta para sua pergunta seja sugerir uma abordagem completamente diferente.
whuber
Desculpas por publicar e executar - agradeço a todos por sua ajuda e voltarei a acompanhar isso o mais rápido possível.
Rosie Bell
wuber, eu queria fazer isso para poder dividir nossos DEMs em pedaços gerenciáveis ​​para criar contornos e costurar novamente. Provavelmente não é a maneira mais eficiente de fazê-lo.
Rosie Bell
1
Qual software está usando? ArcGIS, QGIS, etc?
Chad Cooper

Respostas:

6

Uma coisa que noto é que seu terceiro parâmetro é uma saída codificada (C: / data / lidar). A maneira como ele é gravado agora percorre cada um dos seus recursos e sobrescreve a saída toda vez, mas como você pode não ter permitido a substituição automática de arquivos, isso pode ser um problema. Tente criar um nome de saída exclusivo para cada iteração:

#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")

Além disso, não tenho certeza de que você pretendia colocar as saídas na pasta C: / data chamada lidar ... observe que o terceiro parâmetro no clipe é o caminho completo da sua varredura de saída, não uma pasta na qual ela será colocada Se você não especificar uma extensão no nome do caminho de saída e colocá-las em uma pasta padrão, ela será uma grade; portanto, agora o seu programa está tentando criar um novo conjunto de dados da grade chamado 'lidar' no C: / pasta de dados.

pé azul
fonte
Muito obrigado mbenedetti, foi a falta de um nome único para cada novo clipe que foi um dos meus problemas. Eu também tinha algo estranho acontecendo com o fato de ser um gdb - quando tentei isso usando shapefiles em uma pasta, funcionou bem. Então eu provavelmente não tinha o caminho correto do arquivo (não estava acostumado a trabalhar com o gdb).
Rosie Sino
5

para futuros candidatos: aqui está uma versão modificada do script da ferramenta de rasterização USGS que não requer nada acima do nível de licença do 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')
Jeremiah Poling
fonte
4

Algumas idéias:

  1. Experimente a Raster Split Tool , disponível como ferramenta de script no USGS (consulte o código-fonte em anexo).
  2. Para recortar / ladrilhar raster simples, use a ferramenta ArcGIS 10.1 integrada chamada Split Raster. Você pode dividir as rasters pelo número de peças ou pelo tamanho da peça.

insira a descrição da imagem aqui

Código fonte da ferramenta 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)
Aaron
fonte
Existe um limite de tamanho de arquivo para a Raster Split Tool? Minha varredura é de 20,8 GB! Eu continuo recebendo "falha ao executar". Todos os arquivos têm a mesma projeção. Rodando o ArcMap 10.1
Derelict
Você sempre pode executá-lo em modo de fundo geoprocessamento 64 bit: resources.arcgis.com/en/help/main/10.1/index.html#//...
Aaron