Avoiding edge effects in chi analysis

Posted on

Analyses that use upslope area usually demand that catchments are completely covered by the DEM. Values of upslope areas may be too low if catchments are cut along DEM edges, and so are estimates of  discharge. How can we avoid including these catchments into our analyses?

In one of my previous posts on chi analysis, I showed a rather long code to include those catchments that have 20% or less of their outlines on the DEM edges. Here, I’ll be more strict. I’ll remove those parts of the stream network that have pixels on DEM edges in their upslope area. By doing this, we make sure that drainage basins are complete which is vital for chi analysis or estimating discharge.

Here is the code:

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM);
S  = STREAMobj(FD,'minarea',1e6,'unit','map');

I  = GRIDobj(DEM,'logical');
I.Z(:,:) = true;
I.Z(2:end-1,2:end-1) = false;

I  = influencemap(FD,I);

The influencemap function takes all edges pixels and derives those pixels that they would influence downstream. Again, this is going to be a logical GRIDobj where true values refer to pixels that might potentially be affected by edge effects. Let’s remove those pixels from the stream network S using the modify function.

S2 = modify(S,'upstreamto',~I);
D  = drainagebasins(FD,S2);
imageschs(DEM,[],'colormap',[1 1 1],'colorbar',false)
[~,x,y] = GRIDobj2polygon(D);
hold on
Edge-unaffected drainage network (blue) and its drainage basins (green). The stream network that is potentially affected is shown as gray solid lines.

Ok, now you have a drainage network that won’t be affected by edge effects for further analysis.


An alternative to along-river swath profiles

Posted on Updated on

Today, I want to show another application of TopoToolbox’s new smoothing suite that Dirk and I have recently released together with our discussion paper in ESURF. I will demonstrate that using different methods associated with STREAMobj enable you to create plots that are quite similar to swath profiles calculated along rivers. Let’s see how it works.

Swath profiles require a straight or curved path defined by a number of nodes. Each node has a certain orientation defined by one or two of its neighboring nodes. Creating swath profiles then involves mapping values from locations that are perpendicular to that orientation and within a specified distance. We will use a simplified version of mapping values lateral to the stream network which is implemented in the function maplateral. maplateral uses the function bwdist that returns a distance transform from all pixels (the target pixels) to a number of seed pixels in a binary image. In our case, seed pixels are the pixels of the stream network and bwdist identifies the  target pixels from which it calculates the euclidean distance to the stream pixels. Usually, there are several target pixels for each seed pixel so that we require an aggregation function that calculates a scalar based on the vector of values that are associated with each seed pixel. Unlike swath profiles, however, this approach entails that some seed pixels may not have target pixels  up to the maximum distance. We can thus be pretty sure that some pixels may be missing some of the information that we want to extract. While we cannot avoid this problem using this approach, we can at least reduce its impact using the nonparametric quantile regression implemented in the crs function.

Here is the approach. We’ll load a DEM, derive flow directions and a stream network which is in this case just one trunk river. We use the function maplateral to map elevations in a distance of 2 km to the stream network. Our mapping uses the maximum function to aggregate the values. Hence, we will have for each pixel along the stream network the maximum elevation in its 2 km nearest-neighborhood. The swath is plotted using the function imageschs and the second output of maplateral. Then we plot the maximum values along with a profile of the stream network.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM);
S  = STREAMobj(FD,'minarea',1e6,'unit','map');
S = klargestconncomps(S);
S = trunk(S);

z = getnal(S,DEM);
[zmax,mask] = maplateral(S,DEM,2000,@max);

