Delineation of the watershed for GIS routing tutorial

Setting up ArcGIS Pro to Create GIS Routing Model Input Files

Create the new project

Create a new, blank project in ArgGIS Pro using the blank Map template and called GIS Routing Tutorial.
Go to Insert -> Add Folder to add a connection to the folders containing data required for the GIS Routing Model tutorial
- The tutorial mask included in the GitHUB repository under the folder TrainingMaterials, and
- The more general datasets needed to setup any routing model in the United States, which can be found at /depot/phig/apps/GitHub/GIS-Routing-Model/VIC_Routing_Data.
Go to View -> Catalog Pane to open ArcCatalog in the right-hand panel.
Open the Maps tab.
Right-click the tab heading "Map"
- Select Rename
- Rename this Tab VicGridCellData
Go to Insert -> New Map -> New Map to add a second map.
Right-click on the new map in the catalog pane and rename it to HighResRoutingData.

Creating the mask file for the VIC model domain

Note

IgnoreSkip this section if using the tutorial.asc file as your mask file.

Make sure that you are working in the VicGridCellData map for these steps.
Input data: Load NLDAS 16th degree shapefile (CONUS) + HUC03 shapefile to define the boundary/extent of your study region shapefile.
Pairwise Clip: Clip NLDAS 16th deg to HUC 3 extent = NLDAS_HUC03
Extract by Polygon: Shapefile to raster -- NLDAS_HUC03 to raster file with no extension (confirm the grid spacing)
Raster to ASCII: Use the raster created in previous step to create huc_16th.asc (130701 to 751109) make sure CELL_ID is used.

Establish VIC model domain

This part of the processing is based on the VIC model grid cell resolution (e.g., 1/8\(^{th}\) or 1/16\(^{th}\) degree resolution, unprojected, latitude-longitude coordinates)

Make sure that you are working in the VicGridCellData map for these steps.
Adding the mask File and Renaming: Use the Copy Raster tool to copy the Spatial Model Tutorial cell number file (TutorialMask.asc) to a raster file (VicGridCellNumbers) in the project geodatabase.
If you are working on your [project] you need to create this TutorialMask.asc for your study region as outlines in the previous section.
The new raster should be added to the open project.(Note that there is a gap in cell coverage in the northeast of the model domain, this is an area with no surface outlet, so it would be filled as a sink if left in the domain).
Conditional Evaluation of each cell in the input Raster:
- Select the CON tool from spatial analyst (you will need to activate this extension license, if you have not used it previously on the current computer).
  
  Note
  
  Include instructions here...
- Select the grid cell number file (VicGridCellNumbers) as the input raster file.
- Set clause to be "VALUE > [0]"
- Set TRUE to 1, and leave FALSE empty.
- Output to a new raster in the project geodatabase (GIS Routing Tutorial.gdb) called VicGridCellMask.
- Run the operation.
- Should add a new raster mask with a single value (1) to the project.
Convert Raster to Polygon feature: Select the Raster to Polygon tool (Conversion Tools)
- Set the input raster to the mask file
- Set the field to Value
- Set the output to a file in the geodatabase called VicGridCellDomain
- Uncheck Simplify Polygons
- Run the tool
- This should add a new polygon to the project that outlines the contents of the raster file exactly (if simplify polygon was left check then the resulting polygon will not be an exact match -- delete and start step over).
Use this polygon to clip the USGS gauge station layer to the study domain (it is not projected -- though perhaps the file shared should be reprojected to match)
- Clipping USGS stations for the polygon: Use Clip (Analysis Tool) for the USGS gaging station data, called gagesll_9322_sept30_2011.shp. Input feature is shapefile of USGS gage data, Clip Feature is the VicGridCellDomain.
- Run and save the new file named VicGridGageLocations in the GISRoutingTutorial.gdb
Finally, convert the VIC model grid cell raster file into a polygon file that outlines each VIC grid cell -- this layer will later be reprojected to match the high-resolution routing layers and be used to conduct the zonal statistics required to quantify the high-resolution routing analysis into something that can be used for the coarser resolution VIC model.
- Select the Raster to Polygon tool from Conversion Tools.
- Set the input raster to [VicGridCellNumbers], which has a different cell number in each raster pixel.
- Set the Field to Value
- Save the output feature to the geodatabase GIS Routing Tutorial.gdb (the one created for this project) and name the output feature file VicGridCellNumbers_feature.
- Turn off Simplify polygons.
- Run the tool.
- Change the style to have no fill and a brighter outline color, and the results from this step will look something like the following screenshot where each VIC model grid cell is outlined by a polygon. The attributes for the polygon file match the VIC model grid cell numbers (now called [gridcode] instead of the default Values of the raster file).

