Going back in for loop at the end of a while loop in ArcPy?

Question

I'm quite novice in Python and ArcPy.

I have a river network with Strahler orders (starting from 1). I want to create a new network where the new streams will start from order 1 till the outlet. (so, I'm looking for the affluent of the first stream, then the next affluent, then the next ...) I used a function that compares the coordinates of the point at the end of the stream (XY_end) with the start at the point of the other streams (XY_start) in order to find the next affluent, and it works perfectly. It's certainly not optimal, but it works (after 5 days on it...).

My problem is with the while loop : When I reach the last affluent of the chain, it stops. I'm not able to go back in the for loop to get the next stream with Strahler order = 1.

What did I do wrong and how can I correct it?

The gdb : [zip file of the gdb with input stream network and result][1]

The message I get :

ce n°: 2 d'ordre 1 et son affluent 272 Nouvel affluent [272, u'2', [335323.5, 5180614.5], [335278.5, 5180569.5]] ce n°: 272 d'ordre 2 et son affluent 258 Nouvel affluent [258, u'2', [335278.5, 5180569.5], [335222.5, 5180591.5]] ce n°: 258 d'ordre 2 et son affluent 30 Nouvel affluent [30, u'2', [335222.5, 5180591.5], [335156.5, 5180615.5]] ce n°: 30 d'ordre 2 et son affluent 34 Nouvel affluent [34, u'2', [335156.5, 5180615.5], [335087.5, 5180709.5]] ce n°: 34 d'ordre 2 et son affluent 248 Nouvel affluent [248, u'2', [335087.5, 5180709.5], [335085.5, 5180712.5]] ce n°: 248 d'ordre 2 et son affluent 195 Nouvel affluent [195, u'2', [335085.5, 5180712.5], [335085.2829, 5180712.7171]] Nouvel affluent None Runtime error Traceback (most recent call last): File "", line 101, in File "", line 81, in main File "", line 94, in GetIDaffluent_1 TypeError: 'NoneType' object has no attribute 'getitem'

[![In green the whole chain of new streams with the ID (2) of the first stream (strahler 1][2]][2]

    def main():
        import arcpy
        import os
        arcpy.env.overwriteOutput = True



    fc = r'C:\network.gdb\ZP_Cours_Eau_3D'
        fc_out = r'C:\network.gdb\draft_affluents'

        # # output fields
        fld_ce = "ID_Cours"

        # # create output fc
        sr = arcpy.Describe(fc).spatialReference
        ws, fc_name = os.path.split(fc_out)
        arcpy.CreateFeatureclass_management(ws, fc_name, "POLYLINE", None, "ENABLED", "ENABLED", sr)

        # # add fields in fc
        coord = [["X_start", "!SHAPE.firstPoint.X!"], ["Y_start", "!SHAPE.firstPoint.Y!"], ["X_end", "!SHAPE.lastPoint.X!"], ["Y_end", "!SHAPE.lastPoint.Y!"]]  
        for pt in coord:  
            # # print("adding field {0}...".format(pt[0]))  
            arcpy.AddField_management(fc, pt[0], "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")  
            # # print("  calcing field {0} to be {1}".format(pt[0], pt[1]))  
            arcpy.CalculateField_management(fc, pt[0], pt[1], "PYTHON_9.3" ) 

        # # add fields in fc_out
        arcpy.AddField_management(fc_out, fld_ce, "DOUBLE")

        # # define variables
        # # Trouver la valeur max des ordres de strahler des cours d'eau
        max_ost = arcpy.da.SearchCursor(fc, "strahler", "strahler IS NOT NULL", sql_clause = (None, "ORDER BY strahler DESC")).next()[0]
        # print max_ost
        lst_coord = []
        coord = []
        Cursor = arcpy.SearchCursor(fc)
        # # lister les cours d'eau avec leurs caractéristiques
        for row in Cursor:
            value = row.OBJECTID
            x_start = row.X_start
            y_start = row.Y_start
            x_end = row.X_end
            y_end = row.Y_end
            ost = row.getValue("strahler")
            coord = [value, ost, [x_start, y_start], [x_end, y_end]]
            lst_coord.append(coord)
        # print lst_coord
        lst = []
        IDce = 0
        ost = 0
        pt_start = []
        pt_end = []
        New_ce_1 = []
        index = 0
        #parcourir chacun des cours d'eau
        for ce in lst_coord:
            fld_oid = arcpy.Describe(fc).OIDFieldName
            where = "{0} = {1}".format(fld_oid, ce[0])
            polyline = arcpy.da.SearchCursor(fc, ('SHAPE@'), where).next()[0]
            flds = ('SHAPE@', fld_ce)
            #Si le cours d'eau est de niveau 1, le copier en sortie et chercher ses affluents
            #en attribuant le ID du cours d'eau de niveau 1 à tous les affluents
            with arcpy.da.InsertCursor(fc_out, flds) as curs:
                lst = lst_coord[index]
                IDce = lst[0]
                ost = lst[1]
                pt_start = lst[2]
                pt_end = lst[3]

                if ost == u'1':
                    curs.insertRow((polyline, IDce))
                    New_ce_1 = GetIDaffluent_1(lst_coord, ce)

                    while ce is not None:
                        where_1 = "{0} = {1}".format(fld_oid, ce[0])
                        polyline_1 = arcpy.da.SearchCursor(fc, ('SHAPE@'), where_1).next()[0]
                        curs.insertRow((polyline_1, IDce))
                        print "Nouvel affluent", New_ce_1

                        ce = New_ce_1
                        New_ce_1 = GetIDaffluent_1(lst_coord, ce)
                        #break

            index += 1 

    def GetIDaffluent_1(lst_coord, ce):
        # trouver l'affluent à ce dans la liste des cours d'eau (pt_start affluent = pt_end ce)
        cnt_1 = 0
        aff_1=[]
        for aff_1 in lst_coord:
            aff_1 = lst_coord[cnt_1]
            cnt_1 += 1
            # Comparer les coordonnées end du cours d'eau aux coordonnées start de l'affluent
            if aff_1[2] == ce[3]:
                #aff_1[1] = ce[1]
                print "ce n°:", ce[0], "d'ordre", ce[1], "et son affluent", aff_1[0]
                    # print "ID_ce :", ce[0], "IDaff_1:", aff_1[0],"x,y start:", aff_1[2]
                return aff_1

    if __name__ == '__main__':
        main()
    ```


  [1]: https://www.dropbox.com/s/if94ggz5yduiyap/network.gdb.zip?dl=0
  [2]: https://i.stack.imgur.com/w76D7.png

Answer 1

You can use networkx package to solve your problem:

NetworkX is a Python package for the creation, manipulation, and study of the structure, dynamics, and functions of complex networks.

I had to install it using C:\Python27\ArcGIS10.7\python.exe python -m pip install networkx in windows cmd prompt.

import networkx as nx
import arcpy

arcpy.env.workspace = r'X:\radf\RADF_test.gdb'
fc = r'ZP_Cours_Eau_3D'
fc_out = r'ZP_Cours_Eau_3D_out'

#Create output fc
arcpy.CopyFeatures_management(in_features=fc, out_feature_class=fc_out)
arcpy.TruncateTable_management(fc_out)
arcpy.AddField_management(in_table=fc_out, field_name='pathid', field_type='LONG')
fieldlist = [f.name for f in arcpy.ListFields(fc) if f.type not in ('Geometry','OID') and not f.name.startswith('Shape')]

#Store all original data in a dictionary with oid as key and a list of all attributes and geometry as value
data = {d[0]:d[1:] for d in arcpy.da.SearchCursor(fc,['OID@']+fieldlist+['SHAPE@'])}

#Create directional graph
fields = ['X_start','Y_start','X_end','Y_end','OID@']
L = [row for row in arcpy.da.SearchCursor(fc, fields)]
L = [[(r[0],r[1]), (r[2],r[3]), {'OID':r[4]}] for r in L] #[[1,2,3,4,5].. -> [[(1,2),(3,4),{'OID':5}]...
G=nx.DiGraph()
G.add_edges_from(L)

#If you want to plot the graphs:
#import matplotlib.pyplot as plt 
#nx.draw(G,with_labels = True)
#plt.show()

icur = arcpy.da.InsertCursor(fc_out, fieldlist+['SHAPE@']+['pathid'])

#Iterate over each watersystem/subgraph and list all objectids of paths from startnodes to endnode
pid = 1 #Variable to give each flowpath an id
for watersystem in nx.weakly_connected_component_subgraphs(G):
    startnodes = [x for x in watersystem.nodes() if watersystem.in_degree(x)==0]
    endnode = [x for x in watersystem.nodes() if watersystem.out_degree(x)==0][0]
    for startnode in startnodes:
        sub = []
        a = nx.all_simple_paths(watersystem, startnode, endnode) #Find path from startnode to endnode
        b = [i for i in a][0] #List nodes in path
        for node1, node2 in zip(b,b[1:]):
            sub.append(watersystem.get_edge_data(node1,node2)['OID']) #Extract OID from each edge
        for oid in sub:
            icur.insertRow(data[oid]+(pid,)) #Use dictionary to insert correct row
        pid+=1
del icur

Example of output with 6 startnodes:

Field pathid will hold individual ids for each flowpath:

Example of one graph:

Answer 2

I think I finally got it! I didn't know about reset or continue, and it works now, using "while ... else: continue" . I had to correct other mistakes too. Now I can work on optimisation (dissolve before creating the new streams). Here are the corrections :

def main():
   import arcpy
   import os
   arcpy.env.overwriteOutput = True

   fc = r'C:\RADF_1.gdb\ZP_Cours_Eau_3D'
   fc_out = r'C:\RADF_1.gdb\draft_affluents'

   # # output fields
   fld_ce = "ID_Cours"

   # # create output fc
   sr = arcpy.Describe(fc).spatialReference
   ws, fc_name = os.path.split(fc_out)
   arcpy.CreateFeatureclass_management(ws, fc_name, "POLYLINE", None, "ENABLED", "ENABLED", sr)

   # # add fields in fc
   # coord = [["X_start", "!SHAPE.firstPoint.X!"], ["Y_start", "!SHAPE.firstPoint.Y!"], ["X_end", "!SHAPE.lastPoint.X!"], ["Y_end", "!SHAPE.lastPoint.Y!"]]  
   # for pt in coord:  
       # # print("adding field {0}...".format(pt[0]))  
       # arcpy.AddField_management(fc, pt[0], "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")  
       # # print("  calcing field {0} to be {1}".format(pt[0], pt[1]))  
       # arcpy.CalculateField_management(fc, pt[0], pt[1], "PYTHON_9.3" ) 

   # # add fields in fc_out
   arcpy.AddField_management(fc_out, fld_ce, "DOUBLE")

   # # define variables
   # # Trouver la valeur max des ordres de strahler des cours d'eau
   #max_ost = arcpy.da.SearchCursor(fc, "strahler", "strahler IS NOT NULL", sql_clause = (None, "ORDER BY strahler DESC")).next()[0]
   # print max_ost
   lst_coord = []
   coord = []
   Cursor = arcpy.SearchCursor(fc)
   # # lister les cours d'eau avec leurs caractéristiques
   for row in Cursor:
       value = row.OBJECTID
       x_start = row.X_start
       y_start = row.Y_start
       x_end = row.X_end
       y_end = row.Y_end
       ost = int(row.getValue("strahler"))
       coord = [value, ost, [x_start, y_start], [x_end, y_end]]
       lst_coord.append(coord)
   # print lst_coord
   IDce = 0
   ost = 0
   New_ce_1 = []

   #parcourir chacun des cours d'eau
   for ce in lst_coord:
       fld_oid = arcpy.Describe(fc).OIDFieldName
       where = "{0} = {1}".format(fld_oid, ce[0])
       polyline = arcpy.da.SearchCursor(fc, ('SHAPE@'), where).next()[0]
       flds = ('SHAPE@', fld_ce)
       #Si le cours d'eau est de niveau 1, le copier en sortie et chercher ses affluents
       #en attribuant le ID du cours d'eau de niveau 1 à tous les affluents
       with arcpy.da.InsertCursor(fc_out, flds) as curs:
           ost = ce[1]
           if ost == 1:
               print "cours d'eau initial", ce[0], "ordre de strahler", ost
               curs.insertRow((polyline, ce[0]))
               New_ce_1 = GetIDaffluent_1(lst_coord, ce)
               IDce = ce[0]
               while New_ce_1 is not None:
                   print "Nouvel affluent", New_ce_1
                   ce = New_ce_1
                   New_ce_1 = GetIDaffluent_1(lst_coord, ce)
                   where_1 = "{0} = {1}".format(fld_oid, ce[0])
                   polyline_1 = arcpy.da.SearchCursor(fc, ('SHAPE@'), where_1).next()[0]
                   curs.insertRow((polyline_1, IDce))
               else:
                   continue                    

def GetIDaffluent_1(lst_coord, ce):
   # trouver l'affluent à ce dans la liste des cours d'eau (pt_start affluent = pt_end ce)
   cnt_1 = 0
   aff_1=[]
   for aff_1 in lst_coord:
       aff_1 = lst_coord[cnt_1]
       cnt_1 += 1
       # Comparer les coordonnées end du cours d'eau aux coordonnées start de l'affluent
       if aff_1[2] == ce[3]:
           #aff_1[1] = ce[1]
           print "ce n°:", ce[0], "d'ordre", ce[1], "et son affluent", aff_1[0]
               # print "ID_ce :", ce[0], "IDaff_1:", aff_1[0],"x,y start:", aff_1[2]
           return aff_1

if __name__ == '__main__':
   main()

Answer 3

I didn't read the entire script, but you can use continue to break out of a loop iteration and continue with the next one. Not sure if this is what you're looking for.

    while ce is not None:
        where_1 = "{0} = {1}".format(fld_oid, ce[0])
        polyline_1 = arcpy.da.SearchCursor(fc, ('SHAPE@'), where_1).next()[0]
        curs.insertRow((polyline_1, IDce))
        print "Nouvel affluent", New_ce_1
        if some_condition:
            continue    # skips rest of code in loop, starts back at top of loop
        ce = New_ce_1
        New_ce_1 = GetIDaffluent_1(lst_coord, ce)