Preprocessing of Routing Model Inputfiles

Introduction

This document is provided as part of the complete documentation to the Variable Infiltration Capacity (VIC) model and the Routing model. It presents a short description of how to derive all the necessary input files for the routing model. From this site you will also find information and links to datasets and software that you will need to build all the input files. Some knowledge of basic ArcInfo commands and operations is required, but hints as to which commands to use are provided. Be sure to familiarize yourself with the operation and basic theory behind the routing model (Lohmann et. al. 1998a, b). If you have any questions about this documentation please email Bernt V Matheussen at vicadmin@hydro.washington.edu.

Basin Deliniation

This chapter describes how to find the basin boundaries. Several information sources are used, i.e. atlas maps, digital scanned basin boundaries and digital elevation models (DEMs).

Atlas

The first step in basin delineation should always be to locate the river of inerest in a good atlas. A recommended atlas for rivers in the United States is:
Atlas of River Basins of the United States
U.S. Department of Agriculture
Soil Conservation Service, WA D.C. 20250
Second Edition June 1970
 

It can be obtained from the University of Washington libraries.

The Atlas is used to get a rough feeling for the total area that your river covers. You can also define the region (maxlatitude, minlatitude, maxlongitude, and minlongitude) that covers the entire drainage basin for your river. This box can be used to reduce filesizes by clipping data from large area datasets (like those for the globe, or the entire U.S.) to the smaller region of interest.

ARC/INFO: To make a polygon box that covers your basin area use GENERATE in ArcInfo. You can type in the four corners that you need, project the polygon to the right projection and use this to clip the larger datasets. This polygon will be referred to as your BOX POLYGON troughout of this document.

Delineation Methods

There are two methods by which to delineate the watershed for the river of interest. The first is to extract the basin delineation from compiled databases of those delineations, the other is to use the watershed function in ARC/INFO to work out the flow of water across a DEM.

1. Extracting basin delineations from compiled databases

This sectons describes two databases of watershed delineations, and describes how to extract the boundary of your watershed. This method requires more hands on effort, but can produce good results fairly quickly.
    HUC-GCIP-CD
    There is a CD available that have ArcInfo coverages of all the large basins in the US. This CD is called.
    Hydrologic Units of the Conterminous US
    GEWEX Continental-scale International Project GCIP
    GCIP Reference Data Set (GREDS)
    by: Alan Rea & Joel R Cederstrand
    US.Geological Survey
    Open-File -Report 94 388
     

    Copies of this CD set are available to the local working group, otherwise contact the GEWEX program. The file xportarc/usbasin.gz on the CD-ROM contains the US basin boundaries. Transfer this file from the CD to the UNIX system and import it into ARC/INFO as follows:

    UNIX: gzip -d usbasins.gz to uncompress the file
    UNIX: mv usbasins usbasins.e00 to rename the file
    ARC/INFO: import auto usbasins.e00 basin1 to import the file

    Remember to CLEAN and DEFINE PROJECTION before using the cover in ARC/INFO. The projection for this file is Lambert-Azhimutal, details are given in the documentation on the CD.

    After defining the projection, the file should be reprojected to geographic coordinates (in degrees) so that it matches the coordinates used by the VIC model.

    Next clip the file using the BOX POLYGON created previously. This reduces the file size, and helps eliminate basin delineations which are not part of your watershed.

    Finally, open the clipped file in ARCEDIT. You can now manually select and delete basin boundaries which do not fall within the watershed of your river. To help in this selection process use the map from the atlas, the DEM (use it as the background for your editing session), or any other map image which can aid you in determining where the watershed boundary lies. It may also be possible to use the basin HUC codes to generate Another useful property is if the basin HUC codes that follow the basin polygons you can filter out the polygons which have a specific coding.

    Alternative basin boundaries
    The United States Geological Survey (USGS) has a webpage with an alternative source of basin boundaries for North America rivers. ArcInfo coverages can be downloaded directly from http://edcwww.cr.usgs.gov/;anddaac/gtopo30/hydro/na_basins.html as follows:

    Download the file na_bas.tar from this site.
    UNIX: gtar -xfv na_bas.tar to extract the file from the archive
    UNIX: gunzip na_bas.e00.gz to uncompress the file
    ARC/INFO: import auto na_bas.e00 na_bas1 (Lambert-Azhimutal projection) to import the file into ARC/INFO.

    Major basins in this database have a six digit Pfafstetter code assigned to them. The code for the Arkansas-Red is: Level 1 = 8, level 2 = 82 and 84. This information can be used to filter out the basins you want to include.

