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]
Today’s guest post comes from Philippe Steer, assistant professor at the University of Rennes. Philippe contributed two new functions and I asked him to describe them here to show potential applications. Thanks, Philippe, for this contribution!
I recently needed to compute river network orientation. Indeed river network orientation can reveal distributed tectonics deformation at the vicinity of strike-slip faults (e.g. Castelltort et al., 2012) or can enable to identify developing vs paleo river networks. One way to obtain orientation is first to extract the river segments that are linking confluences to outlets or confluences to channel heads and then to compute the orientation of the resulting segments. The streampoi function allows to extract these feature points of the network such as channel heads, confluences, b-confluences (the points just upstream any confluence) and outlets. To connect these feature points into segment, one can use the propriety that each segment belongs to a unique watershed that is given by the function drainagebasins using as outlets the true river outlets and the b-confluences. By definition, each segment also corresponds to a unique value of strahler order.
Once the river segments are identified, segment orientation is easily obtained. What is interesting is that defining these segments also paves the way to compute several other geometric metrics averaged over the river segment such as segment direct or flowpath length, sinuosity, direct or flowpath slope, mean drainage area and even steepness index (if assuming a value of concavity).
To cleanly do the job, I have written two MATLAB functions:
Here is a piece of code showing you the use of these functions (assuming you have already loaded the DEM of interest):
%% Flow direction and drainage area DEMf = fillsinks(DEM); FD= FLOWobj(DEMf); A = flowacc(FD); %% Extract rivers % Set the critical drainage area A_river=1e7/(res^2); % m^2 W = A>A_river; % Extract the wanted network S = STREAMobj(FD,W); %% River segments % Extract segment and compute their orientations mnratio=0.5; [segment]=networksegment(S,FD,DEM,A,mnratio); % Plot results plotsegmentgeometry(S,segment)
To illustrate the results, I thought the South Island of New-Zealand was a good example. Indeed, as first suggested by Castelltort et al. (2012) the river catchments at the south of the Alpine fault are tilted, possibly suggesting that they have recorded a distributed shear deformation. This tilting is also quite obvious when looking at the river segments in the following figure. Another interesting feature is that the river segments north of the Alpine fault, facing the tilted rivers, display a higher steepness than their south neighbours. This could be due to a long list of forcings, including active tectonics, orographic precipitation gradient or the transition from glacial to fluvial erosion.
Castelltort, S et al. (2012). River drainage patterns in the New Zealand Alps primarily controlled by plate tectonic strain. Nature Geoscience, 5(10), 744-748. [DOI: 10.1038/ngeo1582]
*** 28 Apr 2016: Please note that the weak performance of the new MATLAB JIT is not an issue anymore. Please see comments below. All functions were updated and so I have no reservations using Matlab versions of 2015b and later ***
I like the fact that The Mathworks releases a new version of MATLAB every 6 months. The first thing I usually do when receiving the latest version is to read the release notes. This time, two new features grabbed my attention. First, the introduction of graph and network algorithms is something that I missed so far. Graph theory algorithms (though actually some are hidden in sparse matrix algorithms such as dmperm) were lacking until 2015b (unless you had access to the Bioinformatics Toolbox). I found this lack of graph theory algorithms in the core product strange given that matrix linear algebra and graph theory are quite related to each other. In fact, graph theory is integral to image processing and digital terrain analysis, and thus I am delighted to see this update. I am pretty sure the new algorithms will somehow be incorporated in TopoToolbox sooner or later. These are the good news.
The second feature comprises changes to the JIT (just-in-time) accelerator. Many of you probably have never bothered about this feature, but I did. JIT compilers compile parts of the code at runtime and thus increase the speed at which it is executed. Why did I bother? Well, some of the code at the heart of all flow-related calculations is not vectorized. You’ll find a little chunk of code in nearly each FLOWobj method that loops through the entire stack of indices stored in the FLOWobj. Yes, I know that in many situation this can be vectorized and I had done so before (see the M matrix in TopoToolbox 1). However, letting backslash do the work proved in most situations worse (both concerning memory and computation time) than wisely permuting the matrix and then let loops do the rest. Ok, let numbers speak and show why the JIT update is possibly bad news. I use following piece of code to evaluate computation times in R2014b and R2015b:
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); FD = FLOWobj(DEM,'preprocess','carve'); tic; A = flowacc(FD,ones(DEM.size),ones(DEM.size)); toc
In 2014b, this yields:
Elapsed time is 0.072537 seconds.
In 2015b, it shows:
Elapsed time is 0.420721 seconds.
No worries, I repeated these calculations several times and took the fastest ones of each to print here. Whatsoever, 2015b needs ~6 times longer to evaluate the code. That are bad news. These findings prompted me to run other tests with gradient and FLOWobj. And here, I did find either no or some slight improvement with R2015b.
Ok, what’s next. Well, I cannot accept that the JIT turned so bad for my looped code chunks. The release notes say:
These performance improvements result in significantly faster execution of many MATLAB programs with an average speed-up of 40% among 76 performance-sensitive applications from users. Of these tested applications, 13 ran at least twice as fast and only 1 slowed down by more than 10%.
Am I among the 1.31% (1 out of 76)? Or are there strategies to maintain the speed of TopoToolbox? I consider following: 1. Place these loops in extra functions or subfunctions. 2. Check whether the problem is related to the fact that the indices are stored in an object (e.g., object properties). Before I find the right strategy I will definitely keep both R2014b and R2015b on my computer, and I recommend you do so as well.
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 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.
Image Posted on Updated on
People using TopoToolbox asked me quite frequently about the difference between the two preprocessing options when calculating a flow direction object (FLOWobj): ‘carve’ and ‘fill’.
First and most important: These options do not alter the DEM, they only control how FLOWobj internally preprocesses the DEM to optimize how flow is routed through topographic depressions.
In this first post I’ll describe how the option ‘fill’ works.
Option ‘fill’: The constructor function FLOWobj fills all topographic depressions (function fillsinks). It then identifies all flat sections (function identifyflats) and calculates an euclidean distance transform (function bwdist) within each flat section seeded at their outer rim. Subsequently, the complement of the distance transform is calculated in each flat area. The complement grid is subsequently used as weights in a gray-weighted distance transform (function graydist). This generates an auxiliary surface in each flat area. The surface has the nice property that flow will run through the centerline of flat areas.
So how does this approach to derive flow directions in flat areas compare to other software. Above figure shows how other software handle flow directions in the flat areas. It seems that only TauDEM uses a similar approach as TopoToolbox 2. Other software seem to handle the problem differently and produce linear, somewhat unelegant flowpaths through flat sections.