geoglobe now in MATLAB R2020a

Posted on

In a previous post (here) I have shown how TopoToolbox is able to export data to kml-files which can be directly opened in Google Earth (or other digital globes such as ArcGIS Earth). Indeed, kml is a great way to share geographic data. However, if you want to visually explore your data, the step via kml and Google Earth may not be necessary any longer. Since its latest release, MATLAB’s Mapping Toolbox includes geoglobe, a geographic globe that allows navigating on the Earth surface and to add data. To this end, geoglobe will enable numerous applications that enhance the way we can explore data analyzed and generated using TopoToolbox.

geoglobe opens in a figure created with uifigure.

h = uifigure;
g = geoglobe(h);


Once initiated, you can start with the GeographicGlobe object g. For example, you can change the basemap, which is by default ‘satellite’.

g.Basemap = 'topographic';


There are numerous ESRI basemaps available, as well as some from Natural Earth. Terrain is based on the GMTED2010 model by the U.S. Geological Survey (USGS) and National Geospatial-Intelligence Agency (NGA) and hosted by MathWorks. The DEM has a spatial resolution of 250 m, which is sufficient for many large-scale applications. If you want to add higher resolution data, you can add DTED files with terrain data.

Finally, you can add your own geographic data. Here I show, how you can add stream network data generated with TopoToolbox. I will download some example data from Taiwan using the function readexample:

DEM = readexample('taiwan');
DEM = inpaintnans(DEM);
FD  = FLOWobj(DEM);
S   = STREAMobj(FD,'minarea',5000);
[lat,lon,z] = STREAMobj2latlon(S,DEM);


There are still some slight problems in visualizing the network because the lines remain partly hidden beneath the terrain. I reduced but not fully resolved this effect by adding 10 m to the elevations of the stream network. Notwithstanding, geoglobe is a great utility to visualize topographic data and I am sure there are more functions to come in the future that enable us to interactively work with geographic globes. So far, I particularly like how The Mathworks implemented navigation on the globe which works super smoothly.

Landscapes Live

Posted on

#shareEGU20 is currently setting an example. Virtual conferences are an effective means to bring scientists together and to spur vivid discussions.

To keep this scientific exchange alive during COVID-19, Philippe Steer, Vivi Pedersen, Stefanie Tofelde, Pierre Valla, Charlie Shobe, and me (my role was actually very minor) have initiated Landscape Live, a new remote seminar series focused on sharing exciting geomorphology research throughout the international scientific community.

The remote format allows for free, planet-conscious, and pandemic-proof attendance by anyone who is interested. Talks will take place using the Zoom meeting software.

Talks will happen during the academic year in blocks of 5-6 weekly talks, with long breaks in between. The first block will run from 28th May through 25th June 2020, with the next block taking place during the fall 2020 semester.

Each seminar in the first block will be held on a Thursday at 2 pm GMT (4 pm Central European Time). The first block features a fantastic lineup of speakers as follows:

  • Thursday 28th May at 2 pm GMT: Georgina Bennett, University of Exeter
  • Thursday 4th June at 2 pm GMT: Anneleen Geurts, University of Bergen
  • Thursday 11th June at 2 pm GMT: Liran Goren, Ben Gurion University of the Negev
  • Thursday 18th June at 2 pm GMT: Robert Hilton, Durham University
  • Thursday 25th June at 2 pm GMT: Fiona Clubb, Durham University

Please visit for the most up-to-date information, including the links to each Zoom meeting.

Suggestions for future speakers are welcome; please feel free to send names to any member of the organizing committee. We look forward to seeing you (virtually) at Landscapes Live!

This text is slightly modified from the version written by Charlie Shobe and previously published on the EGU geomorphology blog.


Posted on


Only two days left before #shareEGU20 opens its digital gates. I surely won’t be able to be online during the whole event, but I’ve my personal schedule of displays which I’ll try to view and discuss online inbetween homeoffice, homeparenting, homeschooling, …

If you haven’t a fixed schedule yet, consider getting involved in following displays which I authored or coauthored:

EGU2020-11177 – Divide mobility controls knickpoint migration on the Roan Plateau (Wolfgang Schwanghart and Dirk Scherler), Wed, 06 May, 14:00–15:45 | D1315, DISPLAY

EGU2020-19260 – The TopoToolbox v2.4: new tools for topographic analysis and modelling (Dirk Scherler and Wolfgang Schwanghart), Thu, 07 May, 10:45–12:30 | D818, DISPLAY

EGU2020-5609 – Uncertainties in Chi analysis: implications for drainage network and divide stability (Jens Turowski et al.), Tue, 05 May, 10:45–12:30 | D1132, DISPLAY

