Uso de la herramienta Longest Flow Path de HEC-GeoHMS en bucle con Arcpy

Question

Tengo una clase de características formada por cientos de polígonos, para la que necesito determinar la ruta de flujo más larga utilizando HEC-GeoHMS. Como cualquiera que haya utilizado esta herramienta sabe, no se puede pasar al conjunto de herramientas toda la clase de características si alguna de las cuencas se solapa - no calculará correctamente el camino de flujo más largo para las cuencas que se solapan.

Por lo tanto, estoy tratando de configurar un script que itera a través de cada registro en mi clase de características y genera un resultado de la ruta de flujo más larga para cada registro por separado.

Aquí está mi guión hasta ahora:

# import stuff
import arcpy

# import HEC-GeoHMS toolbox
arcpy.ImportToolbox(r"C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcToolbox\Toolboxes\GeoHMS Tools.tbx")

# set environments
arcpy.env.workspace = arcpy.GetParameterAsText(0)
arcpy.env.extent = "MAXOF"
arcpy.env.overwriteOutput = "TRUE"

DEM = arcpy.GetParameterAsText(1)
if DEM == '#' or not DEM:
    DEM = "agreedem"  # provide a default value if unspecified

Flow_Dir_Grid = arcpy.GetParameterAsText(2)
if Flow_Dir_Grid == '#' or not Flow_Dir_Grid:
    Flow_Dir_Grid = "Fdr"  # provide a default value if unspecified

Watershed_Layer = arcpy.GetParameterAsText(3)
if Watershed_Layer == '#' or not Watershed_Layer:
    Watershed_Layer = "Watershed"  # provide a default value if unspecified

# (Iterate Feature Selection):
# Create a cursor that cycles through all the values in the Shapefile
rows = arcpy.SearchCursor(Watershed_Layer)

count = 1
#
# #Start the loop
for row in rows:
    # get the value of JOIN_ID for the current row
    desc = row.getValue("Name")

    # print(desc)
    input_name = "Watershed_" + desc
    output_name = "LongestFlowPath_" + desc
    arcpy.AddMessage("The layer name that is input to the longest flow path is: " + input_name)
    arcpy.AddMessage("The resulting output is named: " + output_name)

    arcpy.MakeFeatureLayer_management(Watershed_Layer, input_name, "Name='" + str(desc) + "'")

    out_loc = arcpy.env.workspace + "\\" + output_name
    arcpy.AddMessage("Longest path result output location: " + out_loc)

    # Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
    # The following inputs are layers or table views: "Soucook_watershed_HUC10_DEM_LiDAR", "Fdr", "Watershed_1"
    arcpy.LongestFlowpath_geohms(in_dem_raster=DEM,
                                 in_flowdir_raster = Flow_Dir_Grid,
                                 in_subbasin_features = input_name,
                                 out_longestflowpath_features = out_loc)

    arcpy.AddMessage("Longest Flowpath Iter " +str(count) +" successful!")

    # increment the counter
    row = rows.next
    count = count + 1

El conjunto de datos de prueba que estoy utilizando es una clase de rasgo que tiene dos polígonos, y se crean dos clases de rasgo de polilínea, pero en lugar de que cada clase de rasgo sólo tenga una polilínea para el polígono respectivo, ambas tienen ambas. ¿Cómo puedo asegurarme de que cada clase de característica polilínea resultante sólo tiene una única polilínea, para un único polígono, en ella?

Editar 1:

Acabo de probar este código, que crea clases de características separadas, cada una de ellas con un solo polígono, y luego introduce cada clase de características en la herramienta de recorrido más largo. De alguna manera, sigue creando dentro de cada resultado todas las polilíneas de tantos polígonos como hay en la clase original, ignorando de alguna manera el hecho de que cada clase de característica que se lee sólo tiene un polígono.

# import stuff
import arcpy

# import HEC-GeoHMS toolbox
arcpy.ImportToolbox(r"C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcToolbox\Toolboxes\GeoHMS Tools.tbx")

