Landslides hit rivers

Posted on

River networks in steep terrain are frequently affected by landslides. Large landslides may deliver large amounts of debris to the river, dam the river and form lakes in the backwater. Landslide dams may become unstable, burst and generate devastating floods downstream. Some dams are remarkably stable and exist for several thousand years, inhibit upstream migration of knickpoints and protect upstream areas from fluvial incision and landscape denudation. Thus, landslide dams are an important geomorphological form in alpine environments.

The Bhote Kosi landslide is a very recent example of a valley blocking landslide. Please find insights and astonishing pictures in Dave Petley’s landslide blog and at ICIMOD.

In this blog post, I like to respond to a question that I was asked a couple of weeks ago. How is it possible to draw the locations of potential intersections of landslide runout zones with the river network in a plot of a river-length profile. I’ll use an hypothetical spatial distribution of landslides together with the Big Tujunga DEM to show how such a plot can be made using TopoToolbox.

First, let’s load the DEM, calculate flow direction (FLOWobj) and a STREAMobj. For simplicity, streams are simply derived by thresholding the flow accumulation grid. In addition, we derive a stream grid, e.g, a GRIDobj that contains a logical matrix. Nonzero entries in this matrix refer to stream pixels.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD  = FLOWobj(DEM,'preprocess','carve');
S   = STREAMobj(FD,flowacc(FD)>1000);
SG  = STREAMobj2GRIDobj(S);

Note that SG derived by STREAMobj2GRIDobj may not be equal to a grid calculated by the command SG = flowacc(FD)>1000. This is because streams must contain at least two pixels. Simply thresholding flow accumulation may sometimes result in single stream pixels along the DEM edges.

The landslide grid is hypothetical. We simply locate landslide source areas using 100 randomly distributed points. The function randperm is a perfect function to do this.

LS = GRIDobj(DEM);
LS.Z = logical(LS.Z);
LS.Z(randperm(numel(LS.Z),100)) = true;

Again, we calculate an influence map that, based on flow direction, routes the initiation locations downstream. Again, we transform this influence map to a STREAMobj, thus obtaining the potential sediment pathways starting at the landslide initiation zones.

LSi = influencemap(FD,LS);
Sls = STREAMobj(FD,LSi);

Now we need the locations were potential sediment pathways hit the river network. There is no function yet to accomplish this right away. Thus, we need to access the STREAMobj object properties. This provides an opportunity to shortly introduce these.

A STREAMobj contains various properties that contain the network topology and geometry of the river network. You can access the properties and their descriptions using the command

doc STREAMobj

This will open the help browser and will show, among other, the class properties and their descriptions.

IXgrid — [node attribute] linear index of stream nodes into instance of GRIDobj
cellsize — cellsize
distance — [node attribute] distance from outlet (dynamic property)

georef — additional information on spatial referencing
ix — [edge attribute] topologically sorted nodes (givers)
ixc — [edge attribute] topologically sorted nodes (receivers)
orderednanlist — nan-separated index into node attributes
refmat — 3-by-2 affine transformation matrix (see makerefmat)
size — size of instance of GRIDobj from which STREAMobj was derived
x — [node attribute] x-coordinate vector
y — [node attribute] y-coordinate vector

Similarly to FLOWobjs, STREAMobjs contain a downstream, topologically sorted edge list (ix and ixc). ix contains the givers and ixc the receivers. The IXgrid property is a node attribute list and provides access to the linear index into the DEM, from which the STREAMobj was derived.

We now want to identify the edges of the potential sediment network whose receiver nodes are located on the stream network and whose giver nodes are sited on hillslopes. This is accomplished by the two following lines.

IX  = SG.Z(Sls.IXgrid(Sls.ixc)) & ~SG.Z(Sls.IXgrid(Sls.ix));
IX  = Sls.IXgrid(Sls.ixc(IX));

In order to calculate coordinate pairs from the linear indices IX, use the function ind2coord.

[x,y] = ind2coord(DEM,IX);
hold on
plot(Sls,'Color',[.8 .8 .8])
hold off
legend('potential sediment pathways',...
    'stream network',...
    'potental landslide dam locations')


Finally, we can use the linear indices as parameter value in the function plotdz. This will finally plot the longitudinal stream profiles together with the potential dam locations.