Once all of the polygons not in your basin have been deleted, you should select all the remaining polygons and MERGE them into a single basin polygon. Save the merged polygon as a new coverage file which will be refered to later as your TRUE BASIN BOUNDARIES file.

2. Delineating watersheds from Digital Elevation Masks (DEMs)

Traditionally basins have been delineated from DEMs simply by following the instructions in ARC/INFO. However, more recently it has been established that if maps of the actual river network are available they can be used to "burn" the river network into the DEM. This results in significantly improved basin delineations.

To burn the DEM using an ARC line file of the river network do the following:

ARC/INFO: grid01 = LINEGRID converts the line file into a grid file. Use the smae projection and resolution as the DEM.
ARC/INFO: grid02 = CON(ISNULL(river_file), 0, -100) creates a new grid file where the river channel has values of -100, everything else is set to 0.
ARC/INFO: grid03 = DEM + grid02 creates a new grid where the river network has been burned in, so DEM values in the channel are now decreased by 100 (units).

Make sure that the value you use to burn the rivers is big enough to force the proper drainage channels, otherwise watershed may force the basin to drain the wrong way, or include other basins (this is especially true in flat regions).

Once the rivers have been "burned" into the DEM, follow the ARC/INFO instructions to complete the basin delineation.

Fraction File

The fraction file in the routing model is used on the boundaries of the basin deliniation so that a gridcell can have less area than 1.0, or 100 %. This section explains how to produce the fraction file using the TRUE BASIN BOUNDARIES found in the previous section. There are several ways of doing this and only one way is presented here.

DEM GTOPO30

To produce the fraction file, a grid of fine resolution (1-2 km) that covers the area defined by your TRUE BASIN BOUNDARIES cover is needed. You can produce this grid using your TRUE BASIN BOUNDARIES and the POLYGRID command in ArcInfo. Another way of doing this is to clip your basin out from a digital elevation model(DEM) using the TRUE BASIN BOUNDARIES cover. One DEM that can be used in macro scale hydrological modeling is the USGS GTOPO30 DEM. This is a DEM with 1 km resolution and can be downloaded free from the web.

Click here to goto the GTOPO30 DEM from USGS
 

Download the DEM for the areas you need. The files from the website above is usualy named. w140n40_tar.gz, etc.

Use GUNZIP w140n40_tar.gz to uncompress the files

Use TAR -XVF w140n40_tar to get the actual datafiles.

Users of ARC/INFO or ArcView can display the DEM data directly after simply renaming the file extension from .DEM to .BIL. However, if a user needs access to the actual elevation values for analysis in ARC/INFO the DEM must be converted to an ARC/INFO grid with the command IMAGEGRID. IMAGEGRID does not support conversion of signed image data, therefore the negative 16-bit DEM values will not be interpreted correctly. After running IMAGEGRID, an easy fix can be accomplished using the following formula in Grid:

GRID: out_grid = con(in_grid >= 32768, in_grid - 65536, in_grid)

The converted grid will then have the negative values properly represented, and the statistics of the grid should match those listed in the .STX file. If desired, the -9999 ocean mask values in the grid could then be set to NODATA with the SETNULL function.

Look in the *.PRJ files to see which projection to use. Used this information when you PROJECTDEFINE in ArcInfo.

If you have more than one grids downloaded from the web, use the MERGE command in ArcInfo to merge them together