High-resolution data analysis

This part of the processing is based on higher resolution, projected elevation and land use data. It must be projected so that the units of length and velocity needed to develop hydrographs are correctly assigned. Additionally, the data should be higher resolution than the VIC model setup as the GIS routing model setup will quantify the effects of path length, surface roughness and changing channel properties at the higher resolution in a way that can be reflected at the coarser VIC model resolution. At this time, this step typically uses 100 m or 30 m data as they are most common and well developed.

High Resolution Data Setup

This will start by reprojecting the domain polygon to the projection used by the DEM and Land use maps. Then data used for the next steps in the process is clipped to the model domain.

Make sure that you are working in the HighResRoutingData map for these steps.
In the catalog pane, right-click on Databases and select "New File Geodatabase".
- Name the new GDB HiRes GIS Routing Tutorial.gdb
- Save in the same project directory as the GIS Routing Tutorial.gdb folder.
Select Project Data Management Tools
- For the input, select the polygon domain file VicGridCellDomain
- Set the output to the new [HiRes] geodatabase and name it HiResRoutingDomain.
- Click on the world icon to the right of Output Coordinate System.
  - In the Coordinate System window, click on the world icon with a green plus to Add a Coordinate System, and select "Import a Coordinate System..."
  - Navigate to the VIC_Routing_Data/DEM folder and select the 100 meter resolution DEM file (USA_Contiguous_Albers_equal_Area_Conic).
  - Select OK
  - Select OK
- You will see a yellow triangle next to the box labeled Geographic Transformation because the two datasets are using different datums (NAD83 and WGS84), but you can ignore this warning for this process (note that NAD83 and WGS84 are typically within 1 meter of each other so ArcGIS has never provided a transformation algorithm).
- Select Run
- This should import the domain polygon with the USA_Contiguous_Albers_Equal_Area_Conic projection.
Repeat for the USGS gage dataset: VicGridGageLocations -> HiResGageLocations
Extract cells of raster corresponding to the area of the mask:
- Open the Extract by Mask from Spatial Analyst Tools
- Set the input raster to the [100 meter] DEM (DEM_100m.tif) provided in [VIC_Routing_Data]/DEM.
- Use the pull-down menu next to the "Input raster or feature mask data" box to select the [HiResRoutingDomain] polygon created in the previous step.
- Make sure that the Output raster is directed to the [HiRes] GDB and named "RedRiverDEM_100m".
- Click Run
- The clipped DEM will now be added to the HighResRoutingData map.
  - You may not be able to see it, however, as it will be drawn behind the polygon which by default is drawn as a filled object. You can either change the style of the domain polygon to be an outline with no fill, or turn it off to see the DEM.
  - In this figure, the domain polygon style has been changed to a thick back line with no fill.
Repeat the clipping process for the 100 meter NLCD landuse map NLCD_2011_100m.tif provided in the VIC_Routing_Data/Landuse folder to create the new layer RedRiverNLCD2011_100m.
Extract stream network to the clipped domain: Start the Pairwise Clip (Analysis Tools)
- Set the Input Features to the file USA_Rivers_Streams.shp in the VIC_Routing_Data folder.
- Set the Clip Features to the Hi-res domain polygon
- Set the output to be stored in the [HiRes] GDB and be named RedRiver_Rivers_Streams
- Click Run
- The stream network file should now be overlayed on the DEM and land use map, and be clipped to the VIC model simulation domain.
  
  Warning
  
  This layer is not in the same projection as the DEM, so you will need to reproject the line features from the source project to that used by the DEM file. Do this by starting the Project (Data Management Tools) tool.