While writing this post, I realized that this calculation is not as straightforward as I had thought before. I think it would be much better to put all this in an own function so that linear indices of potential dam locations are directly obtained by

IX = funname(S,Sls);

What could be a good function name to replace funname?


Export STREAMobj to table

Posted on

In a last post I described how to export a STREAMobj to a shapefile. In this post, I’ll show how to write a STREAMobj to a delimited ascii text file or excel file.

Again, let’s get the data first and derive a STREAMobj.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM,'preprocess','carve');
A = flowacc(FD);
S  = STREAMobj(FD,A>500);
S = klargestconncomps(S);

hold on
hold off
stream network and digital elevation model
stream network and digital elevation model

Now we need the function STREAMobj2XY. Admittedly, some of the functionality of this function does not open up by the function name itself and I am thinking of renaming the function or writing a new function. So what does this function?

[x,y] = STREAMobj2XY(S);

With one input argument, it returns two vectors with x- and y-coordinates that you can use to plot the planform pattern of the stream network.

stream network
stream network

The fact that several lines are drawn with only two vectors is accomplished by separating the vectors with nans.

Now, what happens if you want an additional vector that contains the z-coordinate and possibly the distance from the outlet. That’s easy, simply add the data to the function STREAMobj2XY.

[x,y,z,d] = STREAMobj2XY(S,DEM,S.distance);

So x and y are again exported. The z vector is taken from the DEM and the variable d is calculated from an node attribute list (nal) which we retrieve as object property .distance. Hence, you can supply STREAMobj2XY with GRIDobj (which must have the same geometric properties as the DEM from which the FLOWobj and STREAMobj have been derived) and node attribute lists. Node attribute lists are obtained from some STREAMobj methods such as STREAMobj/gradient or STREAMobj/streamorder.
Finally, if you like to export this data to a textfile or excel file, simply concatenate the data to a matrix and write it to the disc using xlswrite or dlmwrite. As an example, here I write the data to a text file with tab-separated columns.

data = [x y z d];

New paper out in The Cryosphere Discussions

Posted on Updated on

Meltwater transfered from the surface of glaciers to the bed determines to a large extent ice velocities. Understanding the interaction between supraglacial and subglacial hydrological systems is thus crucial for an understanding of the dynamics of glaciers and their respond to climate change. In a recent paper first-authored by Caroline Clason, we present the results of spatially-distributed modelling for prediction of moulins and lake drainages on the Leverett Glacier in south-west Greenland. Please read and discuss. It’s open for discussion.

BTW, find some excellent pictures of Greenland’s supra-glacial lakes and moulins here.

Clason, C.C., Mair, D., Nienow, P., Bartholomew, I., Sole, A., Palmer, S., Schwanghart, W. (2014): Modelling the transfer of supraglacial meltwater to the bed of Leverett Glacier, southwest Greenland. The Cryosphere Discussions, 8, 4243-4280. [DOI: 10.5194/tcd-8-4243-2014]

New paper out in Icelandic Agricultural Sciences

Posted on Updated on

My interest in statistics, data analysis and numerical modelling has frequently spurred collaborations outside the field of geomorphology. In a recent endeavor of colleagues from University of Basel, the Agricultural University of Iceland and the Soil Conservation Service of Iceland, we used field work and data analysis to investigate below-ground biomass distribution of mountain birch in southern Iceland. Please find the paper published in Icelandic Agricultural Sciences here. It’s open access!

Hunziker, M., Sigurdsson, B.D., Halldorsson, G., Schwanghart, W., Kuhn, N. (2014): Biomass allometries and coarse root biomass distribution of mountain birch in southern Iceland. Icelandic Agricultural Sciences, 27, 111-125.

Geosimulator at the University of Potsdam

Video Posted on

Here is a short video documentation of a teaching project at the University of Potsdam where I have been involved. Have fun watching (it’s in german)

Export STREAMobj to shapefile

Posted on Updated on

I was recently asked, how to export a STREAMobj to the disk so that it is possible to edit and analyse it using different software. Clearly, the answer depends on what kind of software you want to use. Here, I describe how to export a STREAMobj to shapefile for working in a GIS. In the next post, I will show how to export STREAMobj to a nan-separated table.

Ok, let’s get the data first and derive a STREAMobj.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM,'preprocess','carve');
A = flowacc(FD);
S  = STREAMobj(FD,A>500);
S = klargestconncomps(S);