EGU2020-5900 – Evaluating the effect of variable lithologies on rates of knickpoint migration in the Wutach catchment, southern Germany (Andreas Ludwig et al.), Tue, 05 May, 08:30–10:15 | D1114, DISPLAY

EGU2020-8811 – Why do shelf-incising submarine canyons form? – Insights from global topographic analyses and regression trees (Anne Bernhardt and Wolfgang Schwanghart), Fri, 08 May, 08:30–10:15 | D1102, DISPLAY

EGU2020-3737 – Illuminating the speed of sand – quantifying sediment transport using optically stimulated luminescence (Jürgen Mey et al.), Fri, 08 May, 08:30–10:15 | D899, DISPLAY

New paper out: Divide mobility controls knickpoint migration

Posted on Updated on

In recent years, there has been a quite fierce debate about how landscapes evolve in response to lateral dynamics of river networks. These dynamics include laterally shifting rivers, their expansion or contraction in upstream and downstream direction, and the mobility of catchment divides. In fact, as we seek to gain insight into changes in climate and tectonics from the analysis of river networks, we often make the assumption that theses river networks are static, that their spatial configuration remains stable.

In our paper published now in Geology (Schwanghart and Scherler, 2020), we argue that this assumption must be made cautiously. We revisited the Parachute Creek basin, Colorado, US, which has been studied by Berlin and Anderson (2007). The site is a real natural laboratory, because it has a uniform climate as well as (sub)horizontally uniform bedrock which makes it possible to analyze the phenomenon of knickpoint migration into a relict landscape in isolation of other factors that govern knickpoint retreat.

View towards the west on the Parachute Creek basin. Imagery from Google Earth.

In this study, we applied the stream-power incision model to infer present day locations of knickpoints in river profiles of the Parachute Creek basin (see Figure above). The knickpoints emanated from a base level drop at the Upper Colorado River about ~8 Mio years ago. What we realized when looking at the misfit of predicted locations and actual locations, was that there is one subbasin where knickpoints consistently travelled further than our model would predict. Now, you would expect some randomness in knickpoint locations, but the systematic spatial pattern grasped our attention. Looking more closely at the DEM revealed that this subbasin shares some lengths of its divide with the plateau margin, a steep cliff that goes down to the Colorado River. That this margin is actively retreating became clear from numerous beheaded valleys visible in the DEM and Google Earth imagery. A hypothesis was quickly formulated: Knickpoints in this subbasin had migrated further than we expect, because part of the subbasin’s area had been lost to cliff retreat.

In the next few blog posts, I will talk a bit more about our approach to calculate the area lost, and how we actually could even gain some constraints on the timing of area loss. Until then, you may also check Dirk’s and my papers on divide mobility just published in ESURF (Scherler and Schwanghart, 2020a,b). We’ll write here about this two-part paper in due time.


Berlin, M. M. and Anderson, R. S.: Modeling of knickpoint retreat on the Roan Plateau, western Colorado, Journal of Geophysical Research: Earth Surface, 112(F3), doi:10.1029/2006JF000553, 2007.

Scherler, D. and Schwanghart, W.: Drainage divide networks – Part 1: Identification and ordering in digital elevation models, Earth Surface Dynamics, 8(2), 245–259, doi:, 2020a.

Scherler, D. and Schwanghart, W.: Drainage divide networks – Part 2: Response to perturbations, Earth Surface Dynamics, 8(2), 261–274, doi:, 2020.

Schwanghart, W. and Scherler, D.: Divide mobility controls knickpoint migration on the Roan Plateau (Colorado, USA), Geology, doi:10.1130/G47054.1, 2020. Supplementary material can be found here.

Calculating Hack’s Law using TopoToolbox

Posted on

Hack’s Law describes an empirical relationship between river length and drainage area (Hack 1957). The functional relationship is a power function with the equation L = c A^h where L is the length of the longest stream from the outlet to the divide, A is the drainage area above a particular locality, c is a constant, and h is the scaling exponent (see figure).

Measuring stream length and drainage area (Hack 1957).

In general, h is slightly below 0.6 (Rigon et al. 1996). Today, I show how to calculate the parameters of Hack’s Law. First, I will demonstrate how the parameters of Hack’s Law are derived for a single catchment. In the second part, I will then apply the technique to calculate the parameters for catchments draining Taiwan.

Big Tujunga catchment

In the first example, I will use the DEM of the Big Tujunga catchment, which is part of the TopoToolbox distribution. While Hack’s Law appears to hold for any point inside a basin (Rigon 1996), I am showing here how to derive it from set of pixels that comprise the longest path from the divide to the catchment outlet. As usual, I derive flow directions (FLOWobj) and the stream networks (STREAMobj). Then, I extract the largest basin and subsequently the longest stream in this basin that extends from the divide to the outlet.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD  = FLOWobj(DEM);
DIST   = flowdistance(FD);
S   = STREAMobj(FD,'minarea',1000);
S   = klargestconncomps(S);
D   = drainagebasins(FD,S) > 0;
DIST = clip(DIST,D);
[~,IX] = max(DIST.Z(:));
S   = STREAMobj(FD,'channelheads',IX);
imageschs(DEM,D,'falsecolor',[.0 0 0],'truecolor',[1 1 1], ...
hold on
hold off


In a second step, I calculate drainage area and distance. The function flowacc computes the flow accumulation (in pixels) for the entire grid, and the function getnal extracts the grid values for each node in the river network. In addition, I calculate flow distance using the function distance. By default, the function calculates the distance from each node to the outlet, but there are also other options. Here, we use the maximum distance from the channelhead.

A   = flowacc(FD)*DEM.cellsize^2;
a   = getnal(S,A);
d   = distance(S,'max_from_ch');

Now we can plot the two node-attribute lists versus each other.

xlabel('Area [m^2]')
ylabel('Stream length [m]')


Now lets use these data to fit Hacks Law. I am using nlinfit from the Statistics and Machine Learning Toolbox for this purpose which requires initial estimates of the parameters. The choice of these initial values can be challenging and way-off values may render nlinfit unable to converge to a good solution.

b = nlinfit(a,d,@(b,X) b(1)*X.^b(2),[0.002 0.6]);

b is a two-element vector and contains the parameters of Hack’s Law, where b(1) is the constant and b(2) is the scaling exponent. Here, b(1) is 0.0035 and b(2) is 0.8346. Having these parameters, we can now obtain estimates for distances for any value of drainage area.

hold on
fplot(@(X) b(1)*X.^b(2),xlim)
hold off
xlabel('Area [m^2]')
ylabel('Stream length [m]')


You may note that nlinfit returns an exponent different from Hack’s value of 0.6 and there may be different reasons for this. First of all, the basin may have a particular form because it extends on the relict landscapes of the Chilao Flats and features a series of knickpoints that suggest that it experiences a transient response to accelerated uplift. Second, the river contains both mixed bedrock-alluvial sections in the upper part, and a purely alluvial lower part. Finally, our data violates some of the assumptions that underly the inference of parameters, i.e. a identically and independent distribution (i.i.d.). There may be ways to at least partly account for these violations by random sampling of our data, but I don’t want to go into too much detail here.


In the second example, we will use the SRTM-3 DEM of Taiwan, which you can download using the readexample function. In this example, we calculate the parameters in Hack’s Law from multiple drainage basins.

DEM = readexample('taiwan');

We first inpaint missing values in the DEM and use FLOWobj to derive flow directions. In a second second step, we then calculate the trunk stream for each basin, again up to the divide. This is a bit tricky, if you have multiple drainage basins. Here, I show how to do this using FLOWobj2cell and cellfun. I omit many small drainage basins that occur along the coast of Taiwan by setting a minimum threshold of drainage basin are to an arbitrary size of 1000 pixels.

DEM = inpaintnans(DEM);
FD  = FLOWobj(DEM);
[CFD,~,a] = FLOWobj2cell(FD);
CFD = CFD(a>1000);
[~,IX]  = cellfun(@(fd) max(flowdistance(fd)),CFD);

S   = STREAMobj(FD,'channelheads',IX);
imageschs(DEM,[],'colormap',[1 1 1],'colorbar',false,'ticklabels','nice');
hold on


Again, I calculate upslope area and distance and plot both in semilogarithmic axes.

A  = flowacc(FD)*DEM.cellsize^2;
a  = getnal(S,A);
d  = distance(S,'max_from_ch');
OUTLETS = streampoi(S,'outlets','logical');
a  = a(OUTLETS);
d  = d(OUTLETS);
xlabel('Area [m^2]')
ylabel('Stream length [m]')
hold off


Then, I again use nlinfit to determine the parameters of Hack’s Law.

b = nlinfit(a,d,@(b,X) b(1)*X.^b(2),[0.002 0.6]);

hold on
fplot(@(X) b(1)*X.^b(2),xlim)
hold off
xlabel('Area [m^2]')
ylabel('Stream length [m]')

This time, nlinfit yields 0.5315 as exponent.


Finally, we might want to look at the spatial patterns of residuals from our predicted lengths.

rerr = (d - (b(1)*a.^b(2)))./d;
IXoutlet   = S.IXgrid(OUTLETS);
D = drainagebasins(FD,IXoutlet);
D.Z = double(D.Z);
D.Z(D.Z ~= 0) = rerr(D.Z(D.Z~=0));
D.Z(isnan(DEM.Z)) = nan;
imageschs(DEM,D,'caxis',[-1 1],'colorbarylabel','Rel. error')


The map shows that the residuals (or relative errors) from Hack’s Law. We can see that some catchments, particularly those draining the western slopes tend to have longer distances of the main river than predicted by Hack’s Law. These catchments are often rather elongated which might partly due to their coverage of large alluvial areas and lowlands.


Hack, J. T.: Studies of longitudinal stream profiles in Virginia and Maryland, USGS Professional Paper, 295, 45–97, 1957.

Rigon, R., Rodriguez-Iturbe, I., Maritan, A., Giacometti, A., Tarboton, D. G. and Rinaldo, A.: On Hack’s law, Water Resources Research, 32, 3367–3374, 1996.


Open PhD position: ‘Geochronology of the first Eurasian Ice Sheets’

Posted on

Hardanger Vidda (Photo credit: John Jansen).

Want to work with a great advisory team on the geochronology of the first Eurasian ice sheets? Well, here is your chance. My friend and colleague John Jansen (Czech Academy of Sciences) and his colleagues from Charles University and Aarhus University offer a highly exiting PhD opportunity. You will be working with geochronological techniques (cosmogenic nuclide dating) and numerical modelling to unravel the elusive traces of early Eurasian glaciations. There is definitely a lot to discover here.

Check the job advertisement here and apply!

Handling closed basins

Posted on Updated on

Much of the world’s terrestrial area is comprised of endorheic basins. These basins have no drainage to the oceans. Instead, water in these basins seeps into the ground or evaporates (wikipedia).

Internally drained basins come in different sizes. The Caspian Sea, the largest inland body of the world, has a catchment area of 3,626,000 sqkm. In contrast, dolines – small hollows in karstic terrain – have just a few tens of square meters.

In digital elevation model analysis, internally drained basins can be quite challenging, largely because we often do not know whether a closed basin in our DEM is a true closed basin or an artefact of the data.

When calculating flow directions in TopoToolbox, the default setting is that all closed basins are artefacts. However, when working in dryland areas, we must be cautious with this assumption.  In this post, I will show how to deal with closed basins. In addition, I’ll show how to use MATLAB Live Scripts that help us to interactively explore parameters that control the identification of closed basins.

To get started, please download the Tibetan Plateau example using the readexample function:

DEM = readexample('tibet');

The DEM features numerous closed basins that we identify using the fillsinks function. We then visualize the difference between the filled DEM and the actual DEM, which gives us an idea about the location and depth of internally drained basins.

DEMf = fillsinks(DEM);
'colormap',flipud(hot),'colorbarylabel','Sink depth [m]')


