Modelling the spatial pattern of potential impoundment size from DEMs

For much of my research career I’ve been deeply interested in the techniques that are used to model surface and near-surface drainage patterns using DEMs. Much of this work has focused on DEM processing techniques for handling sinks, such as topographic depressions. Removal of depressions, particularly artifact depressions, from DEMs is a necessary step for most flow-path based hydrological analysis and there is a great deal of research dedicated to developing and testing the methods for this type of flow enforcement operation.

For a long time now, however, I’ve been thinking about the opposite problem. Instead of removing depressions from DEMs, how can we add depressions? In fact, artificially damming a site in a DEM is fairly straight forward and simply involves raising the elevations along a dam site. To me, the more interesting problem is how to locate sites in a DEM that, were a dam or dyke of a specified size inserted, would produce an impoundment (i.e. the artificial body of water behind the dam) of a certain size. In other words, can we measure the impoundment size for each grid cell in a DEM based on specified maximum dam/dyke height and length? It turns out that this is a fairly tricky problem which is why I’ve been interested in it for a couple of years now. An impoundment size raster would be very useful for all sorts of applications. For example, it could be used for siting potential wetland restoration projects, monitoring wetland drainage at the landscape scale (presumably using LiDAR DEMs), modelling beaver habitat (or other related aquatic species), siting tailings facilities, hydroelectric dam siting, and many other applications.

Well I’ve finally figured it out. I present to you Whitebox GAT’s new Impoundment Index tool:

Impoundment Index Tool

Impoundment Index Tool

I feel like the image above should have been slowly revealed with an increasing-tempo drum-line and trumpeted crescendo in the background but sadly these are the limitations of the media. Perhaps you could go to the top of your screen again and slowly scroll downward while humming.

It only took two years to ruminate over this little problem and, once I figured it out, an afternoon to code. In fairness, I’ve had a lot of other things to keep me distracted in that time, like developing an entire GIS for example. For a specified input DEM, the Impoundment Index tool will calculate the size, either the area or the volume, of the impoundment the would result from inserting a dam of a specified maximum height and length. This is a sample of what the output raster looks like:

Sample impoundment index raster

Sample impoundment index raster

As you might expect, the impoundment size is highest in sites with locally convergent topography and extensive relatively low-lying upslope areas. In this case, each grid cell in the raster contains the estimated area in square metres of the impoundment that would result from inserting a 5 m (maximum) high 11 grid cell long dam.

There are two parts to the algorithm. The first component measures the height of a potential dam for each grid cell based on the topographic profile perpendicular to the flow line at the site. The second component to the algorithm performs a type of flow accumulation operation. Where most accumulation operations of this type propagate a single value (e.g. the area upslope) along flow-paths, this algorithm propagates an entire elevation distribution. In that way, it is possible to measure, for each grid cell, the number of grid cells in its catchment that are less than the previously measured maximum dam crest elevation. The really tricky part is how it manages to propagate the elevation distributions to each cell in a computationally and memory efficient manner such that your computer doesn’t explode!

The Impoundment Index tool can even output a DEM where grid cells with an impoundment size within a specified target range have their dams inserted, allowing you to model the potential damming sites:

Modelled impoundments

Modelled impoundments

The tool is still experimental but will hopefully be officially released in Whitebox GAT v. 3.2.3 or whatever the next update becomes. If you can’t wait to play with the Impoundment Index tool, it’s already in the code repository. So you can easily start messing around with it by selecting ‘Update Scripts From Repository’ under Whitebox’s Tools menu. After updating and restarting Whitebox, you’ll find the new tool within the Wetlands Tools and Flowpath Terrain Attributes toolboxes. I’d like to write a paper describing exactly how the tool works in detail in the hopes that other GIS packages will also implement this useful method. Maybe you’ll enjoy this new and exciting function as much as I’ve enjoyed the process of creating it. As always, best wishes and happy geoprocessing.

Advertisements

Whitebox Geospatial Analysis Tools v. 3.2.2 Released

