Basin Delineation
The first step in preparing input files for VIC is to create a list of grid cells contained within your basin. This page explains how to delineate a basin and get a list of grid cells at the desired spatial resolution, using ESRI Arc/Info on Windows. The basin delineation is also used as the basis for parameter files for the routing model.
- Find Basin Boundaries
- Delineate Basin using HYDRO1k and Arc/Info
- Other Delineation Methods
- Useful Arc/Info Commands and Tips
Find Basin Boundaries
This section 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 interest 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 may 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 la rge 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.
Delineate a Basin Using HYDRO1k and Arc/Info
Use Arc/Info GIS together with the known gauge location (e.g. from USGS) or predetermined drainage basin boundaries (e.g. from HYDRO 1k).
Obtain the HYDRO1k dataset
- Navigate to the USGS HYDRO1k page
- Select the continent containing your basin (e.g. North America)
- Select drainage basins
- Download basins dataset as a tar file in Arc/Info export format
Prepare Dataset for Importing into Arc/Info
On LINUX:
- untar: tar -xvf na_bas.tar
- unzip: gzip -d na_bas.e00.gz
On Windows:
- Use winzip or another utility to extract the contents of the file
Start Arc/Info
(On Windows) One way this can be done is:
- Open a command prompt (start menu -> run -> cmd.exe)
- Change directory (cd) to the directory where you have your workspace
- (at the prompt) arc
- (at the prompt) arctools
- a menu should appear
Import the Image
- Import the image into arcinfo:Arc:import auto na_bas.e00 na_bas
- You now have a coverage called na_bas that polygons for multi-leveled drainage basins.
Project the Coverage
Describe this coverage: Arc: describe na_bas
You can see that the projection of this basin is Lambert-Azimuthal. Eventually you will want to change the projection of this basin (to geographic in this example). This can be done at this time or after you have picked out your basin. It is easier to determine which basin you are interested in if you change the projection now, but it takes more time to do it this way.
To change the projection, it is easiest to use "arctools":
- Arc: arctools
- select "Command Tools"
- Edit -> Coordinates -> Project Coverage
- Input Cover: na_bas (or right click in space and select)
- Output Cover: na_bas_geo (name it whatever you want)
- PRJ source: click Define..., Hemisphere, Geographic, enter OK
- Close out all the windows
View this Coverage
You can view the coverage either in the module "ArcEdit" (type ae) or using arctools:
- Arc command: arctools
- Select "Edit Tools"
- File -> Coverage, Open
- In right column, select na_bas_geo
- In bottom left column, select ARC, hit OK
Clean the Coverage
Everytime you perform an operation on a coverage, it needs to be "cleaned" in order to re-establish the polygons.
- In Edit Tools Menu -> Arctools -> Commands...
- ARCEDIT: clean
- ARCEDIT: save
- Quit arctools and open again, this time plotting the polygons (the polygon option was previously not available)
- File -> Coverage, Open
- In right column, select na_bas_geo
- In bottom left column, select POLYGON, hit OK
Highlight Basin of Interest
(NOTE: The basins in HYDRO1k are described with up to six digits (or levels) called the Pfafstetter system. For more information see http://eros.usgs.gov/products/elevation/gtopo30/hydro/readme.html.)
Basin Numbering:
- All basins ending in the numbers 2,4,6,and 8 are the four largest tributaries within a basin at that level.
- All basins ending in the numbers 1,3,5,7, and 9 are the interbasins between the tributaries.
- All basins ending in 0 are closed basins (no outflow).
Highlighting:
- LEVEL 1 is the the largest scale, and LEVEL 6 is the finest.
- In the Edit Tools menu, go to Edit -> Attribute Selection...
- A window titled Logical Expresssion pops up:
- Our objective is to highlight the Columbia River Basin in the Current expression line, enter the following and then hit Apply expression:
- LEVEL1 = 2 (This highlights the Mackenzie River Basin)
- LEVEL1 = 4 (This highlights the Nelson River Basin)
- LEVEL1 = 6 (This highlights the St. Lawrence River Basin)
- LEVEL1 = 8 (This highlights the Mississippi River Basin)
- So, the columbia river is not one of the four largest basins in this dataset.
- Now, try the "interbasins":
- LEVEL1 = 1 (This highlights a large area of Northwestern North America - including the Columbia)
- LEVEL2 = 12 (This highlights the Columbia River Basin)
- LEVEL2 = 14 (This highlights the Fraser River Basin)
- LEVEL2 = 16 (This highlights the Kuskokwim River Basin)
- LEVEL2 = 18 (This highlights the Yukon River Basin)
- Continuing in this manner, smaller basins can be highlighted.
- Highlight the Columbia basin again:
- LEVEL2 = 12
Select Basin of Interest
Once you have the basin highlighted, merge all polygons in basin:
- In the Edit Polygon menu: MERGE
- In the Feature Selection menu: hit the Switch Selected Set button (two arrows)
- In the Edit Polygon menu: DEL
- Now clean the basin as before arctools\command then ARCEDIT: clean
- NOTE: After clean how many polygons you have:
- "Built 2 polygon(s) 1 of which have newly created labels" is GOOD (the basin and the background) , go the the next step
- "Built 3 polygon(s) 1 of which have newly created labels"
- This is NOT good, you have to manually select the basin, then do as previously; delete the rest and then clean until you have 2 polygons remaining
- Manually select the basin: click the arrow on the FEATURE SELECTION menu and select the basin. Then type 9 to get out of the selction mode.
- In the Feature Selection menu: hit the Switch Selected Set button (two arrows)
- In the Edit Polygon menu: DEL
- Clean, 2 polygons? if not,
- Click ALL on FEATURE SELECTION
- In the Edit Polygon menu: MERGE
- Clean (now it should work)
Save as something different (e.g. colum_dem); ARCEDIT: save "newname"
Quit arctools. You now have a delineation for your basin.
Create a Mask File
Load your half degree global mask file using asciigrid (newmask) - call it glo_mask.
- asciigrid newmask glo_mask
- Go into grid, then use the following commands:
- setmask glo_mask
- setwindow glo_mask
- NOTE: glo_mask is the global mask at 0.5 degree with 720 cols and 360 rows, xll=0, yll=-90
- If you are not interested in a global mask just do the following:
- setwindow xmin ymin xmax ymax (look at the arcinfo help)
- then setcell your desired value, 0.125, etc
- Use the polygrid command to convert the delineation into a mask
- basin_grd = polygrid ( basin_delineation , # , # , # , cell_size )
- Quit arctools
Use gridascii to export it: gridascii basin_grd basin_mask.asc
NOTE: As an example, what if you would like to have a cell number starting at 0 instead of 1.25?
Type setmask glo_mask then change the setwindow command: setwindow -181.25 -91.25 178.75 91 .25
Clip Out a DEM for the Basin
As an example, if you have a DEM for North America (from gtopo30 or HYDRO1k) called na_dem, then use the following commands:
Note that DEM in HYDRO1k are rasters. You need grid: imagegrid
- Enter the grid module: grid
- Clip the DEM using the basin delineation:
- gridclip na_dem colum_dem cover colum_del
- View the clipped DEM:
- Bring up the graphics window: display 9999 or &stat 9999
- Set the map extent: mape colum_dem
- Set the color scheme: shadeset rainbow
- Plot the DEM: grids colum_dem # linear
- Plot the coverage: arcs colum_del
- Quit grid: quit
- export your basin DEM to an ascii file: gridascii colum_dem colum_dem.asc
Other Delineation Methods
HUC-GCIP-CD
- There is a CD available that has ArcInfo coverages of all the large basinsin 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 f ollows:
- 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
- Hydrologic Units of the Conterminous US
- Remember to CLEAN and DEFINE PROJECTION before using the cover in ARC/INFO. The projection for this file is Lambert-Azhimutal, deta ils 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 t he 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 ima ge 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 pr operty is if the basin HUC codes that follow the basin polygons you can filter out the polygons which have a specific coding.
Delineating from Digital Elevation Models (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 improv ed 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 isset 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 todrain 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.
Useful Arc/Info Commands and Tips
- Smooth out a grid: outfile = RESAMPLE() in GRID
- Transform a delineation into polygon : outfile = POLYGRID ( delineation , # , # , # , cellsize ) in GRID
- For the preparation of VIC routing files, one needs the clipped DEM. The gridclip command will unfortunately change the size of the DEM window to the often weird one of the delineation (cover file). The trick is to:
- setwindow the DEM window (if the DEM window was correct for you, i.e. the number of cols/rows % number of cells to be aggregated for the new resolution (.25, .5 etc) is correct)
- Use polygrid (cover, # , #, #, dem res)
- Use: clipped_dem = con (cover_grid == 2 , dem)
- Transform a grid into polygons : outfile = GRIDPOLY ( grid , # )
- Get a delineation in ready-to-plot GMT format:
- In Arc/Info, use the UNGENERATE POLY command. If you encounter bad results, use LINE instead
- In LINUX: awk '{ if (NF!=2) print ">" ; else print $0;}' basin.txt >! basin_ready.txt
- In GMT: psxy basin_ready.txt -M -R$COORD $PROJ -K -O -W3 -V >> $OUTF
- The -M allows it to close the polygons.
- You have a GRID and want to know the value : in COMMANDS: cellvalue gridname * then select the cell. Type 9 to stop.
- Manually select a basin (in a basin coverage): in EDIT, FILE/open coverage in arcs or polygon, then click on the arrow. 9 to get out of the mode.
- When doing the delination, make sure there are only 2 polygons remaining after merging
- In GRID:
- mape defines the map extent
- grids plots the grid
- mape * will give a different map extent, and use grids to replot with the new map extent
- arcs or polys will plot the basin delination of the previous grid
- gridclip backgroundmap outfile cover basin_delination
- Quit GRID and then gridascii to export it
- use arcview to create nice maps
- Copy an Arc/Info file from someone else:
- You not only need to copy the directory (if grid), but also the info directory and create a log file. Therefore you need to be in >arc to do it, then >arc: copy "entire path and name" "newname"
- In order to get or sent an Arc/Info or Arc/GIS files to/from somebody, in my experience it is better to use the Interchange File Format that you use by entering EXPORT OPTION(AUTO, COVER, GRID, etc.), FILENAME OUTFILENAME or by IMPORT OPTION(AUTO, COVER, GRID, etc) FILENAME
- lg states the grid in the directory
- lc states the coverage in the directory
- kill name removes the file (including the entire file directory if grid)
- asciigrid imports ascii format (used in UNIX) to grid (ArcInfo/ArcView)
- gridascii exports grid files into ascii files
- Create a buffer for a grid: outgrid = EXPAND (grid, # of cell desired for the buffer, LIST , value1, value2, value3, etc.)
- Look at the help for this command for mode detaild. In particular, one can use FILE instead of LIST
- CON in the grid menu: it transforms a value by another one using a conditionnal statement: outfile = con (isnull(infile),0,infile) for example to change -9999 into 0 in a maskfile
- Note that Arcview reads grid files and coverages. It is easier to transform the projection in Arc/Info than in Arcview, (I found anyway, especially for large files)
- http://gis.washington.edu/phurvitz/arc_lms/ has help for GIS Spatial Analyst and aml language
VIC Administrator Last modified: Tue Aug 11 14:57:28 PDT 2009
