Better slope-area plots

Posted on

Plotting upslope area vs. slope in logarithmic axis has been the de-facto standard approach to river-profile analysis. In theory, if a river profile is in a perfect steady state, then this plot should yield a straight line. The TopoToolbox function slopearea calculates these plots. In addition, the function fits a power law S = k*A^(-theta). However, the function has drawbacks, and I want to address one of these in this post.

Slope-area plots require stream gradient, the first derivative of the river profile. If elevation values of the profile have errors, then this uncertainty becomes even more significant in the first derivative. To alleviate this problem, slope is usually averaged in bins. The function slopearea calculates the width of these bins based on upslope area so that bins are equally spaced in log space. However, this approach entails that bin width increases for larger upslope areas. Thus, for larger upslope areas, gradients are averaged over larger stream distances which may hide some of the structures that we wish to detect.

Another approach is to calculate slope averages over river segments of predefined length. While this is not supported by the function slopearea (yet), I’ll show here how to calculate it. The new function labelreach comes in handy here.

Let’s first load a DEM, derive a stream network.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM,'preprocess','carve');
S = STREAMobj(FD,'minarea',1000);
S = removeshortstreams(S,1000);
DEM = imposemin(S,DEM);
A = flowacc(FD);

Then we calculate node-attribute lists for upslope area and gradient.

a = getnal(S,A)*(DEM.cellsize^2);
g = gradient(S,DEM);

The function labelreach calculates a label for each reach in a STREAMobj. The option ‘seglength’ allows us to subdivide the labels into reaches of ~2000 m length. Note, however, that the function separates reaches at confluences so that segments may not always be exactly 2000 m.

label = labelreach(S,'seglength',2000);

The function accumarray is one of my favorite functions. Here we use it to aggregate the values in the node-attribute lists g and a based on label. As an aggregation function, we use the average, while we express errors as the standard error.

gg = accumarray(label,g,[],@mean);
ggs = accumarray(label,g,[],@(x) std(x)/sqrt(numel(x)));

ag = accumarray(label,a,[],@mean);
ags = accumarray(label,a,[],@(x) std(x)/sqrt(numel(x)));

Finally, let’s plot the results.

plot([ag ag]',[gg+ggs max(gg-ggs,0.0001)]','color',[.7 .7 .7]);
hold on
plot([ag-ags ag+ags]',[gg gg]','color',[.7 .7 .7]);
plot(ag,gg,'o','MarkerFaceColor',[.7 .7 .7])
hold off
xlabel('Area [m^3]')
ylabel('Gradient [m m^{-1}]')
ylim([1e-3 0.6])
Distance-binned slope-area plot of rivers in the Big Tujunga area

Now that looks much better than area-binned slope-area plots. Give it a try…

Stream networks and graph theory

Posted on Updated on


Rivers are organized as networks. They start out as small rills and rivulets and grow downstream to become large streams unless water extraction, transmission losses or evaporation dominate. In terms of graph theory, you can think of stream networks as directed acyclic graphs (DAG). The stream network of one drainage basin is a tree routed at the outlet. To be more precise, this kind of tree is an anti-arborescence or in-tree, i.e. all directions in the directed graph lead to the outlet. Usually, stream networks are binary trees, i.e. there is a maximum of two branches or parents entering each node. However, nodes may sometimes have more than two branches at confluences where three and more rivers meet. All drainage basins form a forest of trees and each tree is a so-called strongly-connected component.

Why bother with the terms of graph theory? Because it will make a couple of TopoToolbox functions much more accessible. The function klargestconncomps, for example, has an awkward name that I adopted from graph theory. What does it do? It extracts the k largest trees from a stream network. The function STREAMobj2cell takes all trees and puts each into a single element of a cell array. FLOWobj2cell does the same for an entire flow network. The computational basis of STREAMobj and FLOWobj relies on topological sorting of directed graphs. Thanks to graph theory, this trick makes flow-related computations in TopoToolbox much faster than in most common GI-systems.