hold on
hold off
stream network and digital elevation model
stream network and digital elevation model

Shapefiles are a popular file format for vector data. If you want to work with your stream network in any geographic information system (GIS), export your STREAMobj to a shapefile. Using MATLAB and TopoToolbox, this requires two commands. First,

MS = STREAMobj2mapstruct(S);

I know the function name may not be intuitive. However, for those of you who have been working with the Mapping Toolbox, mapstructs may be known. Mapstructs are structure arrays that contain point, line and polygon vector data. Each element in mapstructs contains one object which includes X and Y data, and attributes.
You can plot the mapstruct easily with following command (Mapping Toolbox required)

stream network
stream network

Second, use the MATLAB Mapping Toolbox function shapewrite to export the mapstruct to a shapefile. Note that if you don’t have the Mapping Toolbox there may be functions available on the File Exchange that have same or similar functionality.


You will now be able to open streamnetwork.shp in QGIS, ArcGIS, SAGA GIS, or wherever you like.
There are, however, a few drawbacks with this approach. First, we didn’t include attributes; second, we didn’t take control on how our objects should look like. Is there only one multipart object? Are there as many objects as individual streams? Neither nor.


MS =

1×248 struct array with fields:

At the moment, MS has 248 objects. In addition, there are three additional attributes named streamorder, IX and tribtoIX. Hence, what STREAMobj2mapstruct does is to cut the stream network into objects based on their (Strahler) stream order. Where a reach changes its stream order it will be placed in a new object. IX contains for each object a unique identifier and tribtoIX indicates the identifier of the downstream object. If there is no downstream object, tribtoIX will be zero.
Yet, there are cases when you wish to export shorter, equally long stream segments and add attribute data other than stream order. This is possible. Just take a look at the help of STREAMobj2mapstruct or read on. I’ll show how to export a STREAMobj to a shapefile including the attributes channel gradient, drainage area and stream order for 500 m long segments.
So let’s first calculate gradient and streamorder.

G = gradient8(DEM);
so = streamorder(S);

Then, again call STREAMobj2mapstruct. This time use additional parameter name/value pairs.

MS = STREAMobj2mapstruct(S,'seglength',500,'attributes',...
        {'uparea' A @median ...
         'gradient' G @mean ...
         'streamord' so @max});

The first parameter name/value (pvpn) pair should be self-explanatory. It tells the function the length of each line object in the resulting mapstruct. The second pvpn pair is more complicated. It defines the attributes and requires a cell array as value. You can create the cell using curly brackets. The elements of the cell array must have following structure

{attname1 variable1 aggfunction1 attname2 variable2 aggfunction2 …}

where attname1 is the name of the first attribute as it will appear as field name in the attribute table. Note that field names should not be longer than 10 letters. variable1 is a workspace variable and must be either a GRIDobj ( this GRIDobj must have the same geometry as the DEM from which FLOWobj and STREAMobj was derived) or a node attribute list of a STREAMobj (I’ll explain node attribute lists in one of the next posts. They will be a new but old feature in the next release). aggfunction1 is a function that is used to aggregate data available for each vertex in the STREAMobj into the 500 m long segments. This must be a function that takes a vector and returns a scalar such as @mean, @max, or @(x) prctile(x,5). These are anonymous functions. If you are new to anonymous functions, check the MATLAB help.


MS =

605×1 struct array with fields:

You can see that MS now contains 605 objects and includes the desired attributes.
Finally, let’s check whether the object length is always 500 m.

seglength = zeros(numel(MS),1);
for r = 1:numel(MS);
    seglength(r) = nansum(hypot(diff(MS(r).X),diff(MS(r).Y)));

xlabel('length [m]')
Histogram of object lengths
Histogram of object lengths

The segment length is not always equal to 500 m. However, most objects in the mapstruct have length around 500 m. This is because STREAMobj2mapstruct attempts to distribute separators in a way that warrants homogeneous lengths while segments are additionally separated at stream confluences.
Above I have used for-loops to calculate the object lengths (seglength). There is an easier way to do so by calculating stream length as attribute in the function call to STREAMobj2mapstruct. Do you have an idea how to do this? Leave a comment!

Video on our research in Nepal

Video Posted on

I am currently involved in a project that studies the Pokhara Formation in Nepal. Here is a short movie made by Melanie Langpap during our field stay. Enjoy watching!