1

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
  • Yes, the green one is the output I want. In fact, the green line is composed of all the parts (affluents) found from the first stream (strahler 1) till the last one, but they have a new common attribute, so I can dissolve them later.. And yes, I want to create all of the other new streams the same way. for a whole watershed. I know it's a lot of redundance, but I need the whole streams frome start to end. I understand the point of .reset() I tried Cursor.reset() in many ways, but I'm not able to change the result :-/ – Julie Aug 21 at 0:40
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:

enter image description here

Field pathid will hold individual ids for each flowpath: enter image description here

Example of one graph: enter image description here

  • I managed to install the networkx package, but ArCIG uses python 2.7, and networkx needs python 3. I'm not sure about what to do to deal with this issue – Julie Aug 27 at 13:04
  • I used the same pathway : double clic on python.exe in the file C:\Python27\ArcGIS10.7\ then enter -m pip install networkx on the command line. and I get a invalid syntax : "File "<stdin>", line 1 -m pip install networkx SyntaxError: invalid syntax – Julie Aug 27 at 14:28
  • Hi Bera, the code runs perfectly and really faster than mine on the test zone (the one you used) 10 sec. vs 50 sec!! But I have a 'IndexError: list index out of range' on my whole watershed (line 39). I don't understand. – Julie Aug 27 at 15:09
  • Sorry, no. Line 39 is endnode = [x for x in watersystem.nodes() if watersystem.out_degree(x)==0 and watersystem.in_degree(x)==1][0] #Find endnode for subgraph As it worked on the test zone, I add two lines to do the dissolve on "pathid". – Julie Aug 27 at 15:35
  • ? sorry ... not supposed to. Try dropbox.com/s/ur0t0r6dur559wa/RADF_2.gdb.zip?dl=0 – Julie Aug 27 at 16:08
1

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()
0

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)
  • Thanks to both of you, but unfortunately, .reset() or "continue" gives me the same result. It stops at the end of the while loop .... – Julie Aug 21 at 0:42
  • I tried some del too ... same result. I'm not sure to use these different methods the right way, but I'll make some more tries – Julie Aug 21 at 0:55
  • It also looks like the last line shouls be ce = GetIDaffluent_1(lst_coord, ce), not New_ce_1, if I'm understanding correctly. May want to add a check to see if ce is None as well. – recurvata Aug 21 at 14:52
  • In fact, no. I need the New_ce_1 : i'm searching for its XY_start = ce XY end (so, New_ce_1 is the affluent of ce). Then, I use ce = New_ce_1 so my last affluent becomes my "ce", and so on. I don't know if my explanation is clear ? – Julie Aug 21 at 15:32
  • Possibly the error occurs elsewhere in the code, such as: if ost == u'1':, preventing the code from reaching the while loop on subsequent iterations? Your English, and grasp of arcpy, seem excellent as well. – recurvata Aug 21 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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