Graph theory has a lot to offer for geomorphometry and I think that this potential has not yet fully been explored.

Colored longitudinal river profiles

Posted on

Did you ever face the problem that you wanted to plot additional variables together with the two spatial dimensions of a river profile. One such variable could be ksn values, or gradient, or lithological identifiers. I recently had this problem, then coded a solution, and now uploaded an updated version of STREAMobj/plotdz.

Such additional variables must be available as so-called node-attributed lists (nal). These lists are column vectors that have as many elements as there are nodes in the stream network. A STREAMobj S has a couple of such nal as properties, e.g. S.x, S.y or S.IXgrid. S.x and S.y are coordinate vectors that locate each node in a projected coordinate system. S.IXgrid is a vector of linear indices into the DEM from which S was derived. nal are returned by a number of STREAMobj methods such as STREAMobj/gradient, STREAMobj/streamorder, STREAMobj/curvature, and STREAMobj/chitransform, among others. Alternatively, you can derive a nal using getnal. For example, a nal with elevation values is simply obtained by

z = getnal(S,DEM);

Now how can you plot colored river profiles. Line objects as produced by the function plot do not support gradually varying colors. As far as I know, there are three work-arounds. The first approach is to plot each line segment with variable colors using the low-level function line. This, however, would probably take ages for large stream networks with thousands of individual line segments. The second approach is to plot the line as a surface with colored edges. This old-time-hack works quite well and leaves you with all the options that surface objects offer. The third approach is one that I came across when reading through Yair Altman’s undocumented MATLAB blog. He reports on some undocumented line properties that were introduced with MATLAB’s new graphics engine in R2014b. I implemented the surface approach and the undocumented approach and let you choose.

Ok, now let’s see the updated plotdz-function in action:

load exampleDEM
DEM = GRIDobj(X,Y,dem);
FD = FLOWobj(DEM,'preprocess','carve');
A = flowacc(FD);
S = STREAMobj(FD,A>100);
S = klargestconncomps(S);
z = imposemin(S,DEM);
s = streamorder(S);
plotdz(S,DEM,'color',s,'LineWidth',1,'cbarlabel','Stream order')

g = gradient(S,z);
c = chitransform(S,A,'mn',0.45);
plotdz(S,DEM,'color',g,'LineWidth',1,'distance',c,'cbarlabel','Gradient [m/m]')
xlabel('\chi [m]')
Colored longitudinal river profiles using STREAMobj/plotdz

Now that looks quite ok. Also note the option to modify the values of the upstream distance that allows you to easily plot Χ-plots without using the function chiplot.

Hope that functions turns out to be useful for you.

Orientation of stream segments

Posted on Updated on

Today’s guest post comes from Philippe Steer, assistant professor at the University of Rennes. Philippe contributed two new functions and I asked him to describe them here to show potential applications. Thanks, Philippe, for this contribution!

I recently needed to compute river network orientation. Indeed river network orientation can reveal distributed tectonics deformation at the vicinity of strike-slip faults (e.g. Castelltort et al., 2012) or can enable to identify developing vs paleo river networks.  One way to obtain orientation is first to extract the river segments that are linking confluences to outlets or confluences to channel heads and then to compute the orientation of the resulting segments. The streampoi function allows to extract these feature points of the network such as channel heads, confluences, b-confluences (the points just upstream any confluence) and outlets. To connect these feature points into segment, one can use the propriety that each segment belongs to a unique watershed that is given by the function drainagebasins using as outlets the true river outlets and the b-confluences. By definition, each segment also corresponds to a unique value of strahler order.

Once the river segments are identified, segment orientation is easily obtained. What is interesting is that defining these segments also paves the way to compute several other geometric metrics averaged over the river segment such as segment direct or flowpath length, sinuosity, direct or flowpath slope, mean drainage area and even steepness index (if assuming a value of concavity).

To cleanly do the job, I have written two MATLAB functions:

Here is a piece of code showing you the use of these functions (assuming you have already loaded the DEM of interest):

