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.
this is my first post and I want to quickly provide an overview on what you can expect from my posts. I am geoscientist with a major background in geomorphology, GIS, geostatistics and numerical modeling. In the last couple of months Dirk Scherler (Caltech) and I have developed TopoToolbox, a function library for Digital Elevation Model (DEM) analysis in MATLAB. This toolbox has been applied in research topics such as pedometrics, glaciology, landscape evolution modeling and hydrology by researchers and practitioners around the globe. TopoToolbox is also interesting for those amongst you, who study geomorphology, geomorphometry and hydrology and who seek to develop their own models and experimenting with DEMs.
The key to understanding how TopoToolbox works requires some knowledge about mathematics and image processing. I will try my best to cover these topics in posts that should be understood with only little background in these disciplines. Since my interests also go beyond the technical issues I will cover other topics related to DEM analysis as well.
So, hope you enjoy reading. Please leave comments.
Schwanghart, W., Scherler, D. (2014): TopoToolbox 2 – MATLAB-based software for topographic analysis and modeling in Earth surface sciences. Earth Surface Dynamics, 2, 1-7. [DOI: 10.5194/esurf-2-1-2014, Discussion paper: 10.5194/esurfd-1-261-2013]. open access and software available here [DOI: 10.1594/IEDA/100175].
Schwanghart, W., Kuhn, N. J. (2010): TopoToolbox: a set of Matlab functions for topographic analysis. Environmental Modelling & Software, 25, 770-781. [DOI: 10.1016/j.envsoft.2009.12.002]
Download the latest TopoToolbox Version at the Community Surface Dynamics Modeling System (CSDMS) code repository.
TopoToolbox provides a set of Matlab functions that support the analysis of relief and flow pathways in digital elevation models. The major aim of TopoToolbox is to offer stable and efficient analytical GIS utilities in a non-GIS environment in order to support the simultaneous application of GIS-specific and other quantitative methods.
TopoToolbox’ algorithms are largely based on sparse matrix algebra and image processing techniques. Each function is implemented in Matlab-code (*.m-files) and can be easily modified and adapted to specific needs of the user.
TopoToolbox can be regarded as a construction kit and new functions will be added in the future. If you have ideas or comments on how to complement or improve TopoToolbox, if you like to contribute, or if you have any questions or remarks, please contact me: