FLOWobj now supports divergent flows

Posted on Updated on

Divergent flows on hillslopes.

I had this request numerous times, and now it’s finally there. FLOWobj, the object that stores flow directions in TopoToolbox now handles multiple flow directions (MFD). MFD entails that flows are distributed from each cell to all or selected downward neighbors. Thus, flows can diverge and terrain metrics such as flow accumulation better represent the actual flow patterns observed on hillslopes, fans, or debris cones. A previous work-around is now obsolete.

FLOWobj with multiple flow directions can be derived using two different algorithms:

  • ‘multi’: Flow is partitioned to all downward neighbors according to slope.
  • ‘Dinf’: The D infinity approach by Tarboton (1997) that distributes flow to the steepest triangular facet. I have adopted the code by Steve Eddins.

The syntax is easy:

FD = FLOWobj(DEM,'multi');


FD = FLOWobj(DEM,'Dinf');

Note, however, that none of the preprocessing options (fill or carve) are possible. You thus must preprocess the DEM before! One way to do this is shown in the second example in the help of FLOWobj. Also note that not all FLOWobj methods can deal with MFD. There may still be bugs although I have updated several functions to either throw error messages or be able to handle MFD. Let me know if you encounter any issues.


Tarboton, D. G. (1997). A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research, 33(2), 309-319.

Eddins, S. (2016). Upslope area function. Mathworks File Exchange. [link]

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.

Changes to imageschs

Posted on

imageschs is one of the functions that I most frequently use. Admittedly, it has a somewhat non-intuitive name. Yet, the function’s name is simply the concatenation of imagesc, a MATLAB built-in function to display an image with scaled colors, and hs, which abbreviates hillshading. In combination, imageschs displays an image (or GRIDobj) with a shading based on an underlying digital elevation model (DEM).

The syntax of imageschs is quite simple, but there are a number of parameter name/value pairs that you may use to control some additional behavior of the function. The parameter name/value pairs are documented in the help section of the function that can be accessed via

help imageschs

Here I briefly outline some new functionality that I added lately.

1. Median filtering

Smoothing the hillshade has proven successful to reduce the effect of small geographic features on shading, instead emphasizing the larger scale geographic patterns. Many software packages use a median filter and imageschs now does so as well:

ax(1) = subplot(1,2,1);

ax(2) = subplot(1,2,2);

2. Percent clipping

By default, imageschs linearly scales the color range between the minimum and maximum value of the coloring image or GRIDobj. Sometimes, this range may be strongly affected by just a few extreme values. This is particularly true for curvature, that usually contains values close to zero, but often has few values which are very high or low. The linear scaling inbetween these values entails that all patterns characterized by smaller variations become invisible. imageschs can now remove those values using percentage clipping. That means that the scaling applies only to the value range between the X%’s and 100-X%’s percentile of the data.

C = curvature(DEM);

3. Make the hillshade permanent

Have you ever been annoyed that imageschs takes a while if you have to call it repeatedly although you always use the same DEM? Now you can use the option ‘usepermanent’ to reduce waiting time. If this option is set to true, imageschs will store the hillshade as a persistent variable and will make it available if the function is called next time. Of course, this will only work out if the DEM (actually, its hillshade) and the coloring matrix or GRIDobj have the same dimensions. And of coarse, beware that any changes that you make to the DEM in the meantime remain disregarded by imageschs.

DEM = resample(DEM,5);
DEM = crop(DEM);

G = filter(gradient8(DEM),'mean',11);
tic;imageschs(DEM,[],'usepermanent',true); toc
Elapsed time is 3.599666 seconds.
tic;imageschs(DEM,G,'usepermanent',true); toc
Elapsed time is 0.717102 seconds.

Ok, that’s it about recent changes to imageschs. Please leave a comment or report bugs!

First Open Call of Geo.X – PhD and PostDoc positions in “Geo Data Science”

Posted on

Geo.X network

Geo.X, a Berlin-Potsdam based research cluster, advertises 14 positions in a new research school. The focus of this research school is on the interdisciplinary area of Data Science as a methodological research area within the geoscientific research topics Natural Hazards and Risks and Geo-Bio-Interactions. As readers of this blog, I am sure, these positions will certainly be of interest to many of you.

Please find further information here and apply here.

Let’s go dynamic with TTLEM!

Posted on Updated on

As geomorphologists we have to live with the fact that we can rarely observe the processes that shape the Earth’s surface. Instrumental records cover only a minor time span over which many geomorphic processes act, thus challenging our abilities to disentangle and understand the complex interactions of landscape evolution. Instead, much of our work has to rely on the analysis of landforms, their geometry, their assemblage, and their constituting material, as well as a blend of geochemical and numerical dating techniques of ever-increasing sophistication.

Numerical landscape evolution models (LEM) provide a useful approach to this challenge. LEMs are simulation tools that attempt to model erosion, sediment transport and deposition as well as feedbacks with vegetation and land use. They amalgamate our state of knowledge in a set of physically-based mathematical formulations and allow us to test whether these equations and their parameter values are able to generate output that is plausible and consistent with field evidence.

