Crossing divides

Posted on Updated on

Divide mobility and changes in flow network configuration are hot topics these winter days. How mobile are drainage divides, what controls their mobility, which are the involved time scales, and how do shifting divides affect other metrics derived from digital elevation models?

One of the issues addressed by some studies is which metrics are useful to characterize divide movements (Forte and Whipple, 2018, Scherler and Schwanghart, 2019). For example, there are chi maps which allow us to map spatial imbalances in drainage network configuration. However, whether these imbalances actually translate into divide movements, or whether these movements will occur in some distant future remains unclear. Gilbert metrics in turn quantify cross-divide differences in hillslope gradient (or other similar metrics) and thus may provide a more suitable proxy for processes that actually act along or close to divides. To this end, it may often be useful to look at both types of metrics.

Now, you know that I like visualisations. And today, I want to briefly show how to combine hillslope and divide geometry with longitudinal river profiles of two adjacent rivers. I am hoping to provide the computational means to automatically derive such profiles which might prove very insightful.

The example uses the Big Tujunga data which we first load. We derive flow directions and the stream network. We will only use the trunk rivers.

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

In a next step, I am extracting two rivers which are sourced in the so-called Chilao Flats and drain this low-relief relict landscape. Note that extractconncomps is an interactive function and you will need to manually pick the trunk rivers yourself.


Plotting the two trunk rivers as profiles (plotdz) generates two profiles that have their outlets set to the distance value 0. That’s not so intuitive and I actually want that the channelheads should be opposing each other.

So, here is a little trick that I can use to connect the two rivers and thereby generate a path across the divide. First, I convert the STREAMobj of the two rivers to a GRIDobj. Then, I run a gray-weighted distance transform to calculate the least-cost path emanating from the outlet of one of the rivers to all other pixels in the DEM. The costs of moving from one pixel to another is elevation, but I decrease the costs at the rivers themselves, which forces the algorithm to run along the river network.

GS = STREAMobj2GRIDobj(S);
DEMaux = DEM;
DEMaux.Z(GS.Z) = 1;
IX = streampoi(S,'outlet','ix');
DEMaux.Z = graydist(DEMaux.Z,IX(1));


I now have a new GRIDobj DEMaux. I chose this variable name because DEMaux represents some auxiliary topography (the same approach is actually used to route through flat sections or closed pits in the DEM).

Now I am using this auxiliary topgraphy to calculate a FLOWobj, and a STREAMobj sourced at the outlet of the second channel. The resulting stream starts at the outlet of the second channel, runs upstream, crosses the divide, and then follows the course of the first channel until its outlet.

FDaux = FLOWobj(DEMaux);
Saux = STREAMobj(FDaux,'channelhead',IX(2));

% plot it
imageschs(DEM,[],'colormap',[1 1 1],'colorbar',false)
hold on


Ok, that’s quite neat. Now, in addition, I might want to plot the river and cross-divide profile together with chi-values.

A = flowacc(FD);
c = chitransform(S,A);
[~,zb] = zerobaselevel(S,DEM);
c = c + zb;

C = GRIDobj(DEM);
C.Z(S.IXgrid) = c;
caux = getnal(Saux,C);

hold on
h = colorbar;
h.Label.String = '\chi [m]';
hold off


Ok, that’s what I wanted. Have fun coding!


Forte AM, Whipple KX. 2018. Criteria and tools for determining drainage divide stability. Earth and Planetary Science Letters 493 : 102–117. DOI: 10.1016/j.epsl.2018.04.026
Scherler D, Schwanghart W. 2019. Identification and ordering of drainage divides in digital elevation models. Earth Surface Dynamics Discussions : in review (open discussion esurf-2019-51).

CASCADE Toolbox – A toolbox for river sediment connectivity assessment

Posted on

I am glad to introduce today’s guest bloggers Marco Tangi and Rafael Schmitt, authors of the MATLAB and TopoToolbox based model CASCADE. Here they briefly describe the model and its implementation.