%% Flow direction and drainage area
DEMf = fillsinks(DEM); FD= FLOWobj(DEMf); A = flowacc(FD);  

%% Extract rivers
% Set the critical drainage area
A_river=1e7/(res^2); % m^2
W = A>A_river;
% Extract the wanted network
S = STREAMobj(FD,W);

%% River segments
% Extract segment and compute their orientations
% Plot results

To illustrate the results, I thought the South Island of New-Zealand was a good example. Indeed, as first suggested by Castelltort et al. (2012) the river catchments at the south of the Alpine fault are tilted, possibly suggesting that they have recorded a distributed shear deformation. This tilting is also quite obvious when looking at the river segments in the following figure. Another interesting feature is that the river segments north of the Alpine fault, facing the tilted rivers, display a higher steepness than their south neighbours. This could be due to a long list of forcings, including active tectonics, orographic precipitation gradient or the transition from glacial to fluvial erosion.

Castelltort, S et al. (2012). River drainage patterns in the New Zealand Alps primarily controlled by plate tectonic strain. Nature Geoscience, 5(10), 744-748. [DOI: 10.1038/ngeo1582]

Side branching

Posted on Updated on

The Strahler stream order is a common way to classify rivers. TopoToolbox has the tools (e.g. STREAMobj/streamorder, FLOWobj/streamorder) to easily derive this stream classification from a river network. However, occasionally this information may not suffice to describe the branching behavior of river networks. Instead, we may wish to quantify which river-orders intersect how often with other river orders. E.g., if a first-order stream intersects with another first-order stream, we may wish to denote that stream as ’11’; if a first-order stream intersects with a third-order stream, we denote that one as ’13’. This is the method of Tokunaga (1978) which extends on the Strahler (1957) order system (Pelletier and Turcotte, 2000).

When Arvind – a TopoToolbox user from India – asked me, if it was possible to derive side-branching in TopoToolbox, my first thought was that this is easy: Just take the output of STREAMobj2mapstruct, that – if called with only one output argument – returns an individual line feature (think of shapefiles here) for each river reach with unique stream order. The mapstruct contains the field tribtoIX and streamorder, such that establishing the connection between tributary and receiving stream is trivial. However, I then realized that this won’t work for connections formed by two streams with the same order, e.g., ’11’, ’22’, etc.

Now, how could a solution to this problem look like. Here is my suggested algorithm, but you may come up with a better (or hopefully simpler) one, if you like. Honestly, I initially thought this was much easier to do. So, here is the function (note that I will include the function on github as soon as I have written a proper help).

function br = sidebranching(S,outputformat)

if nargin==1;
    outputformat = 'nalstr';
    outputformat = validatestring(outputformat,{'nalstr' 'nalnum' 'mapstruct'});

% get confluences
confl = streampoi(S,'confluences','logical') | streampoi(S,'outlets','logical');

% calculate stream order but record at the same the branching at
% confluences
s  = +streampoi(S,'Channelheads','logical');
br = cell(size(s));
br = cellfun(@(x) '', br,'UniformOutput',false);

% calculate stream order
for r = 1:numel(S.ix);
    if s(S.ixc(r)) == 0 || s(S.ixc(r)) <  s(S.ix(r));
        s(S.ixc(r)) = s(S.ix(r));
    elseif s(S.ixc(r)) ==  s(S.ix(r));
        s(S.ixc(r)) = s(S.ix(r))+1;

    if confl(S.ixc(r));
        br{S.ixc(r)} = [num2str(s(S.ix(r))) br{S.ixc(r)}];

br(confl) = cellfun(@(x) sort(x,'ascend'),br(confl),'UniformOutput',false);

for r = 1:numel(S.ix);
    if confl(S.ixc(r));
        s1 = br{S.ixc(r)};
        s1 = str2double(s1(1));
        if s(S.ix(r)) == s1
            br{S.ix(r)} = br{S.ixc(r)};
    elseif confl(S.ix(r));
        br{S.ix(r)} = '';

for r = numel(S.ix):-1:1;
    if isempty(br{S.ix(r)})
        br(S.ix(r)) = br(S.ixc(r));

