Watershed Analysis with GRASS

Arcview users, needing to delineate watersheds and stream networks, would do well to choose the extension called  “Arc Hydro” (requires, at least the Spatial Analyst extension from ESRI). This extension introduces the concepts of “Batch Points” and “Adjoint Catchments”.  Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points.

I’d like to demonstrate how to get similar results with GRASS GIS. As an example, we’ll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We’ll assume our starting data includes an elevation raster called “dem” and a line vector called “roads”.  We first create the regular hydrology layers.

#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>

So far, pretty straightforward. There’s abundant information on r.watershed. I’ll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10×10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

When the process finishes we’ll have three new raster maps: the flow direction map, the streams and the catchments. Let’s see what we’ve got so far:

#Convert the steams and catchments to vectors
r.to.vect in=catch out=catchments feature=area
# the stream raster usually requires thinning
r.thin in=str out=str_thin
r.to.vect in=str_thin out=streams feature=line
r.colors dem col=elevation
# Make a hillshade raster for displaying "3D"
r.shaded.relief map=dem shade=dem_shade zmult=1.5
# Now display layers
d.mon start=x0
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2

Now we need to find all the points where streams cross roads. The v.overlay module does not deal with point vectors. Instead we use a trick in v.clean. When cleaning a line vector, all points where lines intersect and no node exists are considered “errors” and can be saved to a new point vector. So by merging the roads and streams vectors, we have lines (streams) crossing lines (roads) with no node at that location. Then we run v.clean to get all those intersection points in a new layer.

# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points
#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12
# Save crossing points to a text file
v.out.ascii in=cross_points out=cross_points.txt format=point fs=space

Once we have the crossing points in a file, we simply run r.water.outlet in a loop to create the watersheds for each cross point.  The module r.water.outlet sets the raster value to 1 in the watershed and 0 elsewhere. We’ll use r.null to force the raster to ‘null’ outside of the watershed. THen we convert each raster to a vector (polygon).

i=0   #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 ))  \
r.water.outlet drain=fdir east=$X north=$Y basin=wshed_$i  \
r.null -q wshed_$i setnull=0 \
r.to.vect -s in=wshed_$i out=wshed_$i feature=area \
g.remove rast=wshed_$i
done < cross_points.txt
echo "Created $i catchments"

We patch together all the area vectors and clean the merged watershed vector

v.patch in=`g.mlist vect pattern=wshed_* separator=,` out=wshed_patch
v.clean wshed_patch out=wshed_clean tool=snap,rmdupl thresh=20

Choose an appropriate snap threshold value based on your region resolution.  Additional manual cleaning may be required. After cleaning, reset the cat values (only for centroids):

v.category in=wshed_clean out=wshed_clean2 opt=del
v.category in=wshed_clean2 out=wshed_final type=centroid opt=add
g.remove vect=wshed_clean,wshed_clean2

Most likely we’ll want to calculate the area for each watershed.

# Add a table with a column for area in sq.km.
v.db.addtable wshed_final col="area_sqkm DOUBLE"
v.to.db map=wshed_final option=area col=area_sqkm unit=k

And finally, we can view the catchments, and their area values:

d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange

13 thoughts on “Watershed Analysis with GRASS

  1. I don’t get how you run the loop in GRASS, is that Python code and if so, how do you get it to work in GRASS? To me it is not intuitive to run r.water.outlet in a loop. Can you please explain to me a bit more how are you doing this and where I must input your code in the new GRASS with Python. Thank you. S

  2. Stephanie :

    I don’t get how you run the loop in GRASS, is that Python code and if so, how do you get it to work in GRASS? To me it is not intuitive to run r.water.outlet in a loop. Can you please explain to me a bit more how are you doing this and where I must input your code in the new GRASS with Python. Thank you. S

    Hi Stephanie:
    No it’s not Python, although you certainly could do the same loop in python.
    The example above is a simple bash script. If you enter the commands in a terminal window on linux (or a GRASS command prompt in windows) it should work. If you have the drainage outlet points in a text file ‘cross_pts.txt’ then the loop structure:
    while read X Y; do
    .... ;
    done < cross_pts.txt

    reads each line separately till the end of the file and allows you to use the X,Y values for each individual run of r.water.outlet.
    HTH,
    Micha

  3. Hi Micha,

    Thanks so much for the quick reply. I have tried to run this in the GRASS command prompt in windows but it continues to return an error. Perhaps I am doing it wrong. I have just copy pasted your code into the Windows prompt window that is part of the GUI. I also tried it in the command line – still did not work for me. Also, I noticed that in my cross_points.txt file I don’t have any headers – is that problem? I was expecting that Grass would output a .txt with headers, but that is not the case. Perhaps I need more patience. Now that I know the code is right for Grass I will keep trying, as I must simply just be doing some step wrong. Thanks again!

  4. Hi Micha,
    I am not sure what version of GRASS you developed the code for, but I am working on Windows 6.4.1 and the code is not working. I have tried in the GUI as well as in the command line and its not working at all. I get a list of errors indicating that I have not defined things correctly and that i is not recognized. The only thing I can think is that currently R does not talk with Grass in this version. Did you require GRASS to talk to R when you wrote this? Any other suggestions would be really appreciated, thanks for your time.

  5. Hi again,

    Sorry, I think I got it! The bash file wasn’t in the same directory as the data – hence it was not liking anything I was asking it do. All set for now and it appears to be running! Fingers crossed.

    Thank you

  6. Hi,
    I am having some trouble with different parts of the script. I often get an error saying that a file cannot be found or created in the mapset or other. I am trying to run a .sh script in the GRASS command window through a virtual ubunto operation system. The commands work when typed into the command window.
    Thanks,
    Tom

  7. Hi Micha,

    I’ll try be more specific. The error generally reads “No such file or directory” at the line which references the text file of X Y co-ordinates. I have tried with many different file paths and extensions in the script. Also with a dummy set of co-ordinates to test. All with no success.

    Thanks,
    Tom

  8. Tom :

    Hi Micha,

    I’ll try be more specific. The error generally reads “No such file or directory” at the line which references the text file of X Y co-ordinates. I have tried with many different file paths and extensions in the script. Also with a dummy set of co-ordinates to test. All with no success.

    Thanks,
    Tom

    And you say that when you run the commands one by one, everything goes as expected? In that case it’s probably some syntax error, or so. Would you care to share the script somewhere so we can have a look, and (maybe) make some suggestions?

    Micha

  9. Hello,

    Thanks for the exercises which prove to be a very good tool to learn QGIS. I faced a problem though: with QGIS 2.0.1, whenever I try to delineate a watershed with the module r.water.outlet, I receive a black screen covering my Grass vector layer (surrounded with a red line). I did several attempts and even used the “coordinate capture” tool to define precisely the coordinates of a stream outlet, but with no result. The raster layer that I am using as a source is “flowdir”. Your help will be very much appreciated. Many thanks.

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>