The CASCADE (CAtchment Sediment Connectivity And DElivery) model is a network model for quantifying and visualizing sediment connectivity. CASCADE was developed for sediment management and geomorphic studies in large river systems  (Schmitt et al., 2016). An improved and exhaustively documented MATLAB implementation of the CASCADE is now publicly available with the CASCADE toolbox (Tangi et al., 2019).

CASCADE represents the transport of sediment from many individual sources of sediment in a river network as an individual cascading process for each source. Sediment transport along each of the individual sediment cascades is governed by local hydromorphic conditions in the downstream river network, the grain size supplied from the source, and the grainsize distributions in each downstream reach. In the CASCADE toolbox, users can choose from common empirical functions for fractional sediment transport. The graph-based tracking of individual cascades allows identifying some key parameters of fluvial sediment connectivity, e.g., the fate of sediment originating from a specific source, as well as sediment provenance, i.e., from which sources sediment is delivered to a downstream reach.

CASCADE is a statistical model designed for numerical efficiency, which enables application to large river networks and stochastic simulations of sediment connectivity (Schmitt et al., 2018a). Input data requirements are flexible to allow using gauged data as well as outputs from e.g., hydrologic models or similar. River network extraction and geomorphic processing from a DEM is mostly performed using Topotoolbox routines, while most other hydro-geomorphological features can be obtained from remotely sensed data. CASCADE has been applied for quantifying basin-scale connectivity in poorly monitored rivers (Schmitt et al., 2018b), evaluating human impacts on sediment budgets (Schmitt et al., 2018a), and basin scale infrastructure planning (Schmitt et al., 2019). By design, CASCADE can be used to represent human interventions on sediment transport (e.g., dams, sediment mining and reintroduction, changes in sediment supply).

Here, we present some example of how the model is applied for the Vjosa River in the Southern Balkan:

A) Extraction of the river network, reach segmentation and assigning hydromorphologic attributes

Sediment routing in the CASCADE model is based on hydro-morphologic attributes and common sediment transport fomulas for mixed-sized sediments. Outputs include total sediment flux as well as deposition and potential entrainment in each reach. Note that deposition and entrainment can be occurring at the same time in one reach for different grainsizes.

B) Results of a baseline run can then be used to interactively evaluate impacts of river management on network-scale sediment transport. E.g., for evaluating impacts of different combinations of planned dams on sediment transport in the Vjosa River.

