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 myusgs@usgs.gov. Thank you for your prompt attention to this matter.
Skip to end of metadata
Go to start of metadata

How to create a drainage point at the most downstream "point" in the catchment. This "point" is being derived with raster data, so it's really the center of a grid cell with the largest value in the flow accumulation grid for that catchment.

The need is to tie catchment characteristic data to the NLDI and give a nice representative point for the data so storing the catchment boundaries with the characteristics won’t be needed.

  • ArcHydro tool
  • Con(zonalmax(cat, fac) == fac, cat) would give you the cells for each catchment.

Tried this out as a test, results are in attachment outlets.zip.

Ran sample code below and got these results ( Green cells are outlets from FAC):

Ran sample code below on all Production Units and they are posted here.


# Name: _Ex_02.py
# Description: Identifies the rate of maximum change
#              in z-value from each cell.
# Requirements: Spatial Analyst Extension

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *
import sys, os, traceback, arcpy
from arcpy import env
from egis import GPMsg, MsgError
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
# Set the workspace environment
wks = r"H:/Outlets/PUs"
env.workspace = wks
arcpy.env.scratchWorkspace = r'H:/Outlets/PUs'
env.overwriteOutput = True
env.qualifiedFieldNames = False
tv = "tv"
tvStat = "tvStat"
tmpTable2 = "temp2.dbf"
try:
        # Set local parameters
        list=['a','b','c','d']
        for id in list:
            # Set the workspace environment
            GPMsg("t","Starting script for "+id)
            env.workspace = wks
            env.overwriteOutput = True
            env.qualifiedFieldNames = False
        # Set local variables
            valGrid = "H:/Outlets/fac/NHDPlusV21_TX_12_12"+id+"_FdrFac_01/NHDPlusTX/NHDPlus12/NHDPlusFdrFac12"+id+"/fac"
            catPath = r"I:/NHDPlus2/PU12/NHDPlusCatchment"
            catGrid = os.path.join(catPath,"cat")
            catShape= os.path.join(catPath,"catchment.shp")
            outName = os.path.basename(valGrid)
            arcpy.env.cellSize = catGrid
            arcpy.env.Extent = valGrid
            arcpy.env.SnapRaster = valGrid
            arcpy.env.mask = valGrid
            outlets = r"H:/Outlets/PUs/out12"+id

            # run raster stats
            GPMsg("t","starting raster stats for region " + id)
# Set local variables
            inZoneData = catGrid
            zoneField = "Value"
            inValueRaster = valGrid
            arcpy.env.cellSize = "I:/NHDPlus2/PU12/NHDPlusCatchment/cat"

# Check out the ArcGIS Spatial Analyst extension license
            arcpy.CheckOutExtension("Spatial")
            GPMsg("t","starting")
# Execute ZonalStatisticsAsTable
            outZSaT = Con((ZonalStatistics(inZoneData, zoneField, inValueRaster, "MAXIMUM") == inValueRaster), inZoneData)
            GPMsg()
            print 'ZonalStats COMPLETE.'
            GPMsg("t","Calculating")
            print 'ZonalStats COMPLETE.'
            GPMsg("t","Processing output...")
            # Set local variables
            inRaster = outZSaT
            outPoint = "H:/Outlets/PUs/outlet12"+id+".shp"
            field = "VALUE"

            # Execute RasterToPoint
            arcpy.RasterToPoint_conversion(inRaster, outPoint, field)

        # Set local variables
            arcpy.MakeFeatureLayer_management ( outPoint, outlets)
            inFeatures = outPoint
            layerName = outlets
            joinField1 = "GRID_CODE"
            joinTable = "I:/NHDPlus2/PU12/NHDPlusCatchment/featureidgridcode.dbf"
            joinField2 = "GRIDCODE"
            arcpy.JoinField_management (inFeatures, joinField1, joinTable, joinField2)
            arcpy.AddXY_management(outPoint)

except MsgError, xmsg:
        GPMsg("e",str(xmsg))
except arcpy.ExecuteError:
        tbinfo = traceback.format_tb(sys.exc_info()[2])[0]
        GPMsg("e",tbinfo.strip())
        numMsg = arcpy.GetMessageCount()
        for i in range(0, numMsg):
                GPMsg("return",i)
except Exception, xmsg:
        tbinfo = traceback.format_tb(sys.exc_info()[2])[0]
        GPMsg("e",tbinfo + str(xmsg))


