• Spatial Analyst.
• TauDEM (Terrain Analysis using Digital Elevation Models) Toolbar.  This is software developed at USU for hydrologic digital elevation model analysis and watershed delineation and may be obtained from http://hydrology.neng.usu.edu/taudem/.  The component object model (COM) technology used in ArcGIS enables integration of third party COM objects into ArcGIS and this capability has been used in this tool.
• Extendability of ArcGIS using VBA programming

Exercise

• 0950107160101\00001\disk01\area01.  This contains the National Elevation Dataset data we will use.
• 16010203.  This contains the National Hydrography Dataset data we will use.

Now follow a similar series of steps to project the nhd coverage to the same coordinate system. The nhd coverage is '16010203\nhd'.  Save it in the same folder as 'loganutm\nhd'.  Leave the 'Save as type' set to coverages.  Click save, next and finish and wait until the coverage is projected.  Dismiss the execution warning error if you get it.  Close ArcToolbox.   For the remainder of this project you do not need an ArcInfo license.  An ArcView licence is sufficient.  [If you are working this exercise from an ArcView system the projected grid is available in the interchange file ned.e00.  ArcView/ArcMap has the capability to project data on the fly, so do not worry about not having the NHD streams in UTM coordinates.  Once they are added to a data frame in the correct projected coordinates they will display correctly.]

Open ArcMap.  Add the data 'loganutm\ned' and 'loganutm\nhd\route.rch' [or the unprojected '16010203\nhd'].  While loading the dataset 'ned' if prompted click OK to build pyramids.  Adjust the symbology to your taste, e.g. blue rivers and a nice shading for the elevation dataset.  Dismiss the warning about the elevation dataset having more than 2048 unique values and therefore no raster table.  We will not need a raster table.  Examine the source tab under layer properties associated with 'ned' to see the number or rows and columns and cellsize associated with this dataset.

Contours are a useful way to visualize topography.  Spatial analyst allows us to create contours.  Before using Spatial Analyst you need to make sure that it is selected.  Do this in ArcMAP by choosing Tools/Extensions ....  Place a checkbox next to Spatial Analyst, then click close. Choose View/Toolbars/Spatial Analyst to display the Spatial Analyst toolbar.   The Spatial Analyst toolbar looks like this.

From the Spatial Analyst menu select Surface Analysis/Contours ...

Use a 50 m interval from the input surface 'ned'.  Browse to set the output features as 'loganutm\ned50m.shp' contours for the ned dataset.  The dialog should look like

Note that the contour dialog gives the data height range (Zmin and Zmax).  Note the lowest and highest elevations values in this grid.  Lets explore the data a bit.  Zoom in on an area of interest and use the identify button  to examine a few grid values.   Logan Canyon is the canyon that drains north to south then turns to the west and leaves the map.  Logan Peak is the mountain to the south of this canyon.  Zoom in on Logan Peak and click on the highest point and note the height of Logan Peak.   Note that the elvation values you have been reading are decimal numbers.  The National Elevation Dataset records elevations using floating point numbers.  You should note that Logan Peak is not quite as high as the highest point on the map relative to the Zmax that you noted. Spatial analyist also provides the capability to query data so we will use this capability to find the highest point on the map.  Select "Spatial Analyst/Raster Calculator ...".

This can be used to assemble arithmetic like operations to perform on different rasters in the map.  Here double click on the layer 'ned' then click on the '>' then enter a number just less than the maximum elevation to assemble the expression '[ned] > 3033' shown.  Click on Evaluate. A new layer named 'Calculation' appears on your map.

In this case the entire map except for three grid cells is the 0 value representing false.  The three grid cells with value 1 representing true have elevation higher than 3033 m.  They are Mount Elmer.  See if you can find them.  The coordinates are 443972, 4640154 if you are having trouble.

