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.

r.stream.extract elev=dtm thresh=25000 stream_rast=streams_25 stream_vect=streams_25 direction=fdir_25 
r.stream.order stream=streams_25 dir=fdir_25 table=stream_order

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.

Before continuing, we’ll have a look at what we’ve got.

v.db.select streams_25 | head 
cat|stream_type|type_code 
1|start|0 
2|start|0 
3|start|0 
4|intermediate|1 
5|start|0

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:

eval `g.gisenv` 
sqlite3 $GISDBASE/$LOCATION_NAME/$MAPSET/sqlite.db 
sqlite> .header on

A quick SELECT statement on the streams_25 table (attached to the stream_25 vector map) will show the same data as above.

sqlite> SELECT * FROM streams_25 LIMIT 5; 
cat|stream_type|type_code 
1|start|0 
2|start|0 
3|start|0 
4|intermediate|1 
5|start|0

Next we do a SELECT on the stream_order table:

sqlite> SELECT * FROM stream_order LIMIT 10; 
cat         stream      next_str    prev_str01  prev_str02  strahler    horton      shreve      hack        topo_dim    length      cum_length  out_dist     straight    fractal 
----------  ----------  ----------  ----------  ----------  ----------  ----------  ----------  ----------  ----------  ----------  ----------  -----------  ----------  ---------- 
1           1           20          0           0           1           1           1           2           38          929.78593   929.78593   7775.778998  676.39042   1.374629 
2           2           4           0           0           1           3           1           1           29          300.43355   300.43355   8053.78109   233.019313  1.289308 
3           3           4           0           0           1           1           1           2           29          3.828427    3.828427    7757.175967  3.605551    1.061815 
4           4           10          3           2           2           3           2           1           28          667.10974   967.543289  7753.34754   549.315938  1.214437 
sqlite> .quit

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:

v.db.connect -o streams_25 driver=sqlite table=stream_order key=cat

If we now run a SELECT on the streams_25 map, we find:

v.db.select streams_25 | head 
cat|stream|next_str|prev_str01|prev_str02|strahler|horton|shreve|hack|topo_dim|length|cum_length|out_dist|straight|fractal
1|1|20|0|0|1|1|1|2|38|929.78593|929.78593|7775.778998|676.39042|1.374629
2|2|4|0|0|1|3|1|1|29|300.43355|300.43355|8053.78109|233.019313|1.28930
3|3|4|0|0|1|1|1|2|29|3.828427|3.828427|7757.175967|3.605551|1.061815
4|4|10|3|2|2|3|2|1|28|667.10974|967.543289|7753.34754|549.315938|1.214437
5|5|19|0|0|1|4|1|1|39|1050.680374|1050.680374|8229.62171|817.765247|1.284819
6|6|12|0|0|1|1|1|2|26|215.865007|215.865007|7098.509878|174.103418|1.239867

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.

r.stream.basins -l dir="fdir_25" stream="streams_25" basin="major_basins"

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
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.

g.region n=30:45N s=30:18N e=35:15E w=34:36E -p
r.mask in=basin_1 
r.mapcalc nekarot_basin=basin_1 
r.mapcalc nekarot_dtm=dtm 
r.shaded.relief nekarot_dtm shade=nekarot_shade

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.

r.mask -r 
g.remove rast=basin_1 
r.to.vect in=nekarot_basin out=nekarot_basin feature=area 
v.select --o ainput=streams_25 atype=line binput=nekarot_basin btype=area output=nekarot_streams operator=overlap

And the result looks like this:

Nekarot basin
Nekarot basin – wadis displayed with strahler order

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>