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