Some methods for deriving variables from digital elevation models for the purpose of analysis, partitioning of terrain and providing decision support for what-if scenarios

By

Peter Caccetta

CSIRO Mathematical and Information Science (CMIS)

email: Peter.Caccetta@csiro.au

24/02/1999

Background

Landform is a significant factor in agricultural production and land management, influencing the relative growth of crops within a paddock as well as larger-scale processes such as salinity, waterlogging and erosion. The potential for using digital elevation models (dems) in the analysis of landscape processes appears large.

Previous research by members of CMIS, W.A., have concentrated on techniques for deriving variables from dems formed from the interpolation of contour data. The contour data are readily available for much of the south-west of W.A. in elevation sampling intervals of either 5m, 10m or 20m. The accuracy of the dems imposes limitations on the sort of variables which can be reliably derived. Those found to be useful include estimates of (flow) slope, upslope area and percentage upslope cleared area. These variables have been used in the analysis of salinity risk and for partitioning the landscape into landform units such as hilltops, ridges and upper slopes, upper valleys, lower valleys and broad valleys. For a review of these and other approaches, see Caccetta (1997).

The advent of relatively high resolution dems (elevation accurate to 2 metres sampled on a 10m easting/northing grid) allows a greater variety of variables to be derived. For instance, variables such as height above (the nearest feature) the nearest stream or salt-affected area may be determined with a much greater accuracy than from a dem constructed from, say, 10m contour data.

This report describes a number of algorithms, based on principles of water flow, which are relevant to dems. The variables presented are by no means exclusive.

The first algorithm considered preprocesses the height data to obtain hydrologically sound dems. The algorithm is relatively simple to apply, computationally efficient and does not require that streams have been captured by some other means in a pre-specified format (for example in descending order – as in ANUDEM). Not requiring streams as input is considered an advantage as the dems often have greater spatial accuracy than many of the streams captured at 1:50000 or 1:100000, which constitutes much of the region of interest.

Given a hydrologically sound dem, a number of variables are described which can be derived directly from the dem. These include water accumulation, relief energy, average flow path length, average slope gradient, flow slope, height above and distance from. Hypsometric curves, stream order, stream density (km/km²) plan and profile curvatures can also be derived, but are not discussed here.

Next are considered variables which may be derived from combinations of the above variables. These variables are referred to here as compound variables, and include wetness index, stream power and universal soil loss equation.

Finally, notes on how to use the programs to generate the variables is given in section 5.

1     Preprocessing of DEMs to Remove Unwanted Minimum

Typically dems have many local minima (depressions / pits), many of which are the result of errors in the dem production process as opposed to a few which are real (for example lakes and soaks). The erroneous minima cause difficulty for algorithms which simulate water flowing across the landscape, and methods for removing them need to be employed. The first step in this process requires an expert to provide information to locate all true depressions. An algorithm which eliminates all remaining false depressions can then be applied. A (raster) dem with all erroneous depressions filled and all pixels having a defined flow direction(s) is sometimes referred to as being hydrologically sound. The process of removing spurious minima is commonly referred to as pit filling.

As the name implies, pit filling is the process of changing the relative elevation of the dem for an area which does not flow such that the modified dem does flow. Generally the elevations of the pixels within a pit are increased until they match the elevation of the lowest boundary pixel which does flow. In this case, a flat area is produced.

Flat areas by definition do not flow, so the next step in the procedure requires that flow directions be assigned. Again different algorithms exist here, but perhaps the best known is that proposed by Jensen and Domingue (1988). For this algorithm, pixels in a flat region are recursively assigned to flow to a neighbouring pixel which does have flow defined. A pixel at the edge of a flat region which has a flow direction defined is used as the seed for this recursive process. In this way, the flat region of the dem which results from pit filling is still flat; however, flow models may be applied using the flow direction map (commonly referred to as a flow net). A limitation of this technique, and its variants (for example Peckham (1995) – RIVERTOOLS), is that pixels are defined to flow in one direction only. Algorithms which route flow in one direction are commonly referred to as D8.

One means of minimising the number of pits encountered on a flow path is to force, during the dem interpolation process, the surface to conform to stream line data where they exist. Here the dem is constructed such that it has a relative minimum at the location of streams (that is streams will sit in the bottom of valleys) and that stream elevations are ordered such that water will run from the stream head to the stream base. This is the approach taken in ANUDEM (see Hutchinson, (1989)). This approach requires that each stream vector needs to be captured in the direction of flow; stream data in W.A. do not typically have this property.

