chi analysis

Working with knickpointfinder

Posted on Updated on

I have learnt from numerous comments here on this blog, that the knickpointfinder function seems to be quite popular. However, working with the output of the function often entails some difficulties. So here is a short post about knickpointfinder.

Let’s load and prepare some data first:

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

In these first four lines of code, I loaded the DEM of the Big Tujunga catchment, derived flow directions (stored in FD, an instance of FLOWobj), and then extracted the stream network S based on a minimum supporting drainage area of 1000 pixels which is about 0.9 sqkm. In the last line, I clipped the stream network to the largest connected stream network, that is in this case the largest drainage basin.

Now we are ready to run knickpointfinder:

[zs,kp] = knickpointfinder(S,DEM,'tol',30,'split',false);
Animated output as displayed by the knickpointfinder function if ‘plot’ is set to true.

The function takes the stream network S, the digital elevation model DEM and then a few parameter name-value pairs. knickpointfinder reconstructs to the actual longitudinal profile by iteratively removing a curvature constraint at those locations that have a maximum vertical offset to the actual profile. After few iterations, the vertical offset decreases and ‘tol’ sets the minimum vertical offset where the algorithm should stop. The value of tol should reflect uncertainties that are inherent in longitudinal river profile data (Schwanghart and Scherler, 2017).

‘split’ enables running the function in parallel. The function splits the stream network into connected components (see the function STREAMobj2cell) and sends each to individual workers. Parallel execution, however, is only faster if there are several drainage basins. Since we have only one in this example, I set ‘split’ to false. There are other options to control verbose and visual output and you can explore these options by typing

help knickpointfinder

Now, knickpointfinder returned two output variables: zs and kp. zs is a node-attribute list that contains the reconstructed river profile. You can plot it using plotdz.

hold on
hold off

kp is a structure array that holds all kind of information.


kp = 

  struct with fields:

           n: 26
           x: [26×1 double]
           y: [26×1 double]
    distance: [26×1 double]
           z: [26×1 double]
      IXgrid: [26×1 double]
       order: [26×1 double]
          dz: [26×1 double]
         nal: [6985×1 logical]

Admittedly, this is a bit messy and there will be some changes in the future. But this will still have to wait a bit. Now what do the structure fields contain? kp.n is 26. Hence, the algorithm has identified 26 knickpoints. kp.x and kp.y contain the coordinates and kp.z  the elevations of the knickpoints. The distance from the outlet is stored in kp.distance. kp.IXgrid contains the linear indices into the DEM. TopoToolbox works a lot with linear indices and thus this output is very crucial. kp.order contains the order with which the knickpoints where identified. However, several knickpoints may have the same order and thus this information is not very meaningful. is the vertical offset between the modelled profile at the time when the knickpoint was identified. Finally, kp.nal is a logical vector and a node-attribute list. It allows indexing into other node attribute lists.

Now, here are a few things you might want to do with the structure array kp.

1. Plot a map

hold on
hold off

2. Plot a longitudinal river profile

hold on
hold off

3. Plot a longitudinal river profile in chi space

A = flowacc(FD);
c = chitransform(S,A,'mn',0.45);
hold on
[~,locb] = ismember(kp.IXgrid,S.IXgrid);
ckp = c(locb);
hold off
xlabel('\chi [m]')

4. Export as shapefile

MS = STREAMobj2mapstruct(S);
MKP = struct('X',num2cell(kp.x),'Y',num2cell(kp.y),'Geometry','Point');

Ok, that’s a basic work flow that you may apply easily to your own DEM. Give it a try.


Schwanghart, W., Scherler, D., 2017. Bumps in river profiles: uncertainty assessment and smoothing using quantile regression techniques. Earth Surface Dynamics, 5, 821-839. [DOI: 10.5194/esurf-5-821-2017]


Some new functions and features in TopoToolbox

Posted on

Let me highlight some functions that I have added recently. I think that some of these new features come unnoticed.

First, there is the function nal2nal. Obviously, this function handles “nals” which are node-attribute lists. Node-attribute lists are column vectors where each element provides some value (e.g. elevation or distance) for every node in the stream network. Thus, nals are directly linked to the network topology of a particular instance of STREAMobj. But what happens if you want to use a nal of one stream network and map it to another one. Well, that can be done with nal2nal.

Here is an example that uses nal2nal. We map chi-values measured along an entire river network to its first-order streams.

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

A = flowacc(FD);
c = chitransform(S,A);
S2 = modify(S,'streamorder',1);
c2 = nal2nal(S2,S,c,0);
hold on
h = colorbar;
h.Label.String ='\chi';
chi values of first-order streams

Another new function is mapfromnal. The function maps values measured at each node of the stream network to its upstream hillslope pixels. I showed an example in my previous post on chi inferiority.

Then, there is stackedplotdz. This function takes a STREAMobj and several node attribute lists and plots them on top of each other. If you have MATLAB 2018b, then the function comes with additional functionalities provided by the function stackedplot.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM);
S = STREAMobj(FD,'minarea',1000);
S = klargestconncomps(S);
S = trunk(S);
z = getnal(S,DEM);
zs = crs(S,DEM,'tau',0.1,'K',10,'split',false);
g = gradient(S,zs);
k = ksn(S,zs,flowacc(FD));
stackedplotdz(S,{z zs g k},…
'ylabels',{'Elevation' 'smoothed elevation' 'gradient' 'ksn'});
Multiple variables plotted on top of each other using stackedplotdz

Finally, I have created a new repository for some example DEMs. Maybe I have seen Big Tujunga too often now. The repository contains relatively small DEMs (about 30 Mb, github doesn’t allow much larger files) that I will use for demonstration purposes in the future. Until now, there are only two DEMs, but their number will increase in the future. The repository can be accessed via the function readexample that you’ll find the IOtools directory.

