Using Longest Flow Path tool from HEC-GeoHMS in loop with Arcpy

Question

I have a feature class made of hundreds of polygons, that I need to determine the longest flow path for using HEC-GeoHMS. As anyone who has used this tool knows, you can't just pass the toolset the entire feature class if any of the watersheds overlap - it will not correctly calculate the longest flowpath for overlapping basins.

So, I am trying to setup a script that iterates through each record in my feature class and generates a longest flow path result for each record separately.

Here is my script so far:

# 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

The test data set I am using is a feature class that has two polygons, and two polyline feature classes are created, but instead of each feature class just having the one polyline for the respective polygon, they both have both. How can I ensure that each resulting polyline feature class only has a single polyline, for a single polygon, in it?

Edit 1:

I have just tried this code, which creates separate feature classes, each having only one polygon, and then feeds each feature class into the Longest Flowpath Tool. Somehow it still creates within each result all the polylines from as many polygons as are in the original class, somehow ignoring the fact that each single feature class it is read only has one polygon.

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

Picture of the two longest flowpaths that result from running one watershed (3789) through the tool, despite the layer only having a single polygon. The original layer it was created from has both polygons (3789 and 3791), and the tool seems to be latching onto this.

[![enter image description here][1]][1]

Edit 2:

Trying a different tack - rather than trying to run the tool using python on a feature class containing hundreds of polygons, I am trying to run the tool in a batch script on hundreds of separate features class, each containing only a single polygon.

In the pic below, Watershed_3 is the original feature class containing hundreds of polygons and the Watershed_split_XXXX are the individual split feature classes containing a single polygon.

if i run the tool on a single of these watersheds, and copy the action from the Results window as a python script, I get this piece of code:

# 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?

So, using that as a basis, I developed this script:

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

This script seems to work sometimes and fail others, and I can't figure out why yet - one time I run it and it proceeds correctly, yet another time I run it it fails on the second feature class. Can anyone verify whether it is possible to run this tool using a python script in support of a batch process?

Edit 3:

The above script seems to work for me as long as: 1. I keep the splitting of watersheds and the generation of longest flow paths separate 2. I run the script within a feature class 3. I explicitly write out the entire input path and output path - except for the iteratively created file name. 4. I used the same forward or backward slash setup as is used in the python snippet as copied form the Results window. 5. I don't let my computer sleep.

Answer 1

I would imagine the problem is with this line the input is also the output!

arcpy.CopyFeatures_management(input_name, input_name)

If it were me I would create your layer object via MakeFeatureLayer and give it a completely different name such "TempLayer". This would then be what is the input into the copy features tool. The output of your copy features tool should probably be out_loc as that is what is the whole path name.

Answer 2

If I run the tool on a single of these watersheds, and copy the action from the Results window as a python script, I get this piece of code:

    # 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?

So, using that as a basis, I developed this the script below seems to work for me as long as:

1. I keep the splitting of watersheds and the generation of longest flow paths separate
2. I run the script within a feature class
3. I explicitly write out the entire input path and output path - except for the iteratively created file name.
4. I used the same forward or backward slash setup as is used in the python snippet as copied form the Results window.
5. I don't let my computer sleep.

Code:

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