Soille and Gratin (1993) proposed an alternative approach to producing a hydrologically sound dem, based on principles taken from the field of mathematical morphology. Their approach combines the pit filling and flow direction specification into a single step. The result of applying the approach is a pit-free relative dem having no (strictly) flat areas; areas which would normally be flat as a result from pit filling have very small elevations added such that flow is defined.

An attractive feature of the algorithm proposed by Soille and Gratin is that it requires no user-defined parameters and runs in O[n] (that is, it is fast), where n is the number of pixels in the dem. After applying the algorithm, the user is guaranteed a model which flows – there are no cases in which it fails to define flow directions. Unlike algorithms such as that proposed by Jensen and Domingue, the approach may also be used with multiple flow algorithms, which have many advantages over single flow models. For discussions relating to single flow algorithms versus multiple flow algorithms, see for example Quinn et al. (1991) and Caccetta (1997) pp 131-136.

Not considered in this document are approaches commonly referred to as process modeling. The characteristic feature of these models is that they require many parameters and boundary conditions to be specified, and in many instances are based on the solution of differential equations. In the context of broad scale mapping, this requirement limits the use of such models to fairly small catchments (or very coarse grids), and typically accurate information regarding the boundary conditions and parameters are not available. Many algorithms based upon process modeling exist. Examples of this form of modeling include MODFLOW (McDonald and Harbaugh, 1988) and TOPOG.

Figure 1  (bottom) Water accumulation algorithm applied to a dem with pits. (top) Water accumulation algorithm applied to a dem with pits removed. Notice that the flow paths (blue) are continuous for the pit-filled version as compared to being discontinuous for the unfilled version.

 

2     Variables derived directly from the dem

Having removed spurious minima from the dem, a number of variables based on overland flow are readily derived. Here multiple flow is assumed and is formulated as given by Quinn et. al. 1991. This is the formulation adopted for TOPMODEL (Beven and Kirkby, 1979); a rainfall-runoff model which bases its predictions on catchment topography, developed by the hydrology and fluid dynamics group in the department of environmental sciences at Lancaster university.

The multiple flow algorithm distributes flow according to the relative slope of downhill pixels. That is, pixels may contribute flow to one or more (sub) catchments. Variables derived using multiple flow algorithms as their basis are thus weighted by the proportion of flow that each constituent pixel contributes to the (sub)catchment in question.

All figures are coloured red, indicating relatively low values , through orange, yellow green, blue and then purple, which signifies relatively high values.

2.1     Upslope area (commonly called water accumulation)

An upslope area image has assigned to each pixel a value corresponding to the catchment area upslope of the pixel. An example of this variable is given in figure 2. From the figure it is observed that hilltops and ridges have relatively low values as compared with valleys which have relatively high values.

Figure 2  Expected upslope area, or commonly referred to as water accumulation.

An example of this variable is depicted in figure 3. From the figure it is observed that break aways and hilltops reminiscent of rocky outcrops appear as purple as compared with other regions.

Figure 3  Average slope of the expected upslope catchment area of a given pixel.

2.3 Average expected relative altitude of the upslope catchment area (or relief potential)

Each pixel is assigned a value equal to the average upslope catchment elevation minus the pixel elevation.

Under gravity, the potential energy of a mass is proportional to the magnitude of its relative elevation.

If converted to kinetic energy, this equates to velocity increasing with relative elevation. If the mass is an idealized confined water column, this equates to column base pressure increasing with the height of the column. It would be of interest to see if this variable is related to water erosion (velocity) or saline seeps (pressure).

An example of this variable is given in figure 4.

Figure 4  Average expected relative altitude of the upslope catchment area of a given pixel.

2.4     Expected flow path length

Each pixel in the image is assigned the upslope (overland) distance that the flow has traveled. An example of this variable is given in figure 5.

Figure 5  This figure depicts the expected flow path length. The variable may be used to identify relative position within a catchment.

2.5     Slope in the direction of flow