- Set the input layer to RedRiver_Rivers_Streams
- Set the output to RedRiver_Rivers_Streams_forRouting
- Set the output projection by selecting the DEM (RedRiverDEM_100m) from the pull-down menu. This should set the output projection to USA_Contiguous_Albers_Equal_Area_Conic and select an appropriate transformation for the datum from WGS 1984 to NAD 1983.
Finally, you need to reproject the VIC model grid cell number polygon layer to match the projection of the high-resolution routing model analysis.
- Select the Project tool from Data Management Tools
- Set the Input Dataset or Feature Class to VicGridCellNumbers_feature (you will likely need to use the folder icon to navigate to the layer in the geodatabase, since the layer is only visible on the VicGridCellData map tab).
- Redirect the output to the [HiRes] GIS Routing Tutorial GDB file, and give it the name HiResVicGridCellNumbers.
- Use the pull down menu next to Output Coordinate System to select the DEM (RedRiverDEM_100m) and use its projection as the output coordinate system.
- Select Run to start the processing.
- This will result in a polygon feature where the projected VIC model grid cells are more rectangular than square (narrower east-west than north-south). This layer will be needed later for computing zonal statistics.

Watershed Delineation

Watershed delineation can be completed by installing the ArcHydro extension to ArcGIS and following the process described "at this link" - broken.

Install the correct version of ArcHydro for ArcGIS Pro version that you hav installed (if on ABE systems, contact Stan Harlow with the ArcGIS Pro version).
Once installed, activate the Add-In from the Add-In Manager under the Project tab in ArcGIS Pro
When activated, a new tab should appear at the top of the current project called Arc Hydro.
Open Toolboxes from Analysis Tab
Select Arc Hydro Pro Toolbox
Open Terrain Prepocessing Workflows
- Select Dendritic Terrain Processing with Imposed Drainage Line and Walls
- Set the Input Hydro DEM to RedRiverDEM_100m
- Set the Stream Feature Class to RedRiver_Rivers_Streams_forRouting
- Leave input Internal Wall Feature Class and Input External Wall polygon Feature Class blank as they will not be used for this processing.
- Leave other settings as their default values
- Click on Run -- it will take some time to complete (around 17 min)
  
  Note
  
  This seems to complete [all of] the Terrain Processing steps, and those layers then populate the existing [HighResRoutingData] project. Note that the tool's page will indicate that some of layers just created will be duplicated resulting in errors (and red x icons) on the tools page. This does not affect the success of the tool just [run, but] does highlight that the tool cannot be rerun without cleaning up the previous outputs.
- Output should look something like this:
Select Watershed outlet station
- Right-click on the layer HiResGageLocations and open the Attributes Table
- Click on the icon to Select By Attributes to open the selection window
- Make sure that the Input Rows is set to the HiResGageLocations layer and that the Selection Type is set to New selection
- In the query table, to the right of the Where statement, select STANAME, begins with, RED RIVER (type it).
- Click the Apply button, which should select three USGS stations on the main stem of the Red River.
- We only want one station for this demonstration, so now click on the Add Clause option.
- The new clause should be set as follows: And, STANAME, contains the text, FARGO.
- Click apply, and you should now have only the northernmost station still selected.
- Click OK
- Under the Feature Layer tab, on the Data tab, in the Selection group, click Layer from Selection.
- A new layer is added to the map or scene. The layer name is the name of the original layer, appended by the word selection (HiResGageLocations selection). Rename the layer to be RedRiverGage.
Delineate the watershed to the selected output [station]
- Return to the Arc Hydro Tools Pro panel and open Watershed Processing.
- Select the Watershed Delineation tool
- Set the Input Batch Point Feature Class to RedRiverGage
  
  Note
  
  Previously had errors for Input Stream Raster and Input Snap Stream Raster, which are both required for watershed delineation, but were apparently not created with the Terrain Preprocessing workflow selected. Setting these to the layer strlnk, which seems to be most relevant (includes only the steam location), helps.
  
  Note
  
  Zoom in and in this case the USGS gauge location is not quite on the stream channel, which caused the tool to only produce a single cell watershed (see below)
  
  Which is quite annoying, and means the the outlet point needs to be edited to have it overlap with the actual stream channel.
