Using Hec-GeoHMS Longest Flowpath Tool within For loop
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!")
output_list.append(output_name)
# 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:**
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.