It is with tremendous pleasure that I announce the release of the latest version of Whitebox Geospatial Analysis Tools, version 3.2.2. It has been several months since the release of the previous version of Whitebox GAT and we have incorporated several significant new features, tools, and bug fixes. The following is a partial list of changes:

  • LAS files can now be rendered based on the Classification, Scan Angle, and GPS Time in addition to Elevation and Intensity. This is on top of several other improvements that have been made to the display of LAS file data in Whitebox GAT.
  • Added Conditional Evaluation tool (if-then-else statement). I really like this tool and now wonder how it was that I managed to get by without it.
  • Added a PickFromList tool, which outputs a raster with the value from a stack
    of rasters based on an input position raster.
  • Added LowestPosition, HighestPosition, PercentEqualTo, PercentGreaterThan, and PercentLessThan tools for working with raster stacks, i.e. lists of overlapping
    raster images.
  • Added a tool to create histograms based on the elevations, intensity, and
    scan angles within a LiDAR (LAS) file. It will also output the percentiles of
    the distribution, e.g. 95th percentile of elevation.
  • Added the AddPointCoordinateToTable tool, which can be used to add the x-y
    coordinates of points within a Point type ShapeFile as fields within its attribute table.
  • Added a tool to filter the points in a LiDAR (LAS) file based on a threshold
    in the absolute scan angle. Currently the output is a shapefile of mass-point shape but eventually we would like to have it write to a new LAS file.
  • The Merge Points Files tool was replaced with the more general Merge Shapefiles
    tool, which works with any ShapeType.
  • Added the FindLowestHighestLocations tool, which will output vector points
    corresponding to the lowest and highest points in a raster. This has already come in handy several times.
  • Added ExtractRasterValuesAtPoints tool for extracting the cell values of each
    image in a raster stack corresponding to a set of vector points. Values are
    placed within the vector attribute table.
  • Added DeleteSmallLakesAndExtendRivers tool, which can be used to remove small lakes (polygons) in a vector drainage network and to extend the associated river
    network (intersecting polylines) into the interior of the lake. I created this tool in response to an interesting question asked over on the GIS Stack Exchange.
  • Added a Long Profile From Point tool, which can generate one or more longitudinal
    profiles for the flowpaths issuing from a set of vector points.
  • Modified the Mosaic With Feathering tool to handle RGB images in addition to
    continuous scale rasters. At the moment, this only works for the nearest-
    neighbour mode. I’m not sure why I didn’t do this earlier.
  • Added Image Stack Profile tool to create a line graphs for a set of vector points
    through an image stack, such as multi-spectral satellite image data. This can be handy for visualizing the spectral signatures of individual pixels.
  • Added a Simple Region Grow tool that will perform a very simple region grow
    segmentation of pixels in an image stack based on a specified threshold. I’d like to continue development in this area and eventually include a full object-based image segmentation.
  • Added parallel implementations of the D8 flow pointer and accumulation algorithms. At the moment this is really an experimental tool that is not intended for widespread use but there is more to come, including parallel versions of all the flow accumulation algorithms.
  • Fixed a bug in the Hillslopes tool.
  • I’ve added a whole suite of tools to the Elevation Residuals toolbox for
    performing multi-scaled topographic position analysis. This includes modified
    tools for calculating difference and deviation from mean elevations using an
    integral image approach that is extremely computationally efficient, even with
    large search windows. It also includes the Local Topographic Position Scale
    Signature and the Maximum Elevation Deviation tools. Combined with the
    ‘customRGBScript.groovy’ script, these functions allow for the creation of
    some spectacular multi-scale topographic position visualizations. See this link
    for examples.

And like a proud father, I can’t resist showing some nice pictures…

Whitebox GAT 3.2.2 screenshot

Whitebox GAT 3.2.2 screenshot

Whitebox GAT 3.2.2 screenshot

Whitebox GAT 3.2.2 screenshot

Whitebox GAT v. 3.2.2

Whitebox GAT v. 3.2.2

One point of note is that Oracle stopped supporting Java on Windows XP some time ago and therefore recent versions of Whitebox GAT no longer function on this platform. It’s time to upgrade!

Please let me know if you have any feedback or questions regarding the new version of Whitebox GAT and I hope you enjoy all the new goodies. As always, best wishes and happy geoprocessing.

UPDATE (April 15, 2015): It would seem that there were some issues with some of the plugin tools written as Groovy scripts that resulted from breaking changes associated with the update from Groovy 2.3 to Groovy 2.4.1. To overcome these issues I regressed the Groovy library linked to by Whitebox GAT to the 2.3 version and now all of the affected tools are working properly. Thank you to the users that alerted me to this issue.

The language of sinks, topographic depressions, and pits…

Over the last several decades, there has been a tremendous confusion over the terminology used to describe the features that are involved in flow-enforcement of digital elevation models (DEMs) applied in surface drainage pattern modelling. In particular, the terms sink, depression, and pit have been used interchangeably and this has caused, in my opinion, a great deal of confusion in both the academic literature and practice. In a manuscript that I am currently drafting, I developed a typology of features involved in the flow-enforcement step, a necessary first-step for any process involving hydrological analysis of topographic data. The following flow-chart describes the relation between these overlapping concepts in what I think is a consistent manner.

A typology of features involved in flow-enforcement

A typology of features involved in flow-enforcement