Confluence Retirement

In an effort to consolidate USGS hosted Wikis, myUSGS’ Confluence service is scheduled for retirement on January 27th, 2023. The official USGS Wiki and collaboration space is now SharePoint. Please migrate existing spaces and content to the SharePoint platform and remove it from Confluence at your earliest convenience. If you need any additional information or have any concerns about this change, please contact Thank you for your prompt attention to this matter.
Skip to end of metadata
Go to start of metadata
Contributors: Roland Viger, Mike Wieczorek, Curtis Price
Notes: will not find closed basins above the seed feature (i.e., flowlines where TOCOMID == 0 or that are above a flowline where TOCOMID == 0). works on one NHDPlus Production Unit at a time. This script does not make basins that cross NHDPlus Production Units. This routine produces a shapefile "x<seed>.shp" in the arcpy.env.workspace. has been designed to isolate the algorithmic "smarts" from all the data handling specifics, like path names and field names that are specific to a particular user or data set. If you spot ways to make more generic, please let us know! 
The second code block,, is example of code for iterating over a set of seeds based on NHDPlus v2 PlusFlow.dbf and a catchments.shp files for a single production routine, running the first code block (presumed to be in a file called ""). Note that it passes in all the pieces needed, including an instance of arcpy.


seed = bottom-most feature identifier (i.e., where to start!)

networkArr = an array of tuples with the (id, idDown) in each element

id = feature identifier FROMCOMID (case sensitive)
idDown = downstream feature identifier TOCOMID (case sensitive)
catchments = ArcGIS feature class of catchment geometry, with features identified by COMID label
arcpy = ESRI arcpy python object
def findUpstreamFlowlines(currComid, comidList, networkArr, id, idDown):
    inFlowlines = networkArr[networkArr[idDown] == currComid][id]
        for inFlowline in inFlowlines:
            print currComid, inFlowline
            comidList = findUpstreamFlowlines(inFlowline, comidList,
                                              networkArr, id, idDown)
    return comidList

def extractUpstreamCatchments(seed, comidList, catchments, arcpy):
    print seed
    comidList1 = ",".join([str(k) for k in comidList])
    whereClause = "FEATUREID IN ({})".format(comidList1)
    arcpy.SelectLayerByAttribute_management(catchments, '', whereClause)
    if not arcpy.ListFields(catchments, 'COMID'):
        arcpy.AddField_management(catchments, 'COMID', 'LONG')
    arcpy.CalculateField_management(catchments, 'COMID', '{}'.format(seed))
    print "Making basin shapefile for COMID {} ".format(seed)
    arcpy.Dissolve_management(catchments, 'x{}'.format(seed),
                              ["COMID"], '', '', '')
    print "Finished making basin for COMID {} ".format(seed)

def main(seed, networkArr, id, idDown, catchments, arcpy):
        comidList = findUpstreamFlowlines(seed, list(), networkArr, id, idDown)
        extractUpstreamCatchments(seed, comidList, catchments, arcpy)

if __name__ == '__main__':
    import sys
    main(sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], 
         sys.argv[5], sys.argv[6])


Here is the "batch" code for running repeatedly. It is responsible for starting arpy, reading the major input data sets, and calling for each see COMID.
def main():
    import upstreamCats
    import arcpy
    import numpy as np
    arcpy.env.overwriteOutput = True
    # feel free to invent some paired lists to deal w/NHPlus' 
    puNum = '15'
    state = 'CO'
    wks = 'M:/BasinMaker/PU' + puNum
    arcpy.env.workspace = wks

    outletComidsDBF = 'M:/Requests/HeadWater_Basins_PU' + puNum + '.dbf'
    outletComids = arcpy.da.TableToNumPyArray(outletComidsDBF, ['COMID'])
    plusFlowTbl =  r'M:/NHDPlus2/PU' + puNum + '/NHDPlusV21_' + state + '_' + \
                puNum + '_NHDPlusAttributes_02/NHDPlusCO/NHDPlus' + puNum + \
    id = 'FROMCOMID'
    idDown = 'TOCOMID'
    catchmentsShp = 'M:/NHDPlus2/PU' + puNum + '/NHDPlusCatchment/catchment.shp'
    arcpy.MakeFeatureLayer_management(catchmentsShp, 'catchments')
    networkArr = arcpy.da.TableToNumPyArray(plusFlowTbl, [id, idDown])
    networkArr = np.array([i for i in networkArr if i[0] > 0])

    for outletComid in outletComids:
        seed = outletComid[0]
        upstreamCats.main(seed, networkArr, id, idDown, 'catchments', arcpy)