The data you have been looking at has was converted at the beginning into UTM coordinates.  ArcMap has the capability to project data on the fly.  Click the add layer button and add the original NED DEM grid '0950107160101\00001\disk01\area01\demgrid' in geographic coordinates.  This should be displayed on the map aligned with the UTM projected grid.  Change the symbology on the geographic coordinates 'demgrid' layer to a classified, rather than stretched color ramp.  This lets you see the outlines of single grid cells.  Note that the projected grid cells appear rectangular.  Square grid cells in geographic coordinates are rectangular when projected into UTM.  You will also notice that the elevation values are slightly different due to the interpolation during projection.

Hillshading provides another nice visualization of topography.  Select Spatial Analyst/Surface Analysis/Hillshade ... Select the input surface 'ned'.  The remaining parameters in the dialog box control the illumination angle and position and vertical exageration.  I like to set the z factor to 10 for a dramatic effect.  Leave the other parameters at their defaults.  Click OK.  You should see an illuminated hillshaded view of the topography.

To turn in.  A layout with a depiction of topography either with contours or hillshade in nice colors.  Include the streams from NHD.  Label Logan Peak and Mount Elmer and indicate their elevation.

Open ArcMAP and add the TauDEM toolbar.  [Click on tools | Customize | Add from file and select the file c:\program files\Taudem\agtaudem.dll]
You should get a toolbar that looks like

This may be docked.

On the TauDEM toolbar from the Grid Analysis menu 'Select Base DEM grid ...'.  You should get a dialog like this:

Adjust the Base DEM Layer to be 'ned' and click OK.  Next, from the Grid Analysis menu select 'Fill pits'

Check at the dialog that opens that the input grid is 'ned' an the output one 'nedfel'.  Ignore the Flow Path Grid items for now.

Click on Compute and wait for a few minutes (or more) for the job to complete.

Each of the remaining command under the grid analysis menu can be run in sequence from top to bottom.  The program will provide suggested file names for the outputs that follow the convention described in the documentation and given below so you can just click compute at each dialog and have the command execute.  [You may learn a bit more about each file by clicking on the file box label in the dialog box associated with each command.  You can change the names of input and output grids if you like by editing in the textboxes or using the Browse buttons.  If you get an 'Error 1' it means that one of the input grid names specified has a problem.  This is commonly caused if commands are run out of sequence.]

Here we will shortcut running each command by selecting the 'Do all preprocessing' command from the 'Grid Processing' menu to have all results generated without any prompts or interuption.  Processing of this will take 5 to 10 minutes for the Logan DEM.  You should answer OK to the prompt to overwrite the existing 'nedfel' file that you created above or else the automatic processing will stop.  A number of output layers will be added to the map.  The name suffixes designate the contents according to the file naming convention in the documentation  and given below.   Examine these grid layers to understand their contents.

To proceed further and delineate watersheds a shapefile containing outlet points needs to be created.  Open ArcCatalog.  Right Click on the folder 'loganutm' where you are working and select 'New/Shapefile...'.  Set the name 'outlet' and set the feature class to point.  Click 'Edit...' to change the coordinate system then 'Import...' and select the 'ned' dataset in the 'loganutm' folder so that the coordinate system is inherited from this.  Click Add, then OK twice.  Add the shapefile 'outlet' to ArcMap.  It has no data yet.  Display the editor toolbar (View/Toolbars/Editor) and select Editor/Start Editing.  Select the folder that contains the shapefile 'outlet.shp', and set this as the target layer.  Set the Editor task to 'create new feature' and use the create new feature button  to carefully locate a point at the outlet of Logan Canyon near the SW corner of the domain.

Use the *src layer as a backdrop to ensure that you are locating a point on a valid stream path. Select Editor/Stop Editing and Save edits.  This is now a one point shapefile.  More points for multiple channel networks can be added if desired.  If you are having difficulty creating a shape file with valid outlet point you may download outlet.shp and use it.

Now select TauDEM/Network Delineation/Select Outlets Shapefile ...  Browse to select the shapefile and click Open and OK. Now select TauDEM/Network Delineation/Do All Network and Watershed Delineation steps.  A channel network and watersheds are delineated the TauDEM default settings that automatically does a constant drop test using an upwards curved drainage area threshold to delineate channels.   The channel network shapefile 'nednet.shp' has an attribute table that includes among other fields the stream order.  This can be used in symbology to indicate stream order by different colors or line thickness.