# set environments
arcpy.env.workspace = arcpy.GetParameterAsText(0)
arcpy.env.extent = "MAXOF"
arcpy.env.overwriteOutput = "TRUE"

DEM = arcpy.GetParameterAsText(1)
if DEM == '#' or not DEM:
    DEM = "agreedem"  # provide a default value if unspecified

Flow_Dir_Grid = arcpy.GetParameterAsText(2)
if Flow_Dir_Grid == '#' or not Flow_Dir_Grid:
    Flow_Dir_Grid = "Fdr"  # provide a default value if unspecified

Watershed_Layer = arcpy.GetParameterAsText(3)
if Watershed_Layer == '#' or not Watershed_Layer:
    Watershed_Layer = "Watershed"  # provide a default value if unspecified

# (Iterate Feature Selection):
# Create a cursor that cycles through all the values in the Shapefile
rows = arcpy.SearchCursor(Watershed_Layer)

input_list = []
output_list = []
lfp_output_list = []

final_output_name = Watershed_Layer + "_LongestFlowPaths_Merged"

count = 1
#
# #Start the loop
for row in rows:
    # get the value of "Name" for the current row
    desc = row.getValue("Name")

    input_name = "Watershed_" + desc
    output_name = "LongestFlowPath_" + desc
    arcpy.AddMessage("The layer name that is input to the longest flow path is: " + input_name)
    arcpy.AddMessage("The resulting output is named: " + output_name)

    arcpy.MakeFeatureLayer_management(Watershed_Layer, input_name, "Name='" + str(desc) + "'")

    out_loc = arcpy.env.workspace + "\\" + input_name
    arcpy.AddMessage("Longest path result output location: " + out_loc)

    # Write the selected features to a new featureclass
    arcpy.CopyFeatures_management(input_name, input_name)

    input_list.append(input_name)
    output_list.append(output_name)
    # increment the counter
    row = rows.next
    count = count + 1

arcpy.AddMessage("Done creating features - now generating longest flow paths.")

for input_watershed, output_nam in zip(input_list,output_list):

    out_loc = arcpy.env.workspace + "\\Layers\\" + output_nam

    lfp_output_list.append(out_loc)

    # Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
    # The following inputs are layers or table views: "Soucook_watershed_HUC10_DEM_LiDAR", "Fdr", "Watershed_1"
    arcpy.LongestFlowpath_geohms(in_dem_raster=DEM,
                                  in_flowdir_raster=Flow_Dir_Grid,
                                  in_subbasin_features=input_watershed,
                                  out_longestflowpath_features=out_loc)

    arcpy.AddMessage("Longest Flowpath Iter " + str(count) + " successful!")

Imagen de las dos trayectorias de flujo más largas que resultan de pasar una cuenca hidrográfica (3789) por la herramienta, a pesar de que la capa sólo tiene un polígono. La capa original a partir de la cual se creó tiene ambos polígonos (3789 y 3791), y la herramienta parece estar enganchada a esto.

[![introduzca la descripción de la imagen aquí][1]]

Editar 2:

Intentando una táctica diferente - en lugar de intentar ejecutar la herramienta usando python en una clase de característica que contiene cientos de polígonos, estoy intentando ejecutar la herramienta en un script por lotes en cientos de clases de características separadas, cada una de las cuales contiene un solo polígono.

En la foto de abajo, Cuenca_3 es la clase de característica original que contiene cientos de polígonos y el Cuenca_de_segmentación_XXXX son las clases de características individuales divididas que contienen un solo polígono.

si ejecuto la herramienta en una sola de estas cuencas, y copio la acción del Resultados como un script de python, obtengo este trozo de código:

# Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
# The following inputs are layers or table views: "Watershed_HUC10_DEM_LiDAR", "fdr"
arcpy.LongestFlowpath_geohms(in_dem_raster="Watershed_HUC10_DEM_LiDAR",
in_flowdir_raster="fdr",
in_subbasin_features="C:/Users/XXXXX/Desktop/XXXXXX/LiDAR_HUC/LiDAR_HUC.gdb/Watersheds_Separate/Watershed_split_3838",
out_longestflowpath_features="C:\Users\XXXXX\Desktop\XXXXXX\LiDAR_HUC\LiDAR_HUC.gdb\Watersheds_Separate\Watershed_split_3838_LFP")