DEM = readexample('taiwan');

How do you keep up to date with new releases of TopoToolbox? Do you use git? If not, I’d recommend it.

X inferiority

Posted on

I just need to finish the year with a quick blog entry… And this one is – once again – about chi maps.

I thought about different ways to display chi maps. chi maps are usually plotted using stream networks. However, to this end, chi maps are not so much about stream networks but rather about drainage divides and how they potentially migrate in the future. In addition, the distance of potential divide migration depends on the relative difference of chi values on either side of the divide, not so much on the absolute values. How could we create a map that reflects this information

Here is an idea: First, we create a chi transform of the river network (see function chitransform). Then we map these chi-values to the entire grid. This can be accomplished with the new function mapfromnal which lets you map values from a stream network to its nearest upstream contributing pixels. Then, we will derive some-pixel wide fringes of chi values around the divides and calculate their spatial average. Dividing the chi values by these spatial averages then quantifies the degree to which some part of the divide is either a victim or an aggressor. Since victims have higher values, I call this metric the chi-inferiority.

Now here is the code:

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD  = FLOWobj(DEM);
S   = STREAMobj(FD,'minarea',500);
S   = removeedgeeffects(S,FD);
[~,zb] = zerobaselevel(S,DEM);
A   = flowacc(FD);

c   = zb + chitransform(S,A);
C   = mapfromnal(FD,S,c);
D   = drainagebasins(FD,S);

for r = 1:max(D)
    I = D==r;
    I = erode(I,ones(21)) & I;
    C.Z(I.Z) = nan;

C1 = dilate(C,ones(21));
C2 = erode(C,ones(21));

C  = C/((C1 + C2)/2);
hold on
hold off
Map of chi-inferiority of the Big Tujunga catchment. Coloring is based on the vik colormap of the scientific colormaps of Fabio Crameri.

Happy New Year everone!

Bayesian Optimization of the mn-ratio

Posted on Updated on

River profiles are concave upward if they are in a dynamic equilibrium between uplift and incision, and if our simplified assumptions of steady uplift and the stream power incision law (SPL) hold. The concavity derives from the SPL which states that along-river gradients S are proportional to upslope area A exponentiated by the negative mn-ratio.

I have mentioned the mn-ratio several times in this blog. Usually, we calculate it using slope-area plots or chi analysis both of which are included in TopoToolbox. However, these methods usually lack consistent ways to express the uncertainties of the mn-ratio. The lack of consistency is due to fitting autocorrelated data which elude a straightforward statistical analysis.

Today, I want to present a new function that uses Bayesian Optimization with cross-validation to find a suitable mn-ratio. While Bayesian Optimization is designed to find optimal values of objective functions involving numerous variables, solving an optimization problem with mn as the only variable nicely illustrates the approach.

Bayesian Optimization finds a minimum of a scalar-valued function in a bounded domain. In a classification problem, this value could be the classification loss, i.e., the price paid for misclassifications. In a regression problem, this value might refer to the sum of squared residuals. The value might also be derived using cross-validation,  a common approach to assess the predictive performance of a model. Such cross-validation approaches might take into account only random subsets of data, which entails that the value to be optimized might not be the same for the same set of input parameters. Bayesian Optimization can handle stochastic functions.

Now how can we apply Bayesian Optimization for finding the right mn-ratio? The new function mnoptim uses chi analysis to linearize long-river profiles. If there are several river catchments (or drainage network trees), the function will pick a random subset of these trees to fit a mn-ratio and then tests it with another set of drainage basins. This allows us to assess how well an mn-ratio derived in one catchment can actually be applied to another catchment. The goal is to derive a mn-ratio that applies best to other catchments.

Now let’s try this using the new function mnoptim. Here is the code that I’ll use for entire SRTM-3 DEM of Taiwan. I’ll clip the stream network to the area upstream of the 300 m contour to avoid an influence of the alluvial low-lying reaches.

DEM = GRIDobj('taiwan.tif');
FD = FLOWobj(DEM);
S  = STREAMobj(FD,'minarea',1e6,'unit','map');
C  = griddedcontour(DEM,[300 300],true);
S  = modify(S,'upstreamto',C);
A  = flowacc(FD);
[mn,results] = mnoptim(S,DEM,A,'optvar','mn','crossval',true);
% we can refine the results if we need
results = resume(mn);
% and get an optimal value of mn:

ans =




Estimated prediction losses for different values of the mn-ratio.

Now this nicely derives the optimal mn-value of 0.415 which is close to the often reported value of 0.45. Moreover, based on the plot, we gain an impression of the uncertainty of this value. In a transient landscape with frequent knickpoints, the uncertainty about the mn-ratio will be probably larger.

Note that mnoptim requires the latest version of MATLAB: 2017b, as well as the Statistics and Machine Learning Toolbox. It also runs with 2017a, but won’t be able to use parallel computing then.

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.

Chimaps in a few lines of code (final)

Posted on Updated on

A couple of weeks ago, I covered a series of blog posts about chi maps. This post links to these posts for future reference.

  • Part 1: Overview and approach to creating chimaps
  • Part 2: Preparing the stream network
  • Part 3: The mathematical basis
  • Part 4: Finding the right m/n ratio
  • Part 5: Plotting the chi-map

In addition, I like to draw your attention to the tool ChiProfiler that has been developed by Sean Gallen from the ETH Zürich. ChiProfiler is based on the TopoToolbox functions and provides a visual interface to stream profile analysis and includes ksn and chi analysis, amongst other. 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.