Prepare a layout map of the Logan River showing the TauDEM delineated channel network and watersheds draining to each stream segment.  On this layout indicate stream orders in different colors or symbols.

The command 'TauDEM/Network Delineation/River Network Raster Upstream of Outlets' enters a dialog that allows you to select different channel network delineation algorithms and their parameters.

Uncheck the upstream of outlets box and click Compute.  Then run the commands 'Network Delineation/Stream order grid and network files', 'Network Delineation/Stream shapefile and watershed grid', and 'Network Delineation/Watershed grid to shapefile'.  Click OK to overwrite existing files (or rename the outputs).  Output will be produced for the entire domain.  Overlay the NHD stream network with streams delineated using the above procedure and examine discrepancies that occur at the north end of this watershed.

Prepare a layout that shows TauDEM delineated streams and watersheds over the domain compared to NHD streams. Use an inset on this layout to identify the problem that occurs at the north end of this watershed where a subwatershed that according to the NHD drains in to the Logan River is in the NED DEM derived channel network diverted to the north.

The command 'Network Delineation/River Network Raster Upstream of Outlets' enters a dialog that allows you to select different channel network delineation algorithms and their parameters.  The automatic threshold selection by drop analysis is only possible when outlets are used.  After the command has been run the accumulation threshold parameter selected by the constant drop procedure is updated.  Report the threshold selected by the constant drop procedure with the DEM curvature based method.

Rerun the command 'Network Delineation/River Network Raster Upstream of Outlets' but select the "Contributing Area Threshold" method and set the constant drop search to be between 50 and 5000.  Report the contributing area threshold selected in this case.

The command 'Network Delineation/Constant Drop Analysis' shows details of the constant drop analysis that can be used to objectively determine drainage density based upon topographic texture and that is invoked automatically by default controlled by a checkbox within the 'Compute Raster of Network Upstream of Outlets' dialog.  The 'Constant Drop Analysis' also reports the drainage density associated with different channel delineation thresholds. Report the drainage density identified from the constant drop procedure with the curvature based and contributing area methods.

Recall that drainage density is Total Channel Length/Drainage area.  Use GIS methods you have learned to determine and report the total channel length and drainage area associated for each of the curvature based channel network and verify the drainage density given by the constant drop analysis.

Select the portion of the NHD river network that drains the Logan River watershed that you delineated.  Determine the total channel length of NHD streams within the Logan River watershed and the drainage density associated with the NHD river network.

It is not always possible to do what is needed using the functionality provided by ArcGIS.  Therefore the capability to extend ArcGIS through programming or scripting capability has been developed using Visual Basic for Applications (VBA).  This is similar to the concept of Macro's in Excel.  Here you will implement a simple VBA Macro to calculate the flow path distance from each point in the watershed to the outlet.  [This is something that can actually be achieved using the flowpath command in ArcInfo command line GRID functionality, and perhaps through the spatial analyst object library but programming it is a good learning exercise.]

From within ArcMAP select Tools/Macros/Macros...

Enter a Macro name: dist and click Create.

Cut and paste the following VBA code into the Visual Basic Editor that opens.  Be careful to paste code above, inbetween and below the stub code

Sub dist()

End Sub

Sub dist()     ' This is the subroutine that will do most of the work.
' Dimensioning variables
    Dim Foldername As String
    Dim OutletShapefilename As String
    Dim pGridname As String, dGridname As String
    Dim i As Long, j As Long, k As Long

' Specifying inputs.  YOU WILL NEED TO CHANGE THE FOLDERNAME BELOW TO REFLECT THE LOCATION OF YOUR WORK
    Foldername = "E:\Academics\GIS_in_WR\session19_dem_hydro\spawnex"
' YOU WILL NEED TO CHANGE THE FILENAMES BELOW TO REFLECT THE NAMES YOU ARE USING
    OutletShapefilename = "outlet"  ' THIS IS THE SHAPEFILE CONTAINING THE OUTLETS
    pGridname = "spawnp"             ' THIS IS THE D8 FLOW DIRECTIONS FILE CREATED BY THE COMMAND 'D8 FLOW DIRECTIONS'
    dGridname = "spawndout"        ' THIS IS THE NAME OF THE GRID FILE THAT WILL BE OUTPUT