Fíjate en las barras inclinadas que hay entre la entrada y la salida. ¿Importa eso?

Así que, usando eso como base, desarrollé este guión:

import arcpy
# import HEC-GeoHMS toolbox
arcpy.ImportToolbox(r"C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcToolbox\Toolboxes\GeoHMS Tools.tbx")

#set environments
arcpy.env.workspace=arcpy.GetParameterAsText(0)
arcpy.env.extent="MAXOF"
arcpy.env.overwriteOutput="TRUE"

# Local variables:
DEM = arcpy.GetParameterAsText(1)    
Flow_Dir_Grid = arcpy.GetParameterAsText(2)

featureList = arcpy.ListFeatureClasses(wild_card="Watershed_split_*")

lfp_list = []

for feat in featureList:

    name = feat + "_LFP"
    in_feat = r"C:/Users/XXXXX/Desktop/XXXXXX/LiDAR_HUC/LiDAR_HUC.gdb/Watersheds_Separate/" + feat
    out_loc = "C:\\Users\\XXXXX\Desktop\\XXXXXX\\LiDAR_HUC\\LiDAR_HUC.gdb\\Watersheds_Separate\\" + name
    arcpy.AddMessage("Out:" + out_loc)

    # Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
    arcpy.LongestFlowpath_geohms(in_dem_raster=DEM,
                                 in_flowdir_raster=Flow_Dir_Grid,
                                 in_subbasin_features=in_feat,
                                 out_longestflowpath_features=out_loc)

    lfp_list.append(name)

arcpy.AddMessage("Done creating features.")

Este script parece funcionar a veces y fallar otras, y todavía no puedo averiguar por qué - una vez lo ejecuto y procede correctamente, pero otra vez lo ejecuto y falla en la segunda clase de característica. ¿Puede alguien verificar si es posible ejecutar esta herramienta utilizando un script de python en apoyo de un proceso por lotes?

Edita 3:

El script anterior parece funcionar para mí, siempre y cuando: 1. 1. Mantengo la división de las cuencas hidrográficas y la generación de las rutas de flujo más largas por separado 2. Ejecuto el script dentro de una clase de característica 3. 3. Escribo explícitamente toda la ruta de entrada y la ruta de salida, excepto el nombre del archivo creado iterativamente. 4. Utilicé la misma configuración de la barra diagonal hacia delante o hacia atrás que se utiliza en el fragmento de python copiado del Resultados ventana. 5. No dejo que mi ordenador duerma.

Answer 1

¡Me imagino que el problema es que con esta línea la entrada es también la salida!

arcpy.CopyFeatures_management(input_name, input_name)

Si fuera yo, crearía su objeto de capa a través de MakeFeatureLayer y le daría un nombre completamente diferente como "TempLayer". Esta sería la entrada de la herramienta de copia de características. La salida de su herramienta de copia de características probablemente debería ser out_loc ya que ese es el nombre de la ruta completa.

Answer 2

Si ejecuto la herramienta en una sola de estas cuencas, y copio la acción de la ventana de resultados como un script de python, obtengo este trozo de código:

    # Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
    # The following inputs are layers or table views: "Watershed_HUC10_DEM_LiDAR", "fdr"
    arcpy.LongestFlowpath_geohms(in_dem_raster="Watershed_HUC10_DEM_LiDAR",
    in_flowdir_raster="fdr",
    in_subbasin_features="C:/Users/XXXXX/Desktop/XXXXXX/LiDAR_HUC/LiDAR_HUC.gdb/Watersheds_Separate/Watershed_split_3838",
    out_longestflowpath_features="C:\Users\XXXXX\Desktop\XXXXXX\LiDAR_HUC\LiDAR_HUC.gdb\Watersheds_Separate\Watershed_split_3838_LFP")

**Notice the differently facing slashes between the input and output!** Does that matter?

