Strahler stream order in GRASS

In hydrology, a stream network is composed of segments or “reaches” which are arranged in a hierachy. There are several systems of ordering the stream reaches, the most popular of which is the Strahler or Horton number. GRASS GIS offers, alongside the watershed delineation tool r.watershed (discussed here), a set of addons for stream network analysis. We’ll examine how to use these addons, and how to use strahler ordering to improve the visual effect of a stream network map.

We begin with an elevation raster call ‘dtm’. This could be data from the ASTER GDEM program or SRTM tiles. Zoom in to the general area of interest and (as always) set the GRASS region to match that area.

g.region -p n=30:54N s=30N e=35:30E w=34:36E

First we run two r.stream.* commands to create the stream network, and a database table of stream ordering.

The r.stream.extract command uses the dtm as input, together with a threshold value which we determine. The threshold represents the number of raster cells which will become a minimim drainage catchment in the resulting stream network. You need to know the size (resolution) of the cells in your dtm raster in order to do the arithmetic for this threshold value. In the example above, the dtm (based on ASTER GDEM data) has a resolution 30 meters, so 25,000 cells is about 20 sq. km. Then the r.stream.order command uses the newly created streams_25 map of streams to create a table of stream ordering. The command places the table into the default database used by GRASS. In this case sqlite.

So the vector stream map has only three columns – the cat value, and a name and code indicating if the stream segment is a headwater segment, or further down in the hierarchy. To examine the stream_order table, we need to have a look right in the sqlite database, so:

So this second table has all the good stuff – next stream cat number, previous streams, the ordering in four different systems, as well as length, cumulative length for each stream reach and so on. We want to attach this table to the streams_25 map. This is accomplished with v.db.connect. The ‘-o’ flag indicates that we are overwriting the original table connection. And the ‘key=cat’ parameter says that the cat values in the geometry match the ‘cat’ column in the sqlite table. Now we quit sqlite, and back in GRASS we issue:

All the stream ordering and other data is now part of the streams_25 vector map. We next run the r.stream.basins command to create a map of the basins. The ‘-l’ flag is set so that only major basins (with a outlet at the edge of the current region) are delineated.

Now we’re ready to make use of the data and maps we’ve created. We create a shaded relief, and overlay the basins with transparency. In addition, In the GRASS GUI, you can display a line vector at varying widths, taking the width from an attribute column. We’ll use the strahler order values to make stream vector lines at the top slopes of the basin appear thin, and those nearer to the basin outlet will be thicker. Here’s the result.

Faran basin – wadis displayed with strahler order

Often you’ll need to isolate one particular sub-basin and display only that basin and it’s stream network. Here’s a way to do that. We first display the stream_25 raster map, and zoom in very close to the outlet of the sub-basin we’re trying to isolate. We need to get the exact X-Y coordinates at the outlet of that basin. In the GRASS GUI Map Display window we can use the “Query raster/vector map” button (the gui version of the r.what command) to get the coordinates by clicking on the drainage point. The coordinates appear in the “Command console” of the Layer Manager window. Copy/paste these coordinates to the terminal and run r.stream.basins a second time as follows:

r.stream.basins dir=fdir_25 coors=35.226820,30.632014 basin=basin_1

Display this map in the Map Display to be sure it’s the basin you need. Now it’s a good idea to reset the region parameters to “close-in” on this basin, in order to avoid large raster maps with NULL cells covering large areas around your isolated basin.

We’ve created new rasters for the isolated basin, the dtm, and the shaded relief, all clipped to the area of the isolated basin by using a raster MASK. After completing those steps the mask is no longer necessary, we can remove it, and the temporary map ‘basin_1′ as well. Then use the regular GRASS commands to create an area vector of the isolated basin, and to clip the original streams vector to this basin.