The CASCADE toolbox is available on the website ( and the GIThub repository ( The accessible material includes a user manual, fully commented code, ready-to-use scripts and example workspaces, as well as a growing list of case studies and use cases.



Schmitt, R.J.P., Bizzi, S., Castelletti, A., 2016. Tracking multiple sediment cascades at the river network scale identifies controls and emerging patterns of sediment connectivity. Water Resour. Res. 3941–3965. [10.1002/2015WR018097]

Schmitt, R.J.P., Bizzi, S., Castelletti, A., Kondolf, G.M., 2018a. Improved trade-offs of hydropower and sand connectivity by strategic dam planning in the Mekong. Nature Sustainability 1, 96–104. [10.1038/s41893-018-0022-3]

Schmitt, R.J.P., Bizzi, S., Castelletti, A.F., Kondolf, G.M., 2018b. Stochastic Modeling of Sediment Connectivity for Reconstructing Sand Fluxes and Origins in the Unmonitored Se Kong, Se San, and Sre Pok Tributaries of the Mekong River. Journal of Geophysical Research: Earth Surface. [10.1002/2016JF004105]

Schmitt, R.J.P., Bizzi, S., Castelletti, A.F., Opperman, J.J., Kondolf, G.M., 2019. Planning dam portfolios for low sediment trapping shows limits for sustainable hydropower in the Mekong. Science Advances, 5, 10, eaaw2175. [10.1126/sciadv.aaw2175]

Tangi, M., Schmitt, R., Bizzi, S., Castelletti, A., 2019. The CASCADE toolbox for analyzing river sediment connectivity and management. Environmental Modelling & Software 119, 400–406. [10.1016/j.envsoft.2019.07.008]

Colored stream networks in Google Earth

Posted on

In my previous posts (here and here), I showed that MATLAB and TopoToolbox feature tools for plotting stream networks on web maps. However, while wmplot – the TopoToolbox function – neatly serves this purpose, it doesn’t easily allow you to share the data easily so that your peers can either plot it themselves, or that the data can be embedded into online maps.

KML (keyhole markup language) is a file format that can be easily exchanged, viewed in Earth browsers such as Google Earth or ArcGIS Earth and imported into any GIS. Some academic publishers also encourage authors to provide kml files to enrich their papers with online maps. Tools to convert stream networks to kml files would thus be nice to have.

MATLAB’s Mapping Toolbox lets you export kml files. Here is an example that takes the Taiwan SRTM (note that you need an internet connection for these examples):

DEM = readexample('taiwan');
DEM = inpaintnans(DEM);
FD  = FLOWobj(DEM);
A   = flowacc(FD);
S   = STREAMobj(FD,A > 1000)

[lat, lon] = STREAMobj2latlon(S);

The function STREAMobj2latlon converts a STREAMobj to the nan-punctuated coordinate vectors lat and lon. Nan-punctuated means that individual lines are separated by nans in these vectors. kmlwriteline then lets you write these lines to a kml-file.

However, available MATLAB kml-utilities such as kmlwriteline have limited functionality. These limitations particularly pertain to writing files that combine different geometries, colors and symbols. Gladly, solutions exist, and here we adopt the tools provided by the kml-toolbox which you should download and install before running the next code.

In my last post, I asked if someone could help with linking TopoToolbox to the kml-toolbox, and Kai Hu, PhD student in Geological Sciences at the University of Wisconsin-Madison, came up with a nice solution:

c   = chitransform(S,A,'a0',1);
[~,zb] = zerobaselevel(S,DEM);
c   = c+zb;
[lat, lon, c_latlon]=STREAMobj2latlon(S,c);
%convert the nodes of stream networks into lats and lons and their
%corresponding Chi values at each of the nodes.

fh = figure;
ax = gca;
z  = lon*0;
h  = surface([lon lon],[lat lat],[z z],[c_latlon c_latlon],...
        'linewidth',2);%This surface plot is drawn from plotc function.

k = kml('Stream network colored by dChi')

You should have a kmz-file named ‘Stream network colored by dChi.kmz’ in your current folder by now. You can open the file using Google Earth.

Kai’s solution is based on the trick used in the function STREAMobj/plotc which applies the function surface to plot colored stream networks. The function k.transfer than transfers the figure axes to a kml file as an image overlay. This is actually pretty cool, because you can prepare your plot in MATLAB including all kinds of features (e.g. points or contours) (note that you must use axes with geographic coordinates) and then view it in Google Earth.

The drawback of Kai’s solution is that the stream network is exported as a rasterized image which may not look nice at different zoom levels. Moreover, you cannot access, modify, or query the displayed features. To overcome these problems, I now wrote the function STREAMobj2kml. Internally, the function uses STREAMobj2mapstruct, which generates individual line segments with attribute values. Optional parameters let you define the length of these segments as well as the attribute (either a node-attribute list or GRIDobj). In addition, there is an aggregation function that summarizes values available for each stream node to a scalar value for each stream segment. Here is how you call the function:


Note, however, that this is just a start to make great web maps. If you call STREAMobj2kml with one output argument, the function will return a kml-object that can be complemented by additional visualizations such as contours, text, or scatter plots. I’ll certainly show an example in one of my future posts.

Colored stream networks on web maps

Posted on

In my last blog entry, I showed that the TopoToolbox function wmplot makes it easy to plot stream networks on web maps. In a comment to this post, KH asked whether there is a way to color the networks as you would do using the function STREAMobj/plotc.

I searched the documentation of wmline, the function wrapped by wmplot, but apparently, you can only choose one color per line. This means that we need to cut the network into pieces and separately plot these pieces with individual colors. Gladly, there is already a function to support this approach: STREAMobj2mapstruct.

However, iteratively plotting individual lines in a web map takes time. I don’t know the cause of the problem, but plotting a line takes long, irrespective its size. A seizable solution is to merge lines with similar feature attributes. This can be done nicely with the geoshape and mapshape objects. These are pretty handy objects and I think they will become even more powerful in the future if Mathworks and the Mapping Toolbox add functionalities that better support coordinate systems and transformations. As I will likely use these objects more often from now on, I have written a new function called STREAMobj2shape.

Calling STREAMobj2shape, however, is a bit tricky. The syntax can be quite long and I know that many users are rather reluctant to read help files. Thus, I also edited wmplot so that it takes additional arguments that control the coloring.

wmplot now makes it very easy. Here is how TopoToolbox enables plotting a chi map of Taiwan in a web map just in a few lines of code. Plotting takes a while, so be patient and have fun watching the rather unwanted animation of headward growing river networks.

DEM = readexample('taiwan');
DEM = inpaintnans(DEM);
FD  = FLOWobj(DEM);
S   = STREAMobj(FD,'minarea',1000);
A   = flowacc(FD);
c   = chitransform(S,A);
Chi-colored stream network of Taiwan.

KH also came up with the idea that these web maps should be sharable with others. Unfortunately, I think that this is currently not possible. If you have ideas, please let me know.

I am currently also exploring ways to make colored kml or kmz files. I came across the kml toolbox, which comes in very handy. If anyone is interested in co-developing such a tool, that’ll be great.

Stream networks on web maps

Posted on

Did you know that TopoToolbox has a function that lets you plot stream networks on the MATLAB web map? Did you even know that you can access web maps with MATLAB?

Since 2016a, MATLAB has the function webmap. The function opens a web map browser that provides access to a list of available base maps such as World Street Map, USGS imagery, and many more.

TopoToolbox uses this functionality in the function wmplot (web map plot). I use this function quite frequently because it allows me to compare derived stream networks with the river courses visible in satellite imagery.

In the inside, the function wmplot is rather terse. The function reprojects the coordinates back to geographic coordinates (see function STREAMobj2latlon, (the STREAMobj must have a coordinate system!)) and then plots the nan-separated vectors using wmline.

Now here is an example:

DEM = readexample('taiwan');
DEM = inpaintnans(DEM);
FD = FLOWobj(DEM);
S = STREAMobj(FD,'minarea',10000);

However, how do you plot drainage basins? There is no extra function to do this, but a few lines of code will do.

D = drainagebasins(FD,S);
MS = GRIDobj2polygon(D);
[lat,lon] = cellfun(@(x,y) minvtran(S.georef.mstruct,x,y),{MS.X},{MS.Y},'UniformOutput',false);
[MS.Lat] = deal(lat{:});
[MS.Lon] = deal(lon{:});
G = geoshape(MS);
wmline(G,'color',[.5 .5 .5])

That’s it!

TerraceM-2 – MATLAB software for the analysis und modelling of coastal and lake terraces

Posted on Updated on

Todays blog post is written by my colleage Julius Jara-Muñoz (University of Potsdam). Julius is the author of TerraceM, a software for the analysis and modelling of marine and lacustrine terraces.

TerraceM-2 is a Matlab® graphic-user interface to analyze and model marine and lacustrine terraces. TerraceM was originally designed to obtain accurate morphometric assessments of terrace morphology. The new TerraceM-2 version includes novel routines to map the elevation and spatial distribution of terraces, to model their formation and evolution, and to estimate fault-slip rates from terrace deformation patterns. TerraceM-2 uses nested Topotoolbox functions (Schwanghart and Scherler, 2014) to perform basic analyses, such as swath profile extraction, calculation between topographic parameters, and DEM visualization, significantly improving the processing speed and mapping capabilities of TerraceM-2.

TerraceM-2 comprises three modules, a) Mapping tools (MAP) b) Landscape evolution model (LEM), and c) Dislocation Model (DIM), designed to be used both as independent and complementary GUI’s. The MAP module is designed for mapping marine and lacustrine terraces using swath profiles and surface classification models. The LEM module is designed to model terraces formation simulating the effect of coastal erosion and vertical deformation, allowing comparing modelled and measured terraces in order to determine uplift and erosion rates. The DIM module uses the spatial distribution of terrace elevations or surface uplift, comparing them with fault dislocation models to determine fault slip rates.