TTLEM, the TopoToolbox Landscape Evolution Model, is the latest addition to TopoToolbox. Spearheaded by Benjamin Campforts from KU Leuven, we have developed a LEM that allows us to simulate how mountains grow, rivers incise and hillslopes respond to tectonic and climatic forcing. We placed particular emphasis on implementing numerical models that minimize numerical diffusion. To achieve this, Benjamin has adopted a higher order flux limiting total volume method that is total variation diminishing (TVD-TVM) (Campforts and Govers 2015) to solve the partial differential equations of river incision and tectonic displacement.

In our recent manuscript under discussion in ESurf, we show that using these methods is more than just an exercise in numerical modelling. First-order approximations often smooth knickpoints in uncontrollable ways that impact on derived catchment wide erosion rates. Numerical diffusion strongly affects lateral tectonic displacement, thus restricting its simulation to models that use irregular grids and lagrangian approaches. The TVD-TVM approach solves for these issues by reducing numerical diffusion to a minimum, and thus offers a regular-grid based model with wide applications in tectonic geomorphology.

Needless to say that you can directly analyze and visualize the output using TopoToolbox. The implementation comes with a set of examples that you can directly run from the command line. Give it a try and let us know what you think!


Campforts, B., Govers, G., 2015. Keeping the edge: A numerical method that avoids knickpoint smearing when solving the stream power law. Journal of Geophysical Research: Earth Surface 120, 1189–1205. doi:10.1002/2014JF003376

Campforts, B., Schwanghart, W., Govers, G. (2016): Accurate simulation of transient landscape evolution by eliminating numerical diffusion: the TTLEM 1.0 model. Earth Surface Dynamics Discussion, in review. [DOI: 10.5194/esurf-2016-39]

Himalayan hydropower faces higher uncertainty closer to glacial lakes

Posted on

Glacial lakes (blue circles) and hydropower projects (yellow = existing, red = planned/under construction) along the Himalayan Arc (modified from Schwanghart et al., 2016).

The Himalayan nations’ thirst for energy is increasing. Hydropower is among the most important energy sectors and has experienced massive growth rates during the last years. This growth is expected to further increase in the coming years. More than 400 Himalayan hydropower projects (HPP) are currently registered or at validation for the Clean Development Mechanism (CDM) alone, a global environmental investment and credit scheme to enable countries with emission-reduction commitments to implement emission-reduction projects in developing countries.

In our now published paper in Environmental Research Letters, we show that – as opportune sites along major rivers are already occupied – numerous planned or currently constructed HPP expand into river reaches higher upstream and closer to glacial lakes. These lakes, dammed behind moraines, may sporadically burst upon dam failure and create so-called glacial lake outburst floods (GLOFs). We mapped more than 2000 lakes dammed by terminal or lateral moraines along the Himalayan Arc and calculated outburst scenarios for each using a physically based dam breach model. Since many of the variables of this model are impossible to constrain from remote sensing data alone, we adopted a probabilistic approach and used a Monte-Carlo simulation to derive entire distributions of outburst scenarios. Using TopoToolbox, we determined the flow paths of potential GLOFs from each lake and used these paths as input to a 1D flood propagation model that we fed with the results of the Monte-Carlo simulation. This allowed us to track how uncertainties about lake bathymetry and dam characteristics propagate downstream, and showed that these uncertainties remain significant, for most lakes, until a distance of 80 km downstream. This is the distance within which more and more HPPs are planned or currently constructed.

Clearly, our approach is based on a series of assumptions and limitations, and I encourage you to read the full paper (open access) for a detailed discussion of these. We conclude that the current massive expansion of HPP closer to glaciated areas entails that our lack of knowledge about the stability and geometry of glacial lakes will make it increasingly difficult to estimate potential flood magnitudes and thus compromise reliable risk assessment and safety planning of new HPP projects.

Schwanghart, W., Worni, R., Huggel, C., Stoffel, M., Korup, O. (2016): Uncertainty in the Himalayan energy-water nexus: estimating regional exposure to glacial lake outburst floods. Environmental Research Letters, 11, 074005. [DOI: 10.1088/1748-9326/11/7/074005]

Local-scale flood mapping using LiDAR

Posted on

Flood maps are increasingly important for managing flood risks and land use in riverine areas. Much of the efforts on flood mapping have concentrated on large and extensive floods. Localised flooding, however, has received less attention, although those local-scale inundation patterns are often prolonged and frequent in low-land areas and can incur high economic losses due to crop, property and infrastructure damage.

In a recent paper first-authored by Radoslaw Malinowski, we have used full-waveform airborne laser scanning to map inundation patterns of the Norrea River, a low-land river in Denmark. Our particularly focus was on understanding how water effects the backscattered energy at variable water coverages in areas with interspersed vegetation. We show that the backscattered energy negatively correlates with water coverage. We exploited that relation together with machine learning techniques combining radiometric and geometric data to obtain a LiDAR-pointcloud based flood map with a very high mapping accuracy of ~89%. Our results demonstrate the potential of using full-waveform LiDAR data to detect water surfaces on floodplain areas with limited water surface exposition due to vegetation canopy.

Malinowski, R., Höfle, B., König, K., Groom, G., Schwanghart, W., Heckrath, G. (2016): Local-scale flood mapping on vegetated floodplains from radiometrically calibrated airborne LiDAR data. ISPRS Journal of Photogrammetry and Remote Sensing, 119, 267-279 [DOI: 10.1016/j.isprsjprs.2016.06.009].