if __name__ == '__main__':

Code for making basins that are only within a NHDPlus version 2 production unit.

Roland Viger Mike Wieczorek Curtis Price
won't find closed basins above the seed feature (i.e., flowlines where
TOCOMID == 0 or that are above a flowline where TOCOMID == 0)

myNetworkFeatureClass = ESRI feature class Flowlines with PlusFlow.dbf attached (COMID in Flowline joins to FROMCOMID in dbf file)
id = feature identifier COMID (case sensitive)
idDown = downstream feature identifier TOCOMID (case sensitive)
seed = bottom-most feature identifier (i.e., where to start!)
import numpy as np
import sys
import os
import traceback
import arcpy
from arcpy import env
tv = "tv"
wks is the workspace where you want the output basin files to be made
wks = r"Q:/Requests/Gilliom/BasinMaker/PU02"
env.workspace = wks
# Create the search cursor
rows = arcpy.SearchCursor(r"Q:/Requests/GIlliom/HeadWaterBasins/HeadWater_Basins_PU02c.dbf")

# Call to read the first row
row =

# Start a loop that will exit when there are no more rows available
while row:

    # Do something with the values in the current row
    print row.COMID

    env.overwriteOutput = True
    here's the brains. recursive/calls itself.
    def findUpstreamFlowlines(currComid, comidList):


        inFlowlines = networkArr[networkArr[idDown] == currComid][id]
            for inFlowline in inFlowlines:
                comidList = findUpstreamFlowlines(inFlowline, comidList)

        return comidList

    below is just the set up. get numpy array of network's relevant fields, specify
    start point, find upstream features, spit out results.
    # these four vars could be command line args
    myNetworkFeatureClass = r'Q:/NHDPlusV2/NHDPlusMA/NHDPlus02/NHDPlusAttributes/PlusFlow.dbf'

    id = 'COMID'
    idDown = 'TOCOMID'
    at this point I have this hard coded in, seed is the COMID of the Headwater basin to be made,
    I need to make a loop that reads in the CONUS HeadWater file of COMIDs
    seed = row.COMID

    # networkARR is the arrayof upstream COMIDs that gets dumped into a list which I temporarily call BasinMaker.pickle.
        networkArr = np.load('Q:/Requests/Gilliom/BasinMaker/NHD02.pickle')
        import arcpy
        myNetworkFeatureClass = r'Q:/Requests/Gilliom/BasinMaker/NHD02.shp'
        networkArr = arcpy.da.TableToNumPyArray(myNetworkFeatureClass, [id, idDown])

    # create empty container for results
    comidList = list()

    #do the work
    findUpstreamFlowlines(seed, comidList)

    # spit out results into a list that is read into the SelectByLAyer command
    print comidList
    comidList1 = ",".join([str(k) for k in comidList])
    arcpy.MakeFeatureLayer_management("Q:/NHDPlusV2/NHDPlusMA/NHDPlus02/NHDPlusCatchment/catchment.shp", "lyrUp")
    where = "FEATUREID IN ({})".format(comidList1)
    arcpy.SelectLayerByAttribute_management("lyrUp", "", where)
    # join cat shape attributes to stats, write out
    arcpy.CalculateField_management(tv,"COMID", "{}".format(seed))
    arcpy.Dissolve_management("lyrUp", "B{}".format(seed), ["COMID"], "", "", "")

    print "Finished making basin for COMID {} ".format(seed)

    # Call to move on to the next row
    row =
  • No labels


  1. Unknown User (

    This is really helpful. I will try to get my WBD climbing algorithm similarly abstracted.

    BTW Here's how I have settled on laying out my GP tools to they can easily called from Python as a script tool:

    import arcpy
    def Nifty function():
        """Isolated nifty function"""
    def MyTool(...):
       """tool code"""
    if __name__ == "__main__":
        # ArcGIS Script tool interface: Get arguments, call tool
        argv = tuple(arcpy.GetParameterAsText(i)
            for i in range(arcpy.GetArgumentCount()))



    1. It would be interesting to flesh out the differences between using the built-in sys.argv[] and arcpy.GetParameterAsText(), and issues with either approach, perhaps on its own page.