Coast Line Curiosities

Coast line area along Delaware Bay.

 

Blue lines denote flow lines "With Digitized" which means they will have catchments associated with them. \

Notice the coast line is included in this set.

Flow lines with catchments (in brown).

 

Black dots denote outlets generated from the script posted above. Notice the entire coast line in some instances have multiple outlets. How should we deal wit these instances or do we need to?

CATSEED grid in purple and green.

 

SPARROW has a file of downstream nodes on flow lines but they don't always coincide with the fac. Blue is CATSEED, black dots are fac-derived outlets, red dots, end nodes of flow lines.


  • No labels

19 Comments

    1. what about catchments draining to a sink? Point might not be on perimeter. Is this okay?
    2. what about catchments with multiple outflows, such as those draining to ocean or other waterbodies? okay to just take the zonalmax(demfac)?
    1. Good points. This is the caveats and weirdness of NHDPlus. QAQC would definitely have to be done.

    2. IMHO, in both these cases, since the single point is more or less arbitrary and clearly not physically representative. The lowest point in the catchment is as good as any. As long as the point is IN the catchment, I'm happy.

      1. So I think Dave's point about what he needs should be in the problem statement. This product does not necessarily need to be valid from a more rigorous analytical perspective. 

        Dave, what about determination of the downstream catchment? Is this part of the problem definition? 

  1. optimization question: if the flowlines were burned into the DEM, I'm presuming that a raster representation of the flowlines exists. If so, might use this for zonalmax() and con() since it'll have many fewer cells but include the desired point.

    1. This ran pretty fast as is! But this is something to consider. The "burning in" is for catchment delineation only and not for stream generation, so the question becomes would the "burned in" raster capture the right FAC cell?

      1. I disagree that the burned-DEM is only for stream generation. I think it's part of the coordinated set of products that defines the "NHDPlus" way. The burned-DEM is "reality" if you've signed on to use that data set, isn't it?

    2. Unknown User (cprice@usgs.gov)

      In HR with have a an "swnet" raster to use as a mask in use cases like this.  But, that said, this zonal calculation is pretty fast. Maybe not worth gilding the lily.

      I was thinking one could use the flow direction raster to "look downstream" like I did with the WBD Pour Point tool to identify downstream COMID (and sinkpoints) There's no ambiguity here as the catchments and flow derivatives are in synch!). However, I bet this information could be very quickly located using the NHDPlus flow table, especially if indexes are built.

       

      1. I like the idea of using fdr to look downstream from whatever point/cell gets selected.

        The point about index-building is interesting, but is it part of Dave's asked for right now? Regardless, seems like this kind of information would be known at some point of the NHDPlus creation process and why not just write it out?

      2. Unknown User (cprice@usgs.gov)

        It occurred to me that using catseed may work as well. Fewer cells to look at that way.

        1. what's catseed? (I've heard of bird seed...)

        2. Catseed is the "rasterized  flowlines and the only VALUE is the GRIDCODE of the catchment. This was used to generate catchments.  The entire flowline is the "seed" for catchment generation and it does not represent the outlet. I think the way to go is use the FAC and modify the script here to include COMID and XY, an easy addition.

          1. How to deal with coast lines can be a problem. CATSEED may come in handy but in some instances the CATSEED grid is the coast line.

  2. What about just making a point 10% of the way up the flowline vector? Going w/out any raster would probably be really fast!

    1. Unknown User (cprice@usgs.gov)

      I think you'd likely get into trouble where the vector flow and the raster flow don't exactly line up due to discretization effects. The gridcell flow is perfectly in synch with the catchments that are derived from it.

      1. wouldn't the burning process preclude this? I thought that was exactly what it was for.

        1. Unknown User (cprice@usgs.gov)

          The flow lines are not modified at all, the raster flow is derived from it at resampled resolution. So no, the burning process sometimes leads to some weird stuff. What we call "curiosities" in the NHDPlus training course.

          1. Huh. Guessing the burn buffer is smaller than I would have guessed. Seems like this could be improved somehow to assess local topography and compensate as necessary--but that's probably way out of scope.

            Back to Dave's request, would this kind of weird stuff still be a problem as long as we got 1) a decent visual representation for cartographic presentation, and 2) accurate attribute value of downstream catchment identity?

            1. I'm not worried about networking with this. We have the node table for that, right? I'm looking for a lat/lon that is IN the catchment and a reasonable location to call the catchment's outlet.