switch outputformat
    case 'nalnum'
        br = cellfun(@(x) str2double(x),br,'UniformOutput',true);
    case 'mapstruct'
        MS  = STREAMobj2mapstruct(S);
        xy  = zeros(numel(MS),2);
        for r = 1:numel(MS);
            xy(r,1) = MS(r).X(1);
            xy(r,2) = MS(r).Y(1);

        [~,locb] = ismember(xy,[S.x S.y],'rows');
        for r = 1:numel(MS);
            MS(r).sidebranch = br{locb(r)};

        br = MS;

Ok, now let’s prepare some data:

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM,'preprocess','carve');
S = STREAMobj(FD,'minarea',1000);
S = klargestconncomps(S);

We’ll now call above function with the outputformat set to ‘mapstruct’. This will return the same structure array as the function STREAMobj2mapstruct, except that this mapstruct contains another attribute named sidebranch.

MS = sidebranching(S,'mapstruct');

Finally, let’s plot the results. We plot the data with line widths adjusted to streamorder and coloring to the type of side branch.

[unbr,~,loc] = unique({MS.sidebranch});
c = lines(numel(unbr));
hold on
for r = 1:numel(MS);
    h(loc(r)) = plot(MS(r).X,MS(r).Y,'Color',c(loc(r),:),'LineWidth',MS(r).streamorder);
Side branching
Side branching derived for the largest river network in the Big Tujunga DEM.

So, this was definitely not as easy as initially thought. Any ideas how to simplify this are welcome.


Pelletier, J.D., Turcotte, D.L. 2000: Shapes of river networks and leaves: are they statistically similar? Phil. Trans. R. Soc. Lond. B, 355, 307-311.

Strahler, A.N. 1957: Quantitative analysis of watershed geomorphology. Am. Geophys. Union Trans. 38, 913-920.

Tokunaga, E. 1978: Consideration on the composition of drainage networks and their evolution. Geogr. Repl. Tokyo Metro. Univ. 13, 1-27.

Calculating the Width Function

Posted on

Recently I was asked how to calculate the width function of a stream network. My reply was that there is no out-of-the-box solution in TopoToolbox but that it can be easily derived using the functions in TopoToolbox.

First of all, what is the width function of a river network? According to Lashermes and Foufoula-Georgiou (2007), the “width function of a river network is a one-dimensional function which summarizes the two-dimensional branching structure of the river network. It represents the distribution of travel distances through the network and, under the assumption of constant flow velocity, the probability distribution of traveltimes. Thus its significance for understanding the hydrologic response of basins and the scaling characteristics of streamflow hydrographs is important”.

The width is thus nothing but the frequency distribution of distances from the outlet. So, what we need is to calculate the distance from the outlet, which is easily done for stream networks in TopoToolbox because this distance is a property of the STREAMobj.

So, here is a quick solution:

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD  = FLOWobj(DEM,'preprocess','carve');
S = STREAMobj(FD,'minarea',1000);
S = klargestconncomps(S);
d = S.distance;
[n,binc] = hist(d,50);
bar(binc,n,1,'EdgeColor','none','FaceColor',[0.5 0.5 0.5]);
xlabel('x [m]'); ylabel('counts');



Lashermes, B., Foufoula-Georgiou, E., 2007. Area and width functions of river networks: new results on multifractal properties. Water Resources Research 43, W09405–W09405.

Demo on stream network modification

Posted on Updated on

Stream network modification
Stream network modification

A couple of months ago I posted a text on stream network modification which listed some of the ways to modify stream networks in TopoToolbox. These modifications include trimming, extractions, trunk river derivation, etc. Still, many of these functionalities remain somewhat shrouded given that they are all listed in the help section of the function STREAMobj/modify.m where they can be accessed using parameter name-value pairs. To assist users in understanding and finding these tools I have recently added a demo to TopoToolbox that includes some code and a visual output of the different modifications. Hope that this helps you finding the right tools for your analysis.