'  General ArcGIS objects that need to be initiated
    Dim pDoc As IMxDocument
    Set pDoc = ThisDocument

' Create the workspace factory and open shapefile feature class
   Dim pWSF As IWorkspaceFactory    ' This says that the name pWSF will refer to a WorkspaceFactory Object
   Dim pFWS As IFeatureWorkspace    ' This says that the name pFWS will refer to a Feature Workspace (in this case a folder containing feature datasets)
   Set pWSF = New ShapefileWorkspaceFactory    ' This instantiates (creates a new instance of) a WorkspaceFactory object
   Set pFWS = pWSF.OpenFromFile(Foldername, 0)   ' The function OpenFromFile that is part of (a method within) the Workspace Factory object is used to open the given folder as a Feature WorkSpace and assign the result to the named FeatureWorkspace object pFWS.
   Dim pFeClass As IFeatureClass    ' This says that the name pFeClass will refer to a Feature Class (in this case a shapefile).
  Set pFeClass = pFWS.OpenFeatureClass(OutletShapefilename) ' The function OpenFeatureClass is part of (a method within) the FeatureWorkspace object and here returns the feature class object (in this case shapefile) being worked on.
        ' OutletShapefilename is the shapefile name without suffix like "outlet"

'READ THE COORDINATES FROM THE OUTLET SHAPEFILE
' get the number of features
    Dim nFeCount As Integer
    nFeCount = pFeClass.FeatureCount(Nothing)
    Dim dbXX() As Double, dbYY() As Double  ' X and Y coordinates will be stored in double precision arrays
    ReDim dbXX(nFeCount)
    ReDim dbYY(nFeCount)

    Dim pFeature As IFeature ' A feature object variable name
    Dim pPoint As IPoint ' a point object variable name
'Loop through to get the x and y coordinates of each point associated with each feature in the feature class
    For i = 0 To nFeCount - 1
        Set pFeature = pFeClass.GetFeature(i)
        Set pPoint = pFeature.Shape 'if it is a point shapefile
        'to get the coordinate of a point
        dbXX(i) = pPoint.X
        dbYY(i) = pPoint.Y
        MsgBox "Outlet " & dbXX(i) & " " & dbYY(i)   ' This pops up a message box giving the coordinates of each outlet
    Next i

' READ THE FLOW DIRECTION GRID
' Now work with the flow direction grid

'  Raster workspace factory object
'  Dim pWSF As IWorkspaceFactory    '  This statement would generally be needed, but is not here because we are reusing the same workspace factory name used above.
    Set pWSF = New RasterWorkspaceFactory

' Raster workspace object
    Dim pRWS As IRasterWorkspace
    Set pRWS = pWSF.OpenFromFile(Foldername, 0)

' Raster dataset object
    Dim pRsDS As IRasterDataset
    Set pRsDS = pRWS.OpenRasterDataset(pGridname)

' Raster object
   Dim pRs As IRaster
   Set pRs = pRsDS.CreateDefaultRaster

' Raster band collection object.  Rasters may contain multiple bands
   Dim pRsBC As IRasterBandCollection
   Set pRsBC = pRs

' Raster band object
    Dim pRsBand As IRasterBand
    Set pRsBand = pRsBC.Item(0)    ' The first band, i.e. item number 0 in the band collection

' Raster properties object
    Dim pRsProp As IRasterProps
    Set pRsProp = pRsBand

'  Get the size and extent of flow direction grid
    Dim xll As Double
    Dim yll As Double
    Dim csize As Double
    xll = pRsProp.Extent.XMin
    yll = pRsProp.Extent.YMin
    ncol = pRsProp.Width
    nrow = pRsProp.Height
    csize = pRsProp.MeanCellSize.X

' CREATE DISTANCE TO OUTLET GRID WITH NO DATA VALUES
' Spatial reference object for new grid
    Dim pOutSR As ISpatialReference   ' Get the spatial reference of the input grid
    Set pOutSR = pRsProp.SpatialReference