In a second step, we will investigate the properties of the closed basins in more detail. Here we look at two properties. The area of the sinks, and their depth. The function regionprops which is part of the Image Processing Toolbox comes in handy here.

SINK  = DIFF > 0;
stats = regionprops(SINK.Z,DIFF.Z,'Area','PixelIdxList','MaxIntensity');
sinkarea = [stats.Area]*DEM.cellsize^2;
sinkdepth = [stats.MaxIntensity];
xlabel('Area of sink [m^2]')
ylabel('Maximum depth of sink [m]')


There are more than 11,000 sinks in the DEM. Clearly, not all of them are true sinks. I propose that true sinks should have a large areal extent and be rather deep. Hence, we can apply some thresholds to classify sinks into true sinks and artificial sinks. As a first guess, lets assume a combination of an areal extent of 10 km2 and a depth of 10 m as suitable parameters of our sink identification model.

areathres = 10^7;
depththres = 10;
sinkstats = stats(sinkarea > areathres & sinkdepth > depththres);
sinkpixels = {sinkstats.PixelIdxList};

These thresholds returns six closed basins. In a third step, we calculate flow directions in a way that acknowledges the internal drainage of these basins. There are different ways to do this. Here, we use the trick that sets the pixels with the minimum elevation in each sink to nan. There are other ways to do this (see the function FLOWobj), but this approach is actually the fastest.

% A: Here we identify the linear indices of pixels 
% with minimum elevation for each sink
minelevix  = cellfun(@(ix) ix(find(DEM.Z(ix) == ...
% B: Make a copy of the DEM
% C: Set elevations of minimas to nan
DEM2.Z(minelevix) = nan;
% D: calculate flow directions
FD = FLOWobj(DEM2);
% E: calculate flow accumulations
A  = flowacc(FD);
% F: and visualize the results.
hold on
[x,y] = ind2coord(DEM,minelevix);
hold off


Finally, we may not be satisfied with the results. Based on our parameters, there may be to many or to few closed basins. Below is a gif movie that shows how to use the interactive tools in MATLAB Live Scripts to vary the sink-area and sink-depth parameters and instantly see their effects on the resulting flow network.


I hope you enjoyed this post. Applying TopoToolbox in regions with internally drained basins should now be no problem anymore.