Así que, usando eso como base, desarrollé esto el script de abajo parece funcionar para mí siempre y cuando:

1. Mantengo separadas la división de las cuencas y la generación de las vías de flujo más largas
2. Ejecuto el script dentro de una clase de característica
3. Escribo explícitamente toda la ruta de entrada y la de salida, excepto el nombre del archivo creado de forma iterativa.
4. Utilicé la misma configuración de barra diagonal hacia delante o hacia atrás que se utiliza en el fragmento de python copiado del Resultados ventana.
5. No dejo que mi ordenador duerma.

Código:

import arcpy, os, traceback, sys, time
# import HEC-GeoHMS toolbox
arcpy.ImportToolbox(r"C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcToolbox\Toolboxes\GeoHMS Tools.tbx")

# function for error formatting/printing
def showPyMessage():
    arcpy.AddMessage(str(time.ctime()) + " - " + message)

#set environment variables
arcpy.env.workspace=arcpy.GetParameterAsText(0)
arcpy.env.overwriteOutput="TRUE"

# Tool variables:
DEM = arcpy.GetParameterAsText(1)
Flow_Dir_Grid = arcpy.GetParameterAsText(2)

# list of current feature classes
featureList = arcpy.ListFeatureClasses(wild_card="Watershed_split_*")
# empty list of finished feature classes - watersheds where longest flow path has been calculated
lfp_list = []
# filling list of finished feature classes - watersheds where longest flow path has been calculated
for feat in featureList:
    if '_LFP' in feat:
        lfp_list.append(feat)

# printing the previously finished list
arcpy.AddMessage('Previously finished list: ' + ', '.join(lfp_list))        

# subtracting out inputs (and outputs) that already have been finished in a previous run (maybe this script crashed?)
for fin_feat in lfp_list:
    arcpy.AddMessage('Removing ' + fin_feat + ' from processing list.')
    featureList.remove(fin_feat)
    arcpy.AddMessage('Removing ' + fin_feat.replace('_LFP','') + ' from processing list.')
    featureList.remove(fin_feat.replace('_LFP',''))

# printing the previously finished list
arcpy.AddMessage('Updated featureList: ' + ', '.join(featureList) )        

arcpy.AddMessage('Starting batch LFP cals.')

# loop that makes sure the process continues forever - AKA until all inputs have a matched output
while featureList:
    # for each input in the trimmed feature list
    for feat in featureList:
    try:
        # setting up the name
        name = feat + "_LFP"
        # the input and output string need (as far as I can tell) to be pretty much hardcoded in the manner shown below
        in_feat = r"C:/XXXXX/XXXXX/XXXXX.gdb/Watersheds_Separate/" + feat
        out_loc = "C:\\\\XXXXXX\\XXXXXX\\XXXXX.gdb\\Watersheds_Separate\\" + name
        arcpy.AddMessage("Out:" + out_loc)

        # Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
        arcpy.LongestFlowpath_geohms(in_dem_raster=DEM,
                                     in_flowdir_raster=Flow_Dir_Grid,
                                     in_subbasin_features=in_feat,
                                     out_longestflowpath_features=out_loc)

        # adding the new output to the finished output list
        lfp_list.append(name)
        # removing the input for which the LFP was successfully calculated from the input list
        arcpy.AddMessage('Removing ' + feat + ' from processing list.')
        featureList.remove(feat)

    except:

        # error reporting
        name = feat + "_LFP"
        arcpy.Delete_management(in_data="C:\\Users\\Thomas.p.taggart\\Desktop\\Ballestero_Model_Learning\\Soucook_River_LiDAR_HUC\\Soucook_River_LiDAR_HUC.gdb\\Watersheds_Separate\\" + name, data_type="FeatureClass")
        message = "\n*** PYTHON ERRORS *** "; showPyMessage()
        message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
        message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()

arcpy.AddMessage("Done creating features.")

# Merging the finished LFPs
arcpy.Merge_management(lfp_list, "Merged_LFP")

# Deleting the old LPFs

for feat in featureList:
    arcpy.Delete_management(in_data=feat, data_type="FeatureClass")