imageschs(DEM,mask,'truecolor',[1 0 0],...

hold on
legend('River profile','maximum heights in 2 km distance')


This looks quite ok, but the topographic profile has a lot of scatter and a large number of values seem not to be representative. Most likely, those are the pixels whose closest pixels fail to extend to the maximum distance of 2 km. Rather, I’d expect a profile that runs along the maximum values of the zigzag line. How can we obtain this line? Well, the crs function could handle this. Though it was originally intended to be applied to longitudinal river profiles, it can also handle other data measured or calculated along stream networks. The only thing we need to “switch off” is the downstream minimum gradient that the function uses by default to create smooth profiles that monotoneously decrease in downstream direction. This is simply done by setting the ‘mingradient’-option to nan. I use a smoothing parameter value K=5 and set tau=0.99 which forces the profile to run along the 99th percentile of the data. You can experiment with different values of K and tau, if you like.

zmaxs = crs(S,zmax,'mingradient',nan,'K',5,'tau',0.99,'split',0);
hold on
plotdzshaded(S,[zmaxs z],'FaceColor',[.6 .6 .6])
hold off


This looks much better and gives an impression on the steepness of the terrain adjacent to the stream network. Let’s finally compare this to a SWATHobj as obtained by the function STREAMobj2SWATHobj.

SW = STREAMobj2SWATHobj(S,DEM,'width',4000);
xlabel('Distance upstream [m]');
ylabel('Elevation [m]');


Our previous results should be similar to the maximum line shown in the SWATHobj derived profile. And I think they are. In fact, the SWATHobj derived profiles also shows some scatter that is likely due to the sharp changes of orientation of the swath centerline.  While we can remove some of this scatter by smoothing the SWATHobj centerline, I think that the combination of maplateral and the CRS algorithm provides a convenient approach to sketch along-river swath profiles.

Steepness derived from smoothed river profiles

Posted on

Today, I gave a talk in the weekly seminar of Manfred Strecker’s group at our institute. The talk had the title: Introduction to TopoToolbox and its application in tectonic geomorphology. Among other, I provided an example of how to calculate and plot the steepness of river networks (see the new function ksn) and how the new smoothing function crs (Schwanghart and Scherler 2017) is helpful here. I am providing the code here for the DEM of the Big Tujunga catchment, but you can easily adapt the code to your own needs.

First, load the DEM and derive flow directions and the stream network. In addition, calculate a flow accumulation grid that we need for calculating ksn.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM,'preprocess','carve');
A = flowacc(FD);

We then extract the stream network for a threshold upstream area of 1 sqkm.

S = STREAMobj(FD,'minarea',1e6,'unit','map');
S = klargestconncomps(S,1);

To hydrologically correct and smooth the river profile, we use the function crs:

zs = crs(S,DEM,'K',6,'tau',0.1);

And finally, we plot it all

imageschs(DEM,[],'colormap',[1 1 1],'colorbar',false,'ticklabels','nice');
hold on
hx = colorbar;
hx.Label.String = 'k_{sn}';
Map of ksn values in the Big Tujunga catchment.

Which parameter values should be used to smooth the stream profile? The values depend on the amount of scatter in the profiles, their spatial resolution, and the amount of smoothing. Visual cross-checking of the results is thus advisable. In this lecture, I presented the tool crsapp that let’s you identify a suitable set of parameter values. A screenshot of the tool can be seen below.

The crsapp GUI let’s you interactively choose the right parameters for smoothing river profiles.


Schwanghart, W., Scherler, D., 2017. Bumps in river profiles: the good, the bad, and the ugly. Earth Surface Dynamics Discussions 1–30. [DOI: 10.5194/esurf-2017-50]

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…

Chimaps in a few lines of code (5)

Posted on Updated on

Having prepared a stream network and equipped with a reasonable value of the m/n ratio, we are now ready to plot a chimap that visualizes the planform patterns of chi. The main interest in these maps lies in chi values near catchment divides as large differences between adjacent catchment would indicate a transient behavior of drainage basin reorganization.

Some of you might have already experimented with TopoToolbox to create chimaps. Perhaps you became exasperated with the function chiplot that is restricted to calculations with only one drainage basin and has a bewildering structure array as output. The reason for the confusing output of chiplot is that it is fairly old. At this time, I hadn’t implemented node-attribute lists that are now more common with STREAMobj methods.

Realizing this shortcoming of chiplot, I wrote the function chitransform. chitransform is what I’d refer to as a low-level function that solves the chi-equation using upstream cumulative trapezoidal integration (see the function cumtrapz). chitransform requires a STREAMobj and a flow accumulation grid as input and optionally a mn-ratio (default is 0.45) and a reference area (default is 1 sqkm). It returns a node-attribute list, i.e., a vector with chi values for each node in the STREAMobj. Node-attribute lists are intrinsically tied to the STREAMobj from which they were derived. Yet, they can be used together with several other TopoToolbox functions to produce output.

Ok, let’s derive chi values for our stream network:

A = flowacc(FD); % calculate flow accumulation
c = chitransform(S,A,'mn',0.4776);

In the next step, we will plot a color stream network on a grayscale hillshade:

imageschs(DEM,[],'colormap',[1 1 1],'colorbar',false,'ticklabel','nice');
hold on
hold off
Chimap of the Mendocino Triple Junction.

Interestingly, there seem to be some locations with quite some differences in chi values on either side of the divide. “Victims” seem to be rather elongated catchments draining northwest. Let’s zoom into one of these locations.

Detail of the chimap with arrow indicating the expected direction of divide migration.

Are these significant differences? Well, it seems by just looking at the range of values. However, to my knowledge no approach exists that provides a more objective way of evaluating the significance of contrasting chi values and their implications about rates of divide migration. Still, we now have a nice map that can aid our geomorphic assessment together with the tectonic and geodynamic interpretation of the Mendocino Triple Junction.

Unfortunately, I must leave the discussion to you since I am quite unfamiliar with the region. If anyone wants to share his or her interpretation, I’d be more than happy to provide space here. So far, I hope that these few posts on chimaps were useful to you and sufficiently informative to enable you to compute chimaps by yourself. In my next post, I will give a short summary and show with another example that eventually chimaps can be derived really in a few lines of code.

Chimaps in a few lines of code (4)

Posted on Updated on

My last post described the math behind chi analysis. By using the integral approach to the stream power model of incision, we derived an equation that allows us to model the longitudinal river profile as a function of chi. At the end of the last post I stated that this model is a straight line if we choose the right m/n ratio, that is the concavity index. Today, I’ll show how we can obtain this m/n ratio by means of nonlinear optimization using the function chiplot.

Ideally, we find the m/n ratio using a perfectly graded stream profile that is in steady state. I scanned through a number of streams using the GUI flowpathapp and found a nice one that might correspond to Cooskie Creek that Perron and Royden (2013) used in their analysis.

The graphical user interface flowpathapp allows visual exploration of planform river patterns and corresponding profiles.

The menu Export > Export to workspace allows saving the extracted STREAMobj St to the workspace. Then, I use the function chiplot to see how different values of the m/n ratio affect the transformation of the river profile. Setting the ‘mnplot’ option makes chiplot return a figure with chi-transformed profiles for m/n ratios ranging from 0.1 to 0.9 in steps of 0.1 width.

chi-transformed river profiles for different values of the mn-ratio.

Clearly, the chi-transformed profile varies from concave upward for low m/n values to convex for high m/n values. A value of 0.4-0.5 seems to be most appropriate but choosing a value would be somewhat hand-wavy. Thus, let’s call the function again, but this time without any additional arguments.

C = chiplot(S,DEM,A);
chi-transformed profile for an optimal value of the m/n-ratio.

Calling the function this way will prompt it to find an optimal value of m/n. In addition, it will return a figure with the linearized profile. By default, the function uses the optimization function fminsearch that will vary m/n until it finds a value that maximizes the linear correlation coefficient between the chi-transformed profile and a straight line. So what is that value? Let’s look at the output variable C:

C = 

          mn: 0.4776
        mnci: []
        beta: 0.1565
      betase: 2.2813e-04
          a0: 1000000
          ks: 114.8828
          R2: 0.9988
       x_nal: [581x1 double]
       y_nal: [581x1 double]
     chi_nal: [581x1 double]
       d_nal: [581x1 double]
           x: [582x1 double]
           y: [582x1 double]
         chi: [582x1 double]
        elev: [582x1 double]
      elevbl: [582x1 double]
    distance: [582x1 double]
        pred: [582x1 double]
        area: [582x1 double] 

C is a structure array with a large number of fields. is the value we are interested in. It is 0.4776 and corresponds nicely to previously reported m/n ratios although it differs from the value 0.36 found by Perron and Royden (2013). Of course, this value might differ from river to river and you should repeat the analysis for other reaches to obtain an idea about the variability of m/n.

Ok, we have the m/n ratio now. In my next blog we will eventually produce the chimap that we initially set out for.

P.S.: We could have also used slope-area plots to derive the m/n ratio (see the function slopearea). However, along-river gradients are usually much more noisy than the profiles and power-law fitting a delicate matter. I’d definitely prefer chi-analysis over slope-area analysis.

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.