Using Hec-GeoHMS Longest Flowpath Tool within For loop

1741
2
01-26-2017 01:06 PM
Labels (1)
by Anonymous User
Not applicable

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.


 

Tags (3)
0 Kudos
2 Replies
IanMurray
Frequent Contributor

Could you please format your code so it is a bit more legible on the forum?

https://community.esri.com/blogs/dan_patterson/2016/08/14/script-formatting

0 Kudos
by Anonymous User
Not applicable

Edited.  Sorry about that.

0 Kudos