Watershed Basin Delineation

Multiple approaches to the automatic delineation of watershed boundaries, along with the processing associated with the CalWater 2.2.1 dataset. Note that the output from each algorithm is slightly different, and that CalWater does not always represent physiographicaly correct boundaries. The CalWater 2.2.1 original dataset and converted to shapefile are available for download. See attached DEM clip for examples presented below. Matt Perry has written up a nice tutorial on identifying impervious surfaces within the context of watershed delineation.

Watershed basin delineation by the following approaches:

ArcInfo basin() GRID function

multiscale approach suitible for small to medium input grids. raster output.

GRASS r.watershed command

minimum basin size must be defined suitible for small to medium input grids. raster output.

GRASS r.terraflow command

multiscale approach suitible for massive input grids, fast. raster output.

Notes

none of these algorithms can deal with hierarchial watershed decomposition
removal of small vector features does not ensure perfect basin delineation
raw basin output from these algorithms are best further processed some ideas here

watershed 1: 3D comparison of the four methods. calwater (red), basin (black), r.terraflow (white).
watershed 2: 3D comparison of the four methods. calwater (red), basin (black), r.terraflow (white).

Pre-processing in GRASS note that these steps can be done in any GIS

# projection settings: +proj=aea +x_0=0.0 +y_0=0 +lon_0=-96 +lat_0=40.0 +lat_1=20 +lat_2=60.0 +datum=NAD83
# set the region and cut out a sample section of the DEM
g.region res=90 -ap
r.mapcalc "a = ca_90m_dem_meters"
#generate a shaded relief map for viz:
r.shaded.relief map=a shade=a.shade

ArcInfo basin()

# start ArcInfo and import our DEM called 'a_grid'
# note that 'a_grid' was exported from GRASS via: r.out.arc in=a out=a_grid
asciigrid a_grid a float
grid
basins = basin(flowdirection(a))

GRASS r.watershed

# run r.watershed, saving basins in 'w.basins'
# note that a threshold of 200 sets the minimum basin size [update this]
r.watershed --o elev=a threshold=200 basin=w.basins

GRASS r.terraflow

# run r.terraflow, saving basins in 'a.swater'
r.terraflow elev=a filled=a.filled direction=a.dir swatershed=a.swater accumulation=a.acc tci=a.tci
#set cell values of 0 to NULL, these are regions where calculation was not possible
r.null map=a.swater setnull=0

Post-processing of generated basins in GRASS note that this isn't a true merging of sub-basins, just a rough idea of what it might look like!

# vectorize watershed basins from two of the methods:
# r.terraflow
r.to.vect -s in=a.swater out=a_watersheds feature=area
# ArcInfo basin()
r.to.vect -s in=arc.basins out=arc_basins feature=area

# remove basins smaller than approx. 200ac.
# r.terraflow
v.clean --o in=a_watersheds out=b_watersheds tool=rmarea thresh=809371
# ArcInfo basin()
v.clean --o in=arc_basins out=arc_basins_clean tool=rmarea thresh=809371

CalWater note that this comes as an e00 file, and must be converted.

# get from http://casil.ucdavis.edu/casil/gis.ca.gov/calwater/calw221_20041118.zip
# note that this is an e00 file and needs to be converted

# converted with ArcInfo
import AUTO calw221.e00 calwater
# converted with AVCEimport http://avce00.maptools.org/avce00/index.html
~/src/avce00-1.3.0/avcimport calw221.e00 calwater_avc

# ArcMap simple conversion to shapefile and re-projection possible via standard tools
# converted to shapefile polygon

# conversion to shapefile with attributes intact not as simple...

# projection to our spatial ref system via ogr2ogr
# note that this shapfile was created by ArcMap
ogr2ogr -t_srs '+proj=aea +x_0=0.0 +y_0=0 +lon_0=-96 +lat_0=40.0 +lat_1=20 +lat_2=60.0 +datum=NAD83' calwater_aea.shp calwater\ polygon.shp

Visual comparison in GRASS

# cut out small piece of the CalWater vector for analysis
v.in.region out=temp_region
v.overlay ainput=calwater binput=temp_region out=calwater_clip operator=and
g.remove vect=temp_region

## compare the methods:
d.his i=a.shade h=w.basins # r.watershed
d.vect calwater_clip type=boundary col=red width=2 # Calwater
d.vect arc_basins_clean type=boundary col=black width=1 # ArcInfo basin()
d.vect b_watersheds type=boundary col=white width=1 # r.terraflow

Aggregation of CalWater Data

#aggregation of CalWater data: see metadata:
Relevant columms:
      43  HRC           2     2     C     -    Hydrologic Region Code (DWR)
      45  HBPA          2     2     C     -    Hydrologic Basin Planning Area (RWQCB)
      47  RBU           5     5     I     -    Concatenates HR,RB,HU into a single integer
      52  RBUA          6     6     I     -    Concatenates HR,RB,HU,HA
      58  RBUAS         7     7     I     -    Concatenates HR,RB,HU,HA,HSA
      65  RBUASP        9     9     I     -    Concatenates HR,RB,HU,HA,HSA,SPWS
      74  RBUASPW      11    11     I     -    Concatenates HR,RB,HU,HA,HSA,SPWS,PWS
      85  HR            2     2     I     -    Hydrologic Region (1->10)(DWR)
      87  RB            1     1     I     -    Regional Water Qual. Cont. Board (1->9)(RWQCB)
      88  HU            2     2     I     -    Hydrologic Unit (00->~80)(SWRCB)
      90  HA            1     1     I     -    Hydrologic Area (0->9)(SWRCB)
      91  HSA           1     1     I     -    Hydrologic Sub-Area (0->9)(SWRCB)
      92  SPWS          2     2     I     -    Super Planning Watershed (00->~30)(CDF)
      94  PWS           2     2     I     -    Planning Watershed (00->~13)(CDF)
      
# polygons can be dissolved via ArcMap / PostGIS / v.to.rast use=attr -> r.to.vect

regarding watershed decomposition and limitation C/O Helena Mitasova

Dylan,

r.terraflow has not been designed to output landscape decomposition into
watersheds at a given scale/size but its result is used as an input for
hierarchical watershed decomposition - see here

http://terrain.cs.duke.edu/projects.php
(and the papers by Arge et al. and Verdin and Verdin)

However, the method used there can still leave you with small internal
subwatersheds (r.watershed will have that too) that cannot be lumped  
with
any bigger watershed.

Helena