The slope of the surface may be calculated in a number of different ways, and depending on which way it is done, pixels will have different values. The slope in the direction of flow is calculated as the weighted average of the slope to all neighbouring pixels with elevation less than the current pixel. In this way, slope for a pixel wide drainage line, say, is the slope that a flow path observes as opposed to slope produced by common difference kernals.

Figure 6  Average slope in the direction of flow.

2.6     The number of neighbouring downhill pixels (or outflow, convergence index)

Figure 7  The number of neighbouring downhill pixels. This variable appears to be of limited value for the dem considered.

2.7     Expected height above (the nearest feature pixel)

Each pixel in the image is assigned its relative elevation as compared to the nearest feature pixel. Here `nearest’ is defined in terms of overland flow. An example of the variable is given in figure 8.

Flow paths were firstly extracted by stratifying the water accumulation image, these flow paths defined the feature pixels. A height above calculation was then performed giving an estimate of the height above the nearest flow path pixel. Existing salt scalds could also have been used as feature pixels, giving the height above the nearest salt scald. This variable may be of use in salinity prediction.

Figure 8  Height above the nearest feature pixel.

2.8     Expected distance from (the nearest feature pixel)

Each pixel is assigned the value of its distance to the `nearest’ feature pixel. Here `nearest’ is defined in terms of overland flow. An example is given in figure 9.

Flow paths were firstly extracted by stratifying the water accumulation image, these flow paths defined the feature pixels. A distance from calculation was then performed giving an estimate of the distance from the nearest flow path pixel. Existing salt scalds could also have been used as feature pixels, giving the distance from the nearest salt scald. This variable may be of use in salinity prediction.

Figure 9  Distance from the feature pixel. Here flow paths, extracted as surrogates of stream lines, were used as features.

3     Compound Topographic Attributes

This section covers some attributes which are commonly quoted in the literature and may be formed by combining various layers from the previous section. Here we merely introduce them by their commonly used names without any consideration of their validity for the concepts implied by their names. We note that different authors use the same concepts but define the indices differently. The use of Insitu measurements may be useful in determining the appropriate definition.

3.1     Wetness index

The wetness index is the log of the ratio of catchment area and flowslope. Flat areas having a large upslope area will have a relatively high wetness index as compared with steep areas having a relatively small catchment area. An example of this variable is given in figure 10.

Figure 10  An example of wetness index.

3.2     Stream power

Stream power is defined as the product of the upslope area and flowslope. Steep areas having a large catchment area will have a relatively high value for stream power as compared with flat areas having relatively small catchments. An example is given in figure 11.

Figure 11  An example of the stream power variable.

4    Some extensions to the above models

The water accumulation model simulates water flowing across the terrain surface. Obviously, many factors including land cover, soil type and contour banks will affect this flow. Schultz (1994) discusses the use of remote sensing as input to hydrogeological monitoring and suggests that remote sensing input can be used both as a basis for model parameter estimation and model input.

A simple application of this is the modification of the water accumulation algorithm to form estimates of upslope cleared area and percentage upslope cleared area (Caccetta, 1997 pp 234-235). This is achieved by initializing all pixels which have remnant vegetation on them to zero (area), and all others to 1 (unit of area) at the commencement of the algorithm; then the map that results is upslope cleared area. Dividing this by upslope area gives the percentage upslope cleared area. In the context of salinity mapping, catchments which have no clearing (zero percent upslope cleared area) would tend to have little risk of salinity, as compared to a catchment which is completely cleared (100% upslope cleared area). It would be of interest to examine the sensitivity of this variable for mid-range values.

We have implemented the water accumulation algorithm with the following two options:

1. The user may specify a 1-byte land cover relative equivalent area map. This map provides an area weighting for each pixel. For instance, in the preceding example, remnant vegetation is given a value of zero and cleared land a value of one. There may also be other things worth doing, for instance a pixel in a 600mm rainfall zone, say, may be given a value of two compared to a value of one for a 300mm zone. Here we merely speculate, to introduce the concept.
2. The user may specify a 1-byte land cover barrier layer. The barrier layer specifies the percentage of flow that is lost to further flow. For example, a dam may be specified as a 100% barrier, whereas a contour bank may act as a barrier some what less that 100%. Again we merely speculate to introduce the concept.

In regards to salinity prediction, if any of the variables listed above have predictive power, then the fast algorithms which perform these calculations may form the basis of a scenario modeling system, for example, for examining the likely effect of planting previously cleared paddocks to plantation, say.

5    How to generate the above results

The variables discussed in this report were generated using a number of C coded routines. All routines assume 2 byte unsigned integer dems. All files are assumed to have most significant byte first (MSBFirst) format (this is the SUN default). The routines are:

• pitfill2 – his program performs the pit filling and creates a model which is required by the routines which derive variables based on overland flow. The program generates two files: (1) a relative elevation file (REV) and (2) a REV support file (SUP). To note here is that the REF file is not 1 to 1 with the original DEM: it is purely a surface which has the required local flow characteristics.

• mmaccum – this program performs the standard upslope area computation. Three options are provided: (1) standard upslope area (2) upslope area with each pixel initialized to a value specified by an input image (3) the same as two plus the specification of an infiltration map.
• mmaccum_plus – this program calculates upslope slope, upslope relative altitude and flow path length.
• habove4 – this program generates a height above and distance from the nearest feature pixel. Feature pixels are specified as non-zero values in a 1 byte input image.

5.1     pitfill2

    Input file: dem, 2 byte unsigned integer.

    Output files: a 2 byte unsigned integer REV file; a 4 byte unsigned integer SUP file.

5.2     mmaccum

    Option 1

Input files: the REV and SUP files produced by pitfill2.

Output file: 4 byte unsigned integer upslope area file.

Input files: the REV and SUP files produced by pitfill2, a 1 byte unsigned char initialisation layer specifying the relative importance of each pixel (each pixel is worth `x’ units of area).

