1

I have to generate LONGEST FLOW PATH(polyline) within watershed. What is the procedure to generate this ? After generating longest flow length , I will divide polygon with polyline. I have :

  • watershed layer(polygon)
  • stream layer(polyline)
  • Digital Elevation Model (DEM) layer

I am using .NET and python in ArcGIS 10. I have tried "FlowLength()" http://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/flow-length.htm method , but it gives raster output. I tried to convert this raster dataset to polyline. this was my procedure.

I have attached an image to simplify this problem.

enter image description here

2
  • 1
    Have you tried the ArcHydro toolkit?
    – Nxaunxau
    Commented Jul 11, 2016 at 10:34
  • You will have a problem when it comes to cutting your catchment into the 2 halves. Unless the initiation of the river path is generated upstream of the catchment then it will almost certainly be the case that the rivers upstream end does not cross the catchment boundary and you won't be able to do the cut...
    – Hornbydd
    Commented Jul 11, 2016 at 13:45

1 Answer 1

2

Input: enter image description here

Catchments table:

enter image description here

WORKFLOW TO DERIVE LONGEST FLOW PATHS:

arcpy.PolygonToRaster_conversion(in_features="subcatchments", value_field="ID", out_rasterdataset="D:/Scratch/CATCHM")
arcpy.gp.FocalStatistics_sa("CATCHM", "D:/Scratch/RANGE", "Rectangle 3 3 CELL", "RANGE", "DATA")
arcpy.gp.RasterCalculator_sa("""Con("RANGE" == 0,"CATCHM")""", "D:/Scratch/ISLANDS")
arcpy.gp.RasterCalculator_sa("""Con( ~IsNull("ISLANDS"),"dem")""", "D:/Scratch/FENCED")
arcpy.gp.Fill_sa("FENCED", "D:/Scratch/FILLED", "")
arcpy.gp.FlowDirection_sa("FILLED", "D:/Scratch/FDIR", "NORMAL", "")
arcpy.gp.FlowLength_sa("FDIR", "D:/Scratch/FLEN", "DOWNSTREAM", "")
arcpy.gp.ZonalStatistics_sa("ISLANDS", "VALUE", "FLEN", "D:/Scratch/LENMAX", "MAXIMUM", "DATA")
arcpy.gp.RasterCalculator_sa("""Con("LENMAX" == "FLEN","ISLANDS")""", "D:/Scratch/MAXPNTS")
arcpy.DeleteIdentical_management(in_dataset="SOURCES", fields="GRID_CODE", xy_tolerance="", z_tolerance="0")
arcpy.gp.CostPath_sa("SOURCES", "FDIR", "FDIR", "D:/Scratch/CPATHS", "EACH_CELL", "GRID_CODE")
arcpy.gp.StreamToFeature_sa("CPATHS", "FDIR", "D:/Scratch/LongFPs.shp", "SIMPLIFY")

SPLIT CATCHMENTS (DEM cell size = 1*1):

arcpy.FeatureToLine_management(in_features="subcatchments", out_feature_class="D:/Scratch/ARCS.shp")
arcpy.Merge_management(inputs="ARCS;LongFPs", output="D:/Scratch/MERGED.shp", field_mappings="""ID "ID" true true false 10 Long 0 10 ,First,#,ARCS,ID,-1,-1""")
arcpy.ExtendLine_edit(in_features="MERGED", length="15 Meters", extend_to="EXTENSION")
arcpy.FeatureToPolygon_management(in_features="MERGED", out_feature_class="D:/Scratch/splits.shp", cluster_tolerance="", attributes="ATTRIBUTES", label_features="")
arcpy.SpatialJoin_analysis(target_features="splits", join_features="subcatchments", out_feature_class="D:/Scratch/FINAL.shp", join_operation="JOIN_ONE_TO_ONE", join_type="KEEP_ALL", field_mapping="""Id "Id" true true false 6 Long 0 6 ,First,#,splits,Id,-1,-1;ID_1 "ID_1" true true false 10 Long 0 10 ,First,#,subcatchments,ID,-1,-1""", match_option="WITHIN", search_radius="", distance_field_name="")

OUTPUT LABELLED BY COUNT OF PARTS:

enter image description here

As one can see some manual editing of MERGED lines (no split happened to catchment 4) or FINAL (too many splits of catchment 10) will be needed

7
  • If it solved your problem, mark it please
    – FelixIP
    Commented Jul 23, 2016 at 4:14
  • I have tried your code to arcgis python console with my data set. I have done successfully upto : .... .... arcpy.gp.RasterCalculator_sa("""Con("LENMAX" == "FLEN","ISLANDS")""", "D:/Scratch/MAXPNTS") I can not understand your input "SOURCES" at this execution and got exception : arcpy.DeleteIdentical_management(in_dataset="SOURCES", fields="GRID_CODE", xy_tolerance="", z_tolerance="0"). Is "SOURCES" input my "watershed", "DEM" or "Stream" ?? Thanks.
    – Monir
    Commented Jul 23, 2016 at 5:59
  • There is a missing line, sorry. Raster MAXPOINTS has to be converted to points (raster to points) and called SOURCES
    – FelixIP
    Commented Jul 23, 2016 at 6:06
  • Many many thanks . It works perfectly. If I have any inquiry , contact with you later. :)
    – Monir
    Commented Jul 23, 2016 at 6:34
  • I have faced trouble in different data source. I have another data set and when I tried this code statement : arcpy.gp.ZonalStatistics_sa("ISLANDS", "VALUE", "FLEN", "D:/Scratch/LENMAX", "MAXIMUM", "DATA") ; It gives this error : arcpy.gp.ZonalStatistics_sa("ISLANDS", "VALUE", "FLEN", "D:/Scratch/LENMAX", "MAXIMUM", "DATA") Runtime error <class 'arcgisscripting.ExecuteError'>: ERROR 000864: Input raster or feature zone data: The input is not within the defined domain. ERROR 000367: Invalid GP data type . Have you any suggestion ??
    – Monir
    Commented Jul 23, 2016 at 6:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.