1. Here is an example of marine terrace mapping using the TerraceM-2 MAP GUI

2. Here is an example of shoreline angle mapping using the TerraceM_staircase function.

%Add path to topotoolbox folder, TerraceM-2 Maptools, and dem and shapefile.
addpath /Users/TerraceM_2018/Stations/HUEVO
addpath /Users/TerraceM_2018/topotoolbox-master
addpath /Users/TerraceM_2018/TerraceM_Maptools

Box=shaperead('Huevo_clip.shp'); %Shapefile containing the profiles boxes 
DEM=GRIDobj('DEM_huevo.tif'); %DEM generated using Topotoolbox
profnum=2; %Profile number to process

%Create swath function

%Display swath profile 
xlabel('Distance along profile (m)'); ylabel('Elevation (m)');
legend('Min elevation','Mean elevation','Max elevation')

%Mapping a single shoreline angle
Swath profile generated using the TerraceM_makeswath function and mapped using the TerraceM_staircase function.

3. Here is another example of how to create a LEM movie using the TerraceM_anderson function.

%add path to TerraceM-2 LEM and the sealevel folders
addpath /Users/TerraceM_2018/TerraceM_LEM
addpath /Users//TerraceM_2018/TerraceM_LEM/sealevel
sl=load('Bintanja_total.mat');%load a sea-level curve
%parameters in structure format
param.tmax=500;%max time (ka)
param.tmin=0;%min time (ka)
param.dt=100;%bin size of time intervals (yrs)
param.dx=5;%bin size of x axis (m)
param.sill=150;%m, avoid sill by adding an extremely deep value
param.slope_shelf=10;%initial slope model (deg)
param.uplift_rate=1.8; %m/ka
param.initial_erosion=0.3;% initial erosion rate or E0 (m/yr)
param.cliff_diffusion=0.02;% cliff slope difussion (m^2/yr)
param.wave_height=5;%wave height (m)
param.vxx=0;%manually extend x axis if model run out of bounds
param.smth=40;%smooth factor of sea-level curve (yr);%video mode activated
%optional input parameters, ignore for quick video display
param.cumula=0;%cummulative plot deactivated
param.snap=90;%number of video frames
%create axes for video
%call LEM function

