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); imageschs(DEM,log(A),'ticklabels','nice')
For comparison, the single flow direction algorithm returns following flow patterns
FD = FLOWobj(DEM,'preprocess','carve'); As = flowacc(FD); imageschs(DEM,log(As),'ticklabels','nice')
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.
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.
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');
Or call the function modify in this way to select a stream network based on a polygon outline.
Sm = modify(S,'interactive','polyselect');
Do you think the modify function is useful. Are you missing something? Let me know.
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'); imageschs(DEM)
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
imageschs(DEM); hold on plot(S,'w','LineWidth',2)
Ok, that would be a nice walk along the river. Do you have applications in which you need least cost paths? Let me know.
In my last posts I wrote about carving and filling. As noted, both algorithms are implemented in the constructor function FLOWobj and thus are primarily used to derive flow directions through topographic depressions. These options, thus, don’t change the DEM itself. In this post, I will show how you will carve and fill the DEM.
Filling is easy and implemented in most terrain analysis software. A number of different algorithms exist that perform flood-fill operations; in our paper, Dirk and I compared processing speeds to fill an image. We found that the speed at which one and the same DEM is filled varies widely between different software. The Matlab’s image processing toolbox’s function imfill (the function that is called by TopoToolbox function fillsinks) outperforms other software by 8-10 times faster processing speed. However, I like to stress that speed is not the only indicator of performance and that other software may have better ways to manage memory. However, I don’t want to go into detail here…
Here is how to fill a DEM (the Big Tujunga DEM is included in the distribution of TopoToolbox 2)
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); DEMf = fillsinks(DEM);
One way to judge whether filling is a good option is to calculate the “energy” required to modify the DEM. When we fill the DEM, the energy is the total volume added to the DEM.
sum(DEMf.Z(:)-DEM.Z(:))*(DEM.cellsize)^2 ans = 18801000
A volume of ~19×106 m3 is required to force all flow towards the DEM edges. So how about carving? This requires us to calculate a FLOWobj and use it subsequently in a function called imposemin. Imposemin takes a FLOWobj and a GRIDobj and calculates a downstream minimum operator, i.e., it lets minima flow downstream.
FD = FLOWobj(DEM,'preprocess','carve'); DEMc = imposemin(FD,DEM); sum(DEM.Z(:)-DEMc.Z(:))*(DEM.cellsize)^2 ans = 7848000
Carving the DEM requires removing a volume of ~8×106 m3 only. This is less than 1/2 of the energy that filling needs.
Based on these results, I would better carve the DEM to make it hydrologically correct. And, as John Lindsay pointed out in his comment on a previous post, carving is probably the better method in most situations. There may be, however, situations where filling might prove better, e.g., when the DEM has erroneously low valued, single scattered pixels. I had this problem once when working on a DEM in the Himalayas. This is why you find another function in TopoToolbox called elevateminima. This one allows you to raise single pixels to the height of their surrounding pixels.
Another note: Hybrid methods exist that partially fill and breach. I have not implemented these algorithms until now but perhaps I’ll do so in the future.
In my last post, I wrote about the algorithm that FLOWobj uses to derive flow directions in flat areas when choosing the option ‘fill’. This option ensures that flow paths run along centerlines of flat areas – an effect that may be desirable when there are large lakes in the DEM. The other option available is ‘carve’. So what does this option make different?
The preprocess option ‘carve’ works similarly to the ‘fill’ option. The difference, however, is that FLOWobj won’t seek the best path along the centerlines of flat sections. Instead it tries to find the path through sinks that runs along the deepest path, e.g. the valley bottom.
The algorithm works as follows: First, all sinks are identified. The sinks are filled by a flood fill operation (function fillsinks) and the difference between the filled and the original DEM is calculated. Then, the complement (e.g. max(z)-z) is calculated within each sink so that we have a raster where in each sink the smallest values are found in the valley bottoms. Then, this grid is used as weight in a gray-weighted distance transform (function graydist) with the sills of the sinks are used as distance seeds.
Understood? Well, following code and figures might help:
First, let’s load an example DEM and extract a subset that includes a sink in the channel network
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); DEM = filter(DEM,'mean',[13 13]); Z = crop(DEM,[352089 370106]); imagesc(Z); title('DEM')
Identify sink and sill
[I,SI] = identifyflats(fillsinks(Z)); imagesc(I+2.*SI); title('Sill and sink')
Create the cost matrix for the gray-weighted distance transform
C = fillsinks(Z)-Z; C = max(C)-C; C.Z(~I.Z & ~SI.Z) = inf; imagesc(C); title('Complement of sink depth')
Calculate the gray-weighted distance transform
D = GRIDobj(C); D.Z = graydist(C.Z,SI.Z); imagesc(D); title('Gray-weighted distance transform')
The code snippets above are implemented in similar form in the constructor function FLOWobj. Let’s look at a stream network derived for above DEM (Z).
FD = FLOWobj(Z,'preprocess','carve'); imagesc(Z); S = STREAMobj(FD,'minarea',5); hold on plot(S,'k') hold off title('DEM and stream network')
In the last figure, you will notice that the streams tend to run along the deepest cells of the DEM. This is exactly what we wanted to have and what we will need to carve the DEM with a minimum of energy, e.g., with a minimum modification of the DEM. In the next post, I will show how to use the FLOWobj derived with the ‘carve’-option to carve the DEM.
By the way, I have published a paper in ESPL in which I explain carving in more detail (e.g., gray-weighted distance transforms) and how you can use this approach to derive flow patterns in lowland, agricultural landscapes. If you are interested, please read on here.