Output file: 4 byte unsigned integer upslope area file.

Input files: the REV and SUP files produced by pitfill2, a 1 byte unsigned char initialisation layer specifying the relative importance of each pixel (each pixel is worth `x’ units of area), a 1 byte unsigned char infiltration layer (each pixel represents the percentage of flow that leaves to pixel).

Output file: 4 byte unsigned integer upslope area file.

Input files: the REV and SUP files produced by pitfill2.

Output file: 2 byte unsigned integer file.

Input files: the REV and SUP files produced by pitfill2, a 1 byte unsigned char feature file. Features are assumed to be those pixels with non-zero values.

Output file: (options 1 and 2) 2 byte unsigned integer file. (option 3) 4 byte unsigned integer file.

References

Beven K.J. and Kirkby M.J. (1979). A physically based variable contributing area model of basin hydrology. Hydrol. Sci. Bull., 24, 43-69.

Caccetta, P.A. (1997), Remote Sensing, Geographical Information Systems (GIS) and Bayesian Knowledge-Based Methods for Monitoring Land Condition, PhD thesis, School of Computing Science, Curtin University of Technology.

Hutchinson, M.F. (1989), A new method for gridding elevation and stream line data with automatic removal of pits, Journal of Hydrology, 106:211-232.

Jensen S. K. and Domingue, J. O. (1988), Extracting topographic structure from digital elevation data for geographic information system analysis, Photogrammetric Engineering and Remote Sensing 54(11):1593-1600.

McDonald, M.C., and Harbaugh, A.W., 1988, A modular three-dimensional finite-difference ground-water flow model: U.S. Geological Survey Techniques of Water-Resources Investigations, book 6, chap. A1, 586 p.

O' Loughlin E.M (1981). Saturation regions in catchments and their relation to soil and topographic properties. J. Hydrol., 53, 229-246.

Peckham, S. D. (1995), Appendix A: River Network Extraction from Large DEMs, Excerpt from the doctoral thesis of Scott Dale Peckham, University of Colorado, Boulder

Qiuinn, P., Beven K., Chevallier, P., and Planchon, O. (1991), The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models, Hydrological Processes, 5:59-79.

Quinn, P., Beven, K., Chevallier, P., and Planchon, O. (1991), The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models. Hydrological Processes, 5:59-79.

Schultz, G. A. (1994), Meso-scale modelling of runoff and water balances using remote sensing and other GIS data, Hydrological Sciences – Journal des Sciences Hydrologiques, 39(2).

Soille, P. and Gratin, C. (1993), Short Communication: An Efficient Algorithm for Drainage Network Extraction on DEMs, Centre de Morphologie Mathematique, Ecole Nationale Superieure des Mines de Paris, 35, rue Saint-Honore F-773005 Fontainebleau Cedex, France.

Return to Remote Sensing and Image Integration