'  Raster workspace object that has method to create new grid
    Dim pRWS2 As IRasterWorkspace2
    Set pRWS2 = pWSF.OpenFromFile(Foldername, 0)

' Origin object for new grid
    Dim nOrigin As IPoint
    Set nOrigin = New Point
    nOrigin.PutCoords xll, yll

' Create new grid as raster dataset object
    Dim pOutDs As IRasterDataset
    Set pOutDs = pRWS2.CreateRasterDataset(dGridname, "GRID", nOrigin, _
    ncol, nrow, csize, csize, 1, PT_FLOAT, pOutSR, True)

'  new raster indicated by "d" for distance at the beginning of the name
    Dim dRs As IRaster
    Set dRs = pOutDs.CreateDefaultRaster

' new band collection
    Dim dRsBC As IRasterBandCollection
    Set dRsBC = dRs

' new band
   Dim dRsBand As IRasterBand
    Set dRsBand = dRsBC.Item(0)

' new properties
    Dim dRsProp As IRasterProps
    Set dRsProp = dRsBand

' Raw Pixel objects.  These are used to access the raw data in grids
    Dim pRawPixels As IRawPixels, dRawPixels As IRawPixels
    Set pRawPixels = pRsBand
    Set dRawPixels = dRsBand

'  No data value
   Dim noDataValue As Double
   noDataValue = dRsProp.noDataValue

' Set up double point object specifying the extent of the part of the grid to work with.  In this case the entire grid.
    Dim pPnt As IPnt
    Set pPnt = New DblPnt
    pPnt.SetCoords pRsProp.Width, pRsProp.Height

' Pixel block objects to work with,within extent specified by pPnt (in this case the entire grid).
    Dim pPixelBlock As IPixelBlock, dPixelBlock As IPixelBlock
    Set pPixelBlock = pRawPixels.CreatePixelBlock(pPnt)
    Set dPixelBlock = dRawPixels.CreatePixelBlock(pPnt)

' Origin of coordinates to be read from raw pixels into pixel block
  Dim pOrigin As IPnt
    Set pOrigin = New DblPnt
    pOrigin.SetCoords 0, 0

' Read block of flow direction pixels into PixelBlock
    pRawPixels.Read pOrigin, pPixelBlock

' Convert PixelBlock into arrays that can be referenced by rows and columns.  These are dimensioned before the sub statement so that they are global variables accessible to other modules
    pPixels = pPixelBlock.SafeArray(0)
    dPixels = dPixelBlock.SafeArray(0)

' Initialize the distance pixels to no data
    For i = 0 To nrow - 1
        For j = 0 To ncol - 1
           ' col, row
            dPixels(j, i) = noDataValue
        Next j
    Next i

'PREPARE ARRAYS FOR NAVIGATING GRID
'   Set up directions
'  row offsets
    di(1) = 0
    di(2) = -1
    di(3) = -1
    di(4) = -1
    di(5) = 0
    di(6) = 1
    di(7) = 1
    di(8) = 1
'  column offsets
    dj(1) = 1
    dj(2) = 1
    dj(3) = 0
    dj(4) = -1
    dj(5) = -1
    dj(6) = -1
    dj(7) = 0
    dj(8) = 1
'  distances
    For k = 1 To 8
        dd(k) = (((csize * di(k)) ^ 2 + (csize * dj(k)) ^ 2)) ^ 0.5   ' Pythagorous Theorem
    Next k

' CONVERT EACH OUTLET COORDINATE TO ROW AND COLUMN AND START DISTANCE CALCULATION AT IT
    For k = 0 To nFeCount - 1
       i = nrow - ((dbYY(k) - yll) / csize)
       j = (dbXX(k) - xll) / csize
       dPixels(j, i) = 0
       DistCalc i, j   ' This is a recursive subroutine call
    Next k
'  Save the result
    dRawPixels.Write pOrigin, dPixelBlock

'  Add the new raster layer to the document display
    Dim ROutLayer As IRasterLayer
    Set ROutLayer = New RasterLayer
    ROutLayer.CreateFromDataset pOutDs
    pDoc.AddLayer ROutLayer