- To correct the problem, first we need to delete the Watershed and Watershed Point files that were created incorrectly.
- The Watershed and Watershed Point Layers are in the GDB file names after the project, so GIS Routing Tutorial Practice.gdb, not after the current Map. They are within a Folder called Layers.
  - Right-click on the layers Wateshed and WatershedPoint in the left panel, and remove each from the current map.
  - Click on the View tab and then the icon to open the Catalog Pane
  - Open the list of Databases
  - Open the GDB with the same name as the project file you are currently using (e.g., GDB Routing Tutorial Practice.gdb or HiResGISRoutingTutorial.gdb)
  - Open Layers then right click on the Watershed and WatershedPoint entries and delete each file.
  - Close out of the Catalog Pane
- Now we need to edit the outlet point file, so click on the Edit tab, and then on the Modify icon.
  - Click on Move
  - Click on the layer RedRiverGage and make sure that the outlet point is highlighted (is now a yellow circle on my system).
  - Use the mouse to select the outlet point (yellow dot).
  - Then left-click to grab it and drag it to the nearest part of the channel raster (blue square).
  Note
  
  For these sample images, only the layers RedRiverGage, HiResRiversStreams, strlnk and RedRiverDEM_100m are visible so the location of the stream channel is quite clear, as is the location of the streams line which overlays the stream channel raster file.
  - When you are happy with the modification, click on Save
- Close out of the editing window and return to the Watershed Delineation tool
- Set up the inputs to the tool again and Run.
- This will result in a watershed covering about ¼ of the total area in the southcentral part of the study domain (greenish area in the figure).
  
  Note
  
  You might need to change from the default fill color on the Watershed layer to see the results as ArcGIS Pro makes use of many muted and very similar colors by default, which can make it difficult to make out the new layer without turning other layers off or changing the fill color.
- The Red River of the North gage at Fargo, ND is unfortunately the further downstream USGS gage that incorporates only drainage areas within the United States. The domain used for this tutorial matches that for the VIC model spatial simuation tutorial which incorporates more of the Red River basin. For this tutorial, we want the output of the routing model to match with a monitored streamflow location for which we have simulation data for all grid cells.

Using the GIS Routing Model Toolbox

The next set of instructions will walk you through the process of applying the GIS Routing Model toolbox to the provided datasets to yield the output files required to run the GIS routing model using VIC model output files.

If you have not already, make sure that you can access the GIS Routing Model toolbox as follows:
- If you are a member of PHIG and are working on the RCAC cluster computers, there should be an updated version of the toolbox available on the PHIG depot drive under apps/ArcGIS_toolboxes/GISROutingModelTools/ArcGIS_toolbox/GIS Routing Model tools.tbx. The depot drive can be mounted on your [Windows] computer, and the toolbox can be imported directly from the networked directory.
- Otherwise, the toolbox cade can be cloned from the PHIGOrganization GitHUB repository
To add the toolbox to ArcGIS Pro, go to Insert -> Toolbox -> Add Toolbox.
- Navigate to the location of the GIS Routing Model toolbox using the window that opens.
- This will add the toolbox to the project space, so it will appear in the Catalog Pane under Toolboxes.

Routing model processing starts with Step 0 and runs through Step 8. Unless the instructions tell you otherwise, assume that all steps are conducted in the HighResRoutingData map. Work with the added toolbox GIS Routing Model Tools.tbx.

Step 1. Create RoutingLayers Folder

Set the box labeled ArcHydroAnalysis to the folder that contains the project file you built for this tutorial. That folder should contain the project folder, the two GDB files you created and a subfolder called Layers that was created by Arc Hydro. You can pick a different location, but it may not be auto populated correctly for subsequent tools.
You can leave the output Folder Name set to the default, RoutingLayers. This is where the intermediate files from the GIS routing model will be stored. Where to save the files to actual use with the routing model will be defined later.
Successful completion of this step should result in the creation of a folder called RoutingLayers in the project folder.

Step 2. Compute Shreve Stream Order

If you have followed the instructions so far, this menu should mostly self-populate and not leave you with any angry icons. The exception being that Arc Hydro no longer seems to keep the str layer, so now you need to set that to strlnk to continue.
- Run the tool.
- This should add a new layer called ShreveOrder which is a raster layer that will match locations with the stream channel but be populated with the stream order number for each reach of the river system. In this figure, it appears as the grey squares under the stream network feature.

Step 3. Calculate Slope of Land and River

This tool computes the slope of the hillslopes and the river channel.
It should be possible to complete this step simply using the defaults (or by changing str to strlnk).
- Run the tool, which will add the [LndRivSlope] to the current map view. Slopes are likely to be highest in the river channels (Grey lines in the figure).

