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!

Multiple flow directions in TopoToolbox 2

Posted on Updated on

With the release of version 2 of TopoToolbox, multiple flow directions (MFD) are not available anymore, and I got frequent requests whether MFD will be included in the future. Of course, MFD will be included in one of the next releases. Until then, here is a quick work-around. Note, however, that FLOWobjs with MFD are not supported by all functions (e.g., drainagebasins) and may show undesired behavior when used with some functions.

You’ll need TopoToolbox 1.06 or earlier versions which you can download here. Then use following code:

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
[dem,X,Y] = GRIDobj2mat(DEM);
[X,Y] = meshgrid(X,Y);
M = flowdir(X,Y,fillsinks(double(dem)));
FD = FLOWobj(M,'size',DEM.size,'refmat',DEM.refmat);
A = flowacc(FD);
Multiple flow direction
Multiple flow direction

For comparison, the single flow direction algorithm returns following flow patterns

FD = FLOWobj(DEM,'preprocess','carve');
As = flowacc(FD);
Single flow direction
Single flow direction

The SFD algorithm returns only rivers that are one pixel wide. The MFD algorithm may produce section with dispersive flow which might be more desirable on hillslopes, alluvial river sections and fans.

Paper spotlight: IDA: An implicit, parallelizable method for calculating drainage area

Posted on Updated on

In our 2010 paper on TopoToolbox published in Environmental Modelling and Software, we have introduced a matrix algebraic approach to calculate flow accumulation. In a recent paper in Water Resources Research, Alan Richardson and colleagues take this approach a step further by presenting a method that enables implicit, parallizable solving of the flow accumulation equation by a preconditioned iterative solver. Their approach leads to significant speed-ups particularly for processing large datasets on many-core computational systems such as large GPU clusters. In addition, this method is particularly suited in situations when drainage areas need to be calculated very frequently such as in landscape evolution models (LEMs) since upslope areas calculated in a previous time step can be used as initial guess for the next iteration. When drainage patterns stabilize during a simulation, this may speed up such calculations tremendously.

Modify your stream network

Posted on Updated on

One of the functions that I frequently use but which probably only few people know is the function STREAMobj/modify. As the name implies, the function modify lets you manipulate stream network data stored as a STREAMobj. The function allows you to make a number of choices among automated and manual editing procedures. Here is an overview:

Given one or more STREAMobj, you can
– extract streams with a given streamorder. The syntax for extracting all streams with Strahler streamorder less than 5 would be

Sm = modify(S,'streamorder','<5');

– extract a part of the network defined by the distance from the outlet. An example to extract the streams located between 50 and 100 km away from the outlet would be

Sm = modify(S,'distance',[50000 100000]);

– extract a part of the network that is upstream or downstream to regions specified in a GRIDobj. E.g., the syntax for extracting the stream networks above 1500 m is

Sm = modify(S,'upstreamto',DEM>1500);

– extract a part of the network that is tributary to another network. E.g., to extract the streams that are tributary to the trunk stream simply call

Sm = modify(S,'tributaryto',trunk(S));

Finally, there is the parameter ‘interactive’ which may take two options: ‘reachselect’ or ‘polyselect’. Check the two gif movies below how these work. Make following call if you want to manually select a reach

Sm = modify(S,'interactive','reachselect');
STREAMobj/modify with the option ‘interactive’ – ‘reachselect’


Or call the function modify in this way to select a stream network based on a polygon outline.

Sm = modify(S,'interactive','polyselect');
STREAMobj/modify with the option ‘interactive’ – ‘polyselect’

Do you think the modify function is useful. Are you missing something? Let me know.

Find your way through steep terrain

Posted on

One of my favorite functions in the Image Processing Toolbox is graydist. Actually, it was by my request that Mathwork’s Image Processing development team included this function together with bwdistgeodesic in release 2011b (Thanks for this!). The function help says:

graydist(A,mask) computes the gray-weighted distance transform of the grayscale image A. Locations where mask is true are seed locations.

I am using this function frequently, e.g., when routing through flat areas or topographic depressions when calculating a FLOWobj. However, there is a lot more to do with graydist. One nice application is to find a least cost path over a surface. In this post, I will show how to use TopoToolbox to find a way through steep terrain that minimizes the cumulated gradient (up and down) along the path. Thus, this function will be indispensible for planning your next hiking trip into the wild.

Ok, let’s see how to calculate least cost paths using Matlab and TopoToolbox 2. First, load the DEM shipped with TopoToolbox.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');


Then calculate a gradient grid which we will later use as weights into the distance calculation.

G = gradient8(DEM);

Create a GRIDobj with same spatial referencing as G that we will then use to store the distance grid.

D = GRIDobj(G);

Calculate the gray-weighted distance transform seeded at the lower left corner of the DEM

D.Z = graydist(G.Z,DEM.size(1));
imagesc(D); colorbar


Use the distance grid as input into the FLOWobj constructor function

FD = FLOWobj(D);

Calculate the influencemap seeded at the upper right corner

I  = influencemap(FD,prod(DEM.size)-DEM.size(1)+1);

Take this influencemap as stream grid for deriving a STREAMobj

S  = STREAMobj(FD,I);

and plot it

hold on


Ok, that would be a nice walk along the river. Do you have applications in which you need least cost paths? Let me know.