End Sub

'RECURSIVE DISTANCE CALCULATION FUNCTION
Sub DistCalc(i, j)
Dim k As Integer, inb As Long, jnb As Long
For k = 1 To 8  ' for each neighbor
    inb = i + di(k)
    jnb = j + dj(k)
    If (inb >= 0 And inb < nrow And jnb >= 0 And jnb < ncol) Then   ' guard against out of domain
        If pPixels(jnb, inb) > 0 Then   ' guard against no data
            If (pPixels(jnb, inb) - 4 = k Or pPixels(jnb, inb) + 4 = k) Then
            ' Here we have a grid cell that drains back to the grid cell we are at
                dPixels(jnb, inb) = dPixels(j, i) + dd(k)
                DistCalc inb, jnb ' Call the function for that pixel
            End If
        End If
    End If
Next k
End Sub
' END OF CODE

Congratulations!  You have completed a GIS VBA program.  Now the only limit on what you can do is your creativity.

Prepare a layout that shows the distance to the outlet of the Logan River.  Report the longest distance to the outlet from within the Logan River Watershed.

The following default naming convention is suggested and used by the TauDEM software.  Any file names may be used with interactive input, but I suggest sticking to this convention to avoid confusion.  File names are:
nnnnsss
nnnn comprises the name of the dataset. Maximum length is operating system dependent.
sss comprises the suffix used to designate the data type as follows:

no suffix.	Elevation data.
fel	Pit filled elevation data.		produced by Fill pits and nondraining hollows
p	D8 drainage directions.		produced by D8 flow directions
sd8	D8 slopes.		produced by D8 flow directions
ad8	D8 contributing area’s, units are number of grid cells.		produced by D8 drainage area
slp	Dinf slopes.		produced by Dinf flow directions
ang	Dinf flow directions.		produced by Dinf flow directions
sca	Dinf contributing area, units are specific catchment area, i.e. number of grid cells times cell size.		produced by Dinf drainage area
plen	Longest path length to each grid point along D8 directions.		produced by Grid network order, Upslope total flow length, Upslope longest path length function
tlen	Total path length to each grid point along D8 directions.		produced by Grid network order, Upslope total flow length, Upslope longest path length function
gord	Strahler order for grid network defined from D8 flow directions.		produced by Grid network order, Upslope total flow length, Upslope longest path length function
src	Network mask based on channel source rules.		produced by RiverNetwork Raster
ord	Grid with Strahler order for mapped stream network.		produced by River Network Raster
w	Subbasins mapped using subbasinsetup.		produced Create Network and Sub-Watersheds
fdr	Flow directions enforced to follow the existing stream network		produced by Convert Connected Reach Network to Forced Flow Direction Grid
fdrn	Flow directions enforced to follow the existing stream network after cleaning to remove any loops		produced by flood when forced flow directions are used

Summary of answers to turn in

• A layout with a depiction of topography either with contours or hillshade in nice colors.  Include the streams from NHD.  Label Logan Peak and Mount Elmer and indicate their elevation.
• A layout map of the Logan River showing the TauDEM delineated channel network and watersheds draining to each stream segment.  On this layout indicate stream orders in different colors or symbols.
• A layout that shows TauDEM delineated streams and watersheds over the domain compared to NHD streams. Use an inset on this layout to identify the problem that occurs at the north end of this watershed where a subwatershed that according to the NHD drains in to the Logan River is in the NED DEM derived channel network diverted to the north.
• The threshold selected by the constant drop procedure with the DEM curvature based method.
• The threshold selected by the constant drop procedure with the contributing area method.
• The drainage density identified from the constant drop procedure for each of the curvature based and contributing area methods.
• The total channel length and drainage area associated with the curvature based channel network.  Verify the drainage density given by the constant drop analysis.
• The total channel length of NHD streams within the Logan River watershed and the drainage density associated with the NHD river network.
• Layout that shows the distance to the outlet of the Logan River.  Report the longest distance to the outlet from within the Logan River Watershed.