Step 4. Estimate Hydraulic Radius

This tool estimates the hydraulic radius of the stream channel.
The Flow Accumulation and Routing File Path should be left as default values.
The method used to distribute hydraulic radius can be set as follows:
- T2 -- 2-year flood distribution, a reasonable default to get normal flows routed correctly
- T10 -- 10-year flood distribution, a reasonable default for getting normal high flows routed correctly
- T100 -- 100-year flood distribution, a reasonable default for getting extreme high flows routed correctly
- CUSTOM -- used when more is known about the stream channels being routed and more control over final routed streamflow is desired. requires the following additional information:
  - Minimum Flow Depth (Bankfull flow in tributaries [meters]) -- estimate of how deep bankfull flow is in the smaller tributaries within the watershed, must be provided in units of meters.
  - Minimum Drainage Area (used to define start of channels in ArcHydro [sq km]) -- is generally set to the same value used to start channels when delineating the watershed, must be provided in units of square kilometers.
  - Maximum Flow Depth (Bankfull flow at basin outlet [meters]) - estimate of how deep bankfull flow is at the watershed outlet, must be provided in units of meters.
  - Maximum Drainage Area (Drainage area of outlet gage [sq km]) -- is set to the watershed area at the outlet, must be in units of square kilometers.
  Note
  
  Bankfull depths can be estimated from representative USGS gaging stations that provide time series of stage height. And bankfull width from Google maps by measuring distance between banks of the river mostly bare of vegetation.
- For the tutorial use T10.
Once the hydraulic gradient distribution method is set, run the tool.
The resulting raster layer will have hydraulic radius estimated for all channels in the watershed. In the figure below, hydraulic radius ranges from 0.0109441 m to 41.513 m.

Step 5. Compute Manning N

This tool will compute the Manning N value for the hillslopes using the land use map, and the channels using the range of N values provided.
The Shreve Stream Order and Routing File Path should be set correctly by default.
The Land Use Map (Use to apply Manning N values to the land surface) should be selected for the land use data already clipped to the simulation domain.[ ]For the tutorial, you should have a layer called RedRiverNLCD2011_100m that was created earlier in the process.
As that land use file is based on the NLCD 2011 land use classification scheme, you will need to convert the land use class values to those employed by this tool. For the Landuse Remap Table, click on the folder icon to the right of the box and navigate to the location of the GIS Routing Model Toolbox (found previously). Within the ArcGIS_toolbox folder, select Databases -> Routing Table.gdb -> NLCD20112GISrouting and click OK.
The box for Manning N from Land Use Table (Maps Land Use Classes to Manning N values) should be set to the file Landuse2ManningN (if not this table is in the same GDB as the NLCD land use conversion table).
The Manning N value for the channel will be distributed along the channels in the watershed using the Shreve Stream Order. Values for Manning N will be interpolated from the Maximum Manning N (smallest channels) and the Minimum Manning N (largest channels). These values are adequate for most [applications, but] can be adjusted if desired to change the routed streamflow characteristics. Use the default values for the tutorial.
Once all settings are completed, Run the tool.
The resulting raster layer will be populated with Manning N values from 0.8 to 0.02 if you left the channel defaults (Manning N for the hillslopes is controlled by the values in the table not from the maximum and minimum values set for the channel). See how the roads and city affect the distribution of Manning N in and around the Fargo, ND gage.

Step 6. Calculate Flow Velocity

This tool uses the slope, hydraulic radius and Manning N values calculated in the previous steps to calculate the velocity of flow across the high-resolution grid cell.
This tool should be correctly populated with the outputs from the previous tools.
Maximum Velocity and Minimum Velocity in meters per seconds. The default values should work for most cases, but clearly adjusting those values will affect the timing of water arriving at the watershed outlet.
Click Run to apply the tool to your dataset.
Velocities in the final raster layer will range between the maximum and minimum values set for the tool.

Note

Velocities are highest on the slopes along and into the river channel and along roads and urban areas. Tributary channels with higher slopes will also have higher velocity flows.

Step 7. Compute Travel Time in Hours

This tool computes the travel time to the outlet for every grid cell in the domain. Note that this is computed on the full image, so the outlet is defined as where any channel runs out of the domain.
The Flow Direction and Flow Velocity layers should be set correctly.
Click on Run to apply the tool.
The longest travel times should be where the channel is furthest from the outlets of the clipped DEM. In this screenshot, travel time increases as you move further from the location of the main channel.

Step 8. Calculate Watershed Routing Statistics

This tool aggregates that high-resolution spatial data from the routing analysis into summary statistics for each coarser VIC model grid cell.[ ]This is the information that the routing model will [actually use] to build unit hydrographs for the outlet for each grid cell within the watershed. That in turn is used to work out when runoff and baseflow from each grid cell will reach the watershed outlet.
The Travel Time Raster and Watershed Boundaries should both be set correctly their default names from previous steps.
FieldID will default to HydroID (field name from Watershed), which is the default number given to watersheds in the batch delineation program. If you did not use Arc Hydro for watershed delineation, then make sure that the watershed feature file has an attribute column that has a unique integer identifier for each watershed being routed and set that as the FieldID.
Set the Simulation Zones to the projected polygon feature with VIC grid cell numbers, HiResVicGridCellNumbers.
Set the Zone Field to gridcode (field name from HiResVicGridCellNumbers) so that each Zonal field is annotated using the VIC model grid cell number.
The Fraction Threshold for Inclusion is used to exclude grid cells along the edge of a watershed that do not contribute substantially to the overall watershed discharge. The default value, 0.01 or 1%, will keep most grid cells in the routing model. Setting the number higher reduces the number of grid cells that need to be run while reducing the size of the simulated watershed. For this tutorial, leave the default value.
Click Run
This step will not produce a spatial layer, but instead will produce a Zonal statistics table for each watershed being setup for routing. These files will appear in the ZonalTables subdirectory inside the RoutingLayers folder.

Step 9. Export GIS Routing Model Setup Files

Run this tool from the VicGridCellData map panel, not from the HighResRoutingData map.
Add the ZonalFract and ZonalStat tables by using Add Data from the RoutingLayers/Zonal Tables folder, if the tables are not available on the VicGridCelldata map panel.
This tool will generate the required input files for the GIS routing model at the resolution of the VIC model simulation and output those files to the selected directory.
The layer Watershed Polygon (stored in the main GDB -> Layers) should have HydroID (field name from Watershed) and FieldID correctly populated using the default values, if you used Arc Hydro for watershed delineation. Otherwise, change as needed.
For the VIC Simulation Resolution Cell Number Raster, click on the folder and navigate to the project directory and the project GDB folder. Pick the raster file VicGriCellNumbers], which should be in the resolution of the VIC model.
For VIC Simulation Zones, select the layer VicGridCellNumbers_feature from the pull-down menu.
Zone Mapping Field should be correctly set to GRIDCODE.
Routing File Path should be set to the default.
To set the Path for Routing model Inputs, which is where the output of this tool will be written that is also the input files required to run the GIS routing model, click on the folder and navigate to the Project folder. Create a new folder called RoutingModelFiles at the same level as Layers from ArcHydro and RoutingLayers from this toolset.
If you are rerunning this tool to update the GIS routing model input files, then make sure to select Overwrite routing files. If you are processing a lot of watersheds (such as all USGS gaging [stations] locations for the INCCIA) then not selecting the overwrite feature will allow processing to skip over watersheds that were previously correctly processed.
Click on Run.
There will be no visible changes in the Project file, but if you navigate to the selected GIS routing model output folder, you will find five different types of output files (each output file is likely to have multiple file extensions under a common name). All filenames will start with "WS[=" and the ID number associated with the current watershed. Types of files produced include:
- WS*_area = Actual contributing area. This is the area in square kilometers contributing to discharge at the outlet for each VIC simulation cell.
  
  Note
  
  This will only be in square kilometers if the high-resolution watershed analysis was conducted in a projection using meters as the spatial unit. This should be updated in the future to confirm the output file units either by asking the user to select the high-res analysis units or by setting them after the user selects a representative input layer.
- WS*_fract = Fractional contributing area. This is the fractional area in square kilometers contributing to discharge at the outlet for each VIC simulation cell. Must run this through the script ConvertFractFile2Area.py in order to convert the values to square kilometers.
- WS*_mask = Watershed mask. Indicates which VIC simulation cells contribute to discharge for the current watershed.
- WS*_ttmean = Mean travel time. Mean travel time each VIC simulation cell to the outlet for the current watershed. Used in the next step to build the unit hydrographs for each simulation cell to the outlet.
- WS*_ttstd = Standard deviation of travel time. Standard deviation of travel time each VIC simulation cell to the outlet for the current watershed. Used in the next step to build the unit hydrographs for each simulation cell to the outlet.

Building Unit Hydrograph Files

After completing the ArcGIS analysis, the next step is to create unit hydrograph files for each watershed that is going to be routed. The unit hydrograph files is used by the routing model to assess how much total runoff (surface runoff + baseflow) makes it to the routed location each time step based on the travel times calculated using the ArcGIS parameterization steps.

The unit hydrograph files are currently computers using the original Fortran program, iuhcell.f90. DEvelopment has been started on a C-version of the code, MakeUnitHydrographFiles.c, that will function more like the other steps of the routing model. Check the GitHub documentation for updates on code development.

To generate unit hydrograph files from the data layers produced in the previous steps, complete the following steps:

Move the RoutingModelFile outputs to a working repository on scratch/working location to access the scripts/commands.
The command line options for the Fortran code are
```
Usage: iuhcell <mask grid> <mean time grid> <std time grid> <out file> <routing time step ([hrs])>
```
where
- <mask grid> = ws*_mask.txt
- <mean time grid> = ws*_ttm.txt
- <std time grid> = ws*_ttsd.txt
- <out file> = uh_cell.txt
- <routing time step (hrs)> = the time step in hours of the VIC model output (e.g., fluxes) file. Use 24 for dailyh, 3 for 3 hour and 1 for hourly.
The <outfile> format can be .txt (example: uh_cell.txt)

For PHIG members, the code is installed at /depot/phig/apps/linux/iuhcell.

Note

If you are logged into the current PHIG cluster and do not get a usage message when you try using the command iuhcell, you may need to add the path /depot/[phig]/apps/linux to your environment in your personal ~/.bash_profile.

The most up to date PHIG bash profile can be found in the PHIGOrganization GitHub page. Need to provide updated link to repository.

Specific to this tutorial, you need these modifications to the default .bash_profile file:

# define group basepath 
EXECDIR="/depot/phig/apps/" 

# define default application search paths
PATH=$PATH:.:${EXECDIR}/script/VIC:${EXECDIR}/script/public:${EXECDIR}/VIC:${EXECDIR}/linux:${EXECDIR}/script:${EXECDIR}/VICpreprocessing

Run the GIS Routing Model

Once all of the input files have been created, those from ArcGIS Pro that define the watershed, and the creation of the unit hydrograph files, then the GIS routing model is ready for application to any VIC model simulation of the watershed.

Note

Make sure that OUT_RUNFF and OUT_BASEFLOW are defined in a VIC model output file. That output file should be at the time step you want to use for the routing model. Thus if you want to route daily discharge, you need a VIC model output file with runoff and baseflow values for each day (24 hours). To route hourly streamflow, your VIC model output must be written as an hourly time step. Note that the VIC model can produce output files at different time steps, so you could produce one output file with hourly values for routing and another at a daily time step for other analysis. The VIC model can aggregate from finer to coarser time scales (hourly to daily), but to get hourly output the model has to be run using an hourly time step.

The unit hydrograph files must also be generated for the same time step.

The routing model is called GisStreamflowRouting and is installed on the cluster computers for PHIG members. Otherwise, the source code is included in the GIS routing model GitHub repository.

Usage: GisStreamflowRouting <flux path> <discharge time step> <mask grid> <unit hydrograph file> <fract file> <out file>

    GIS routing model will generate streamflow discharge for the
    selected watershed using output from the current VIC model
    simulation.

    Routing model parameter files must be created using
    ArcGIS toolbox before the model is run, and must match the time
    step of the VIC model flux files.

    <flux path> is the path where the flux files [tobe] routed are stored;
    <discharge time step\> is the time step of the output file (should
        match input parameters);
    <mask grid> is the Arcraster file with watershed mask;
    <unit hydrograph file> is the unit hydrograph for each grid cell;
    <fract file> is the Arc raster file with fractional contribution from
        each grid cell; and
    <out file> is the output discharge file.

Note

The represents area contributed (sq km) to routed discharge by each VIC model grid cell. This information is available in the ws*_area.txt file.

But GisStreamflowRouting needs the in the format shown below (first two columns represent CellID and VicID (as these values are not currently used that can be random numbers for this use case). The 3rd column is the name of the flux file that will be searched along the directory that has VIC flux outputs and 4th column is derived from ws*_area.txt which contains the contributing area.

1    1    fluxes_38.65625_-86.21875    33
2    2    fluxes_38.65625_-86.15625    220
3    3    fluxes_38.65625_-86.09375    33
...

Update section to reflect the script nbeeded to build this file!

arcgrid_to_xyz.script ws*_area.txt | awk '( "%i\t%i\tfluxes_%.5f_%.5f\t%f\n", NR, NR, $2, $1, $3 ) }' > fractArea_ws*.txt

Requirement: make sure arcgrid_to_xyz.script added to your bash profile.

Or

Use a python script to create this file

(arcgrid_to_xyz.script ws_area.txt > ws_Area.csv) - gives lat, lon and area
Create a dataframe two add two random columns and third column has 'fluxes_'+lon+'_'+lat and last column is the area. Make sure no headers in this file.

You will obtain a file with date(YYYY-mm-dd, same as the VIC flux output, as first column and streamflow (cms if your GIS work was in the unit of meters). Output will look something like below:

VIC_Routing_output$ head out_flow.csv

01-01-2005    2596.814589
01-02-2005    2947.181141
01-03-2005    2945.052011
01-04-2005    2924.321271
01-05-2005    2907.485796
01-06-2005    2894.613053
...

Now plot them for desired time period.

Using the Red River should mean that the output from the VIC Spatial tutorial can be used as input to the routing model. need to check this

References to check: https://github.itap.purdue.edu/PHIGOrganization/GIS-Routing-Model/wiki/Using-the-GIS-Routing-Model

It can also be completed by using the ArcGIS hydrologic analysis tools directly, as described in the following steps:

Make sure that you are working in the HighResRoutingData map for these steps.
DEM Reconditioning

Restarting Arc Hydro

ArcHydro is a very useful tool for working with hydrologic datasets, but getting things done sometimes requires stepping back in the processing to correct problems with input datasets. The following steps must be done within ArcGIS and the ArcCatalog tool. Deleting files in file explorer will not correctly document the change in the spatial database system and will cause problems within ArcGIS.

In the Content pane on the left side of the ArcGIS Pro map panel, select the layers you need to regenerate
- Arc Hydro raster files are typically lowercase names (e.g., agreedem, str, fdr, fdrstr) while feature files are merged names with first letters capitalized (e.g., DrainageLine).
- Use - to select multiple layers.
Note

If you were successful in running Arc Hydro previously and are not restarting because Terrain Preprocessing failed, you will also likely have some attribute tables associated with the stream and catchment networks that will appear at the very bottom of the Contents list (they will not be drawn directly, but instead are related to other feature files that Arc Hydro produces) – these files also need to be deleted.
Right-click on one of the selected layers and select Remove to remove all selected layers from the current project.

Next, you need to delete the source files or Arc Hydro will not be able to generate new ones using its default naming convention.

Open the Catalog Pane (View -> Catalog Pane)
Open the Databases folder
Open the database with the same name as the current project
Open the Layers folder (if there is no Layers folder, this is not the GDB file that Arc Hydro has been using)
Select the features you need to regenerate or all of them to start over.
Right-click and select Delete
Now look in the main folder of the geodatabase, select and delete files generated by Arc Hydro (e.g., Catchment_FS, DrainageLine_FS) – DO NOT delete any files you created outside of Arc Hydro
Finally, to delete the raster files produced by Arc Hydro, open Folders in the Catalog Pane.
Select the Folder being used for the hydrologic processing, the one where your project and geodatabase files are stored.
Open the Layers folder – if there is no Layers folder, then you are not in the folder being used by Arc Hydro
Select and delete the raster files you want to regenerate or all raster files to restart processing. Note that the “info” folder should not be modified, it will be updated by ArcGIS as it deletes the selected raster files.

Warning

Check the following feature(s) for potential problems in DrainageLine:

OBJECTID=172,469,1152,1163,1180,1189,1244,1390,1388,1553,1275,1416,1371,1486