When you have projected and cleaned the DEM file, use the GRIDCLIP command in GRID to cut out the basin area from the DEM using your TRU BASIN BOUNDARIES.

GRID: gridclip whole_dem true_dem1 COVER true_basin_boundaries

Set all the values in the true_dem1 cover to ones (1)

GRID: true_ones1 = con(true_dem1 > 0,1,0)

Now you have to make a grid that covers a larger area than your real basin boundaries and set all values in this grid to zeros (0).

GRID: gridclip whole_dem true_box1 COVER box_polygon

GRID: true_zeros1 = con(true_dem1 > 0,0,1)

You now have two grids, one with only ones (1) in it, and one with only zeros (0) in it. Merge these two grids using the MERGE command.

GRID: true_ones_zeros = merge(true_ones1,true_zeros1)

Its important to have the grid with the ones first in the parentezis.

Now you must make a grid that have the final resolution you will set up your VIC model at. The grid must cover a larger area than the real basin boundaries and have a uniqe value in each cell. An easy way to do this is to use a spreadsheet. For example for the Columbia River (1/8th degree) a grid with 94 rows and 100 coloumns were made. Starting in the upper left corner the value 1 was put in, and going right and downwards increasing values for each gridcell. Row 1 col 1 have the value 2, row 1 col 2 have the value 2, etc. The area that this grid covers can be the same that your BOX POLYGON. Put an ArcInfo header with rows,cols,xllcorner, ect in the top of the file, and import it to ArcInfo using ASCIIGRID in ARC.

Project this grid to the correct projection and resample it to the same resolution as the DEM. Lets call the resampled file for cube_resampled1.

Now you must use the ZONALMEAN function in grid. This function calculates the mean of the values given in a zone. See HELP ZONALMEAN for more details.

GRID: zmean1 = zonalmean(cube_resampled1,true_one_zeros)

Resample this file to the modelling resolution. 1/8tth degree = 0.125 degrees

GRID: fraction1 = resample(zmean1,0.125)

The grid fraction1 is will now have a value from zero to one (0-1) in all the pixels. This value represents the fraction of the area of each gridcell that is in the basin. Convert this file into ascii in ArcInfo using:

ARC: gridascii fraction1 fraction1.asc

This file (fraction1.asc) is your fraction file and can be used directly as by the routing model. Take alook at it and check that all the gridcells in the middle of the basin have only ones (1) in it, and that the pixels on the boundaries have a number between zero (0) and one (1). The gridcells that are not part of the basin should all be zeros.

Elevation Maskfile

When processing meteorological data an elevation maskfile is needed. This is a file that have the same resolution that your final model will have. Each pixel have the mean gridcell elevation in it. The gridcells that are not in the basin have a defined VOID number in them. This is usualy a zero (0) or the number -9999.

The elevation maskfile is produced the exact same way as the fraction file but instead of using zeros and ones, use real elevations and zeros. Use the real boundaries cover to clip the DEM. Calculate the zonalmean using the cube_resampled1 and the true_dem1. Resample to your resolution and gridascii the file.

Flow Direction File

This section describes in detail how to get the final flownetwork file for the routing model. You will need the fraction1 grid from before as well as the DEM.

Cube Outline

Its now time to use the grid fraction1 that you prepaired in preveious section. Make a new grid of the fraction1 grid that has the same value in all gridcells that are in the basin. The fraction1 file should have the same resolution as you are building the VIC model for. BVM 1999 Arkansas-Red 1/8th degree.

GRID: cube_grid = con(fraction1 > 0,1,0)

Now use the GRIDPOLY command in ArcInfo to make a polygon that is an outline of your basin.

GRID: cube_outline = gridpoly(cube_grid,0.00)

This polygon is the cube_outline of your basin. On the boundaries it will look like cubes or squares, which is because you are using a grid that have the final model resolution.

DEM

Now use the cube_outline and cut out the basin from the DEM.

GRID: gridclip box_dem1 cube_dem1 COVER cube_outline