The integration of TerraceM-2 modules and routines, in addition to user-friendly GUIs, optimizes the time required for mapping and modelling marine and lacustrine terraces for inexperienced users in Matlab® scripting. Furthermore, the TerraceM-2 architecture comprises independent functions allowing experienced users creating custom scripts and workflows beyond the GUI environment.

TerraceM-2 is available at, there you can also find tutorials and examples, updates, and a list of recent publications using TerraceM.


Schwanghart, W., and Scherler, D. (2014). Topotoolbox 2–MATLAB-based software for topographic analysis and modeling in earth surface sciences. Earth Surf. Dyn. 2, 1–7. doi: 10.5194/esurf-2-1-2014

Jara-Muñoz, J., Melnick, D., Pedoja, K., and Strecker, M.R. (2019). TerraceM-2: A Matlab® Interface for Mapping and Modeling Marine and Lacustrine Terraces. Frontiers in earth Science 7(255). doi: 10.3389/feart.2019.00255.

Open PhD position in landscape evolution modelling

Posted on Updated on

Numerical models of landscape evolution (LEMs) are now increasingly sophisticated and able to model the interplay of various geomorphic processes and their controls. Parametrizing and testing these models using empirical data, however, remains challenging.

Markus Fuchs, professor for physical geography at the University of Giessen, Germany, now offers a PhD position. Equipped with a vast array of empirical data, the prospective candidate will aim at tackling some of the pressing methodological issues on model-data comparison that pertain validation and parametrization of LEMs. I have agreed to co-supervise the PhD student. If you are about to finish your Master degree and have skills in numerical computing and statistical analysis, we cordially invite you to apply for this position. If you are teacher or supervise students, then it’ll be great if you forward this open position.

The official job announcement in german and english language can be found here.