The cube_dem1 grid now has 15x15 1 km pixels in each 1/8th degree box. This is also true for the gridcells on the boundary of the basin. If you are using another model resolution you will have a different set of pixels in you gridcells. The reason for this is because some of the programs that will be used later. Needs the same number of 1 km pixels in each model resolution gridcell. Now we will use only the cube_dem1.

Now use the FILL command in GRID to fill all the sinks in the cube_dem1 grid. If the DEM have sinks in it the flowaccumulation routine will not work.

GRID: fill cube_dem1 cube_fill1

GRID: flowdir1 = flowdirection(cube_fill1)

GRID: flowacc1 = flowaccumulation(flowdir1)

ARC: gridascii flowacc1 flowacc1.asc

The file flowacc1.asc contains the accumulated flowvalues in each cell. The outlet of your river have the highest value

Developing Final Flownetwork

Bart Nijssen and Greg O'donnel have developed some simpel algorithms to create a flow network for the basin. A reference on this paper will be put in later. Maybe will put in the whole document on another page and link this one to it.

The six first lines in flowacc1.asc have to be removed. Save these lines in a new file for later use.

Try this: head -6 flowacc1.asc > head1.txt, then remove top 6 lines.

flowacc1 have to be converted into binary long integers with the program called convert.c

convert asci longint flowacc1.asc target1 ROWS COLS

ROWS and COLS should be number of rows and coloumns given in the header of the flowacc1.asc file. Look in head1.txt so see this.

The outfile target1 is then used in the program flowgen.c

Try this: flowgen target1 ROWS COLS flow_dir1 15 15

15 and 15 meens the number of 1 km pixels in one model resolution gridcell. For the Arkansas Red at 1/8th degree, there were 15x15 pixels in one gridcell.

flow_dir1 is your flow direction file.

Plotting and control of flownetwork

Generic Mapping Tools (GMT)

There is a GMT script available to plot the routing network.

MO_GMT.SCR is a GMT script that uses a fortran program called concoord.f. Compile this as

f77 -o concoord concoord.f

Then use: mo_gmt.scr flow_dir1 RESOLUTION (0.125)

A postscript file with the routing network is then produced.

The program is not perfect and some manual corrections must now be done. This can be done by use of maps, and or plot the rivers and do a rough estimate of the flowdirection. It shoudnt matter that much because its mostly the cells on the boundaries that have errors. These cells mostly point in the direction towards the center of the basin.

All the files needed can be found in :

/nfs/norway/home/tempbvm/web/rout/software/

Columbia River Network at 1/4th degree

ArcInfo and C - programs

There is also available other programs which can be used to plot your routing network. These programs are written by BVM 1999 and are ment to be a substitution to the programs by Nijssen and O'Donnell. The programs are written in C and works together with ArcInfo.

Gridnet.c

This program reads the direction file that is produced by the flowgen.c program. It then generates a textfile that can be imported into ArcInfo. Use GENERATE, INPUT filename, LINES The input file to this program have to have only integers in it. Void numbers have to be zeros (0), an example is given below

0 0 0 5
3 3 3 5
1 1 3 5
4 3 3 3

Usage: gridnet filename rows cols resolution > newfile

Sources of stream networks

There are digital coverages of the rivers for the US.

GCIP-Rivers
Download the rivers from
Hydrologic Units of the Conterminous US
GEWEX Continental-scale International Project GCIP
GCIP Reference Data Set (GREDS)
by: Alan Rea & Joel R Cederstrand
US.Geological Survey
Open-File -Report 94 388

Check also the HYDRO1 dataset where you downloaded the DEM

Preprocessing Scripts and Programs

There is a suite of programs (such as the flowgen.c program mentioned above) that have been written to facilitate the preparation of flow direction and fraction files for routing VIC output. These programs and scripts, as well as the README.route_prep file documenting a streamlined approach (specifically designed for using a 30 arc-second digital elevation file and aggregating output for a 1/8 degree VIC application) are available for downloading here.

VIC Administrator

Last modified: Sun Apr 11 18:04:46 PDT 1999