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.

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

Multi-scale topographic position visualizations

I thought that people would enjoy this beautiful map that I have been working on this holiday season. It is a visualization derived from an SRTM digital elevation model based on a multi-scale topographic position analysis. This is one that you really have to click to enlarge to fully appreciate. I have spent hours lost in the detailed galactic colouring of this map!

This beautiful map of eastern Canada and the US was made with Whitebox GAT’s newest multi-scale topographic position tools. (click on image to enlarge)

This beautiful map of eastern Canada and the US was made with Whitebox GAT’s newest multi-scale topographic position tools. (click on image to enlarge)

I’m working on a paper right now in which I describe a new form of local topographic position analysis and the map above is one of the resulting visualizations. It shows prominent features at the small (blue channel), medium (green channel), and large (red channel) scales. A prominent feature is one that is either significantly above or below the surrounding landscape at the scale of interest. There’s actually a fair amount of analysis (and coding!) involved but if you’re really interested and can’t wait for me to finish the paper, send me an email. Here is a similar map but covering parts of British Columbia, Canada:

Mulit-scale topographic position for British Columbia, Canada (click to enlarge)

Mulit-scale topographic position for British Columbia, Canada (click to enlarge)

And this is the multi-scale topographic position map derived from SRTM data for the Northern Territory of Australia:

Topographic position map of the Northern Territory, Australia (click to enlarge)

Topographic position map of the Northern Territory, Australia (click to enlarge)

Personally, I think that these visualizations are remarkable for their ability to characterize the structure of the surface geology of a region but I’m sure that there are many other interesting applications as well. Interpreting the images takes a bit of experience but the following interpretation key can help:

Interpretation Key

Interpretation Key

Leave your comments below and, as always, best wishes and happy geoprocessing.

****UPDATE (May 27, 2015)****

This work was eventually written up as a manuscript and has recently been accepted for publication by the journal Geomorphology. The citation is:

Lindsay, J.B., Cockburn, J.M.H., and Russell, H.A.J. In press. ‘An integral image approach to performing multi-scale topographic position analysis’ Geomorphology, 245, DOI: 10.1016/j.geomorph.2015.05.025.

This article can be downloaded for free until July 18, 2015 from the following link:

http://authors.elsevier.com/a/1R6Uf,3sl3TsZi

And the submitted draft is available here.

Displaying LAS LiDAR point clouds in the Whitebox map area

I enjoy working with LiDAR data whenever I can because of its remarkable topographic detail and unique characteristics. More often then not, I work with LiDAR data interpolated onto a raster grid. Lately, however, I’ve been working with terrestrial laser scanner data and having a means of quickly visualizing the point data itself has become important to my workflows. Did you know that as of the latest release of Whitebox GAT (v. 3.2.1) there is now native support for displaying LAS files, the commonly used standard format for storing LiDAR point clouds, in the map area?

Adding LAS layers to a map area

Adding LAS layers to a map area

The LAS point cloud will be added to the Whitebox map area in the same way that you can overlay other vector or raster geospatial data. Here’s an example of a LAS dataset overlaid on top of a raster hillshade image.

Example LAS point cloud

Example LAS point cloud (click to enlarge)

The display properties, including the point size, colour palette, and display value ranges can be easily manipulated.

LAS display properties

LAS display properties

You can even render the point cloud based on elevation, intensity, class value, scan angle, or the GPS time (This is a newly added feature that will be present in the next public release).

Rendering attribute options

Rendering attribute options

Here’s an example of a LAS point cloud rendered using its intensity data. It’s interesting how much it resembles a fine-resolution orthophotograph, but this image is actually made up of millions of individual points. You need to zoom in to see them all.

LAS file rendered by intensity

LAS file rendered by intensity (click to enlarge)

Visualizing the GPS time and scan angle can be very useful for identifying individual flight lines in airborne LiDAR datasets. Here’s an example of rendering the LAS file displayed above using the GPS time as the display attribute:

GPS time rendering

GPS time rendering (click to enlarge)

It may not be much to look at (unless you really like orange), but the edges of the various overlapping flight lines are quite apparent.

I’ve also recently added a new LiDAR Histogram tool that will allow you to visualize the statistical distribution of elevation, intensity, or scan angle within a LAS dataset, including outputting a table of the percentiles.

LiDAR histogram tool

LiDAR histogram tool (click to enlarge)

Whitebox has certainly become a first-rate open-source platform for manipulating and analyzing LiDAR data. And it’s getting better with every release! Leave your comments below and, as always, best wishes and happy geoprocessing.

If you enjoyed this blog, you may also enjoy “Working with LiDAR data in Whitebox GAT“.

The Nile River Basin from SRTM data

Someone asked the other day whether the cross-platform, free and open-source GIS Whitebox GAT can handle watershed delineation from massive, regional-scale DEMs. They had a particular interest in the Nile River basin. Heck, if you’re going to go big, why not go huge, right? So I decided to give it a try. First off, I used the Retrieve SRTM Data tool to download the approximately 800 SRTM 3-arcsecond (~90 m) tiles that make up the Nile River basin. This required some experimentation because my first attempt at doing so hit the boundary of the basin and I had to give it a second try. The tool downloaded each of the tiles and mosaicked them into a single large DEM. The final DEM was 45,601 rows by 25,201 columns (a little over 1.1 billion grid cells) and was 4.28 GB in size. I then used the new Breach Depressions (Fast) tool to hydrologically pre-process the DEM by removing artifact topographic depressions and flat areas (i.e. cells with no downslope neighbours). I used a D8 flow algorithm to calculate flow directions, perform flow accumulation, trace the flowpaths issuing from Lake Victoria (the White Nile) and Lake Tana (the Blue Nile), and lastly, to delineate the watershed. The result was this map:

The Nile River Basin (click to enlarge)

The Nile River Basin (click to enlarge)

To fully appreciate this amazing map, you need to enlarge it. Just to put a bit of perspective on the scale of this analysis, take a look at this one:

Nile River World Map (click to enlarge)

Nile River World Map (click to enlarge)

All of Europe has an area of approximately 10,180,000 km2 and the Nile River basin has an area of 3,400,000 km2. That is truly vast.

For my initial attempt, the one in which I truncated the watershed, I used my 13 inch Macbook Pro (2.8 GHz dual-core i7, 16 GB RAM, SSD). When I expanded the area, I also moved to my workstation (3.0 GHz 8-core Xeon, 64 GB RAM, SSD) just to speed up the process a little. I even extracted a long-profile for the White and Blue Nile, although I should have converted the distance units to metres:

Long profile for the Nile River, extracted from SRTM 3-arcsecond data (click to enlarge)

Long profile for the Nile River, extracted from SRTM 3-arcsecond data (click to enlarge).

It was the longest river that I have ever plotted a long profile for; of course, it is the longest river so I guess you can’t get much larger than that! I probably should have extracted the river network using a dispersive flow algorithm like Tarboton’s excellent D-infinity, since there are places where the river bifurcates (i.e. the river course splits), even before the delta. Nonetheless, I’m quite pleased with the result. In fact, I was quite surprised at how well the river course, extracted from the 90 m resolution SRTM DEM data, matched a mapped Nile River shapefile that I located:

Mapped vs. extracted river (click to enlarge)

Mapped vs. extracted river (click to enlarge)

Leave your comments below and, as always, best wishes and happy geoprocessing.

Creating beautiful relief models

Recently I stumbled across a couple of very interesting questions over on the GIS Stack Exchange. One of them was about the use of appropriate colour palettes for displaying DEMs, and the other was about colour palettes for hillshade images. This last one in particular got me thinking a lot about effective cartographic display of terrain and introduced me to the work of the famous Swiss cartographer Eduard Imhof. Why is it that every hillshade image that I have ever created is displayed in a greyscale palette?

Greyscale hillshade of the Banff area.

Greyscale hillshade of the Banff area. (Click to enlarge)

I’ve always thought that hillshading is a very effective way of visualizing topographic information, but greyscale hillshade images by themselves are so boring! Of course the answer is that I don’t normally just display a hillshade image. Instead, the hillshade image is displayed transparently over top a DEM, with colour used to relay elevation information and tone used for illumination. This type of composite relief model is very effective for visualizing terrain. But what about using other colours for illumination other than simply greyscale hillshading, like Imhof did so long ago? To do this effectively in Whitebox, you may want to create a custom palette using the Palette Manager (under the Tools menu), one of the most under-used tools in Whitebox:

The Palette Manager tool in Whitebox GAT

The Palette Manager tool in Whitebox GAT

So here is a hillshade image displayed using the custom palette above [simply blending RGB(0, 50, 100) to RGB(255, 240, 170)]:

Imhof inspired hillshade

Imhof inspired hillshade. (Click to enlarge)

I think that’s beautiful. It’s like bathing the landscape in warm, bright sunlight. Of course, we can actually take this a step further. My father is an artist and retired art teacher. I remember when I was growing up, he taught me that the atmosphere is actually coloured when you have enough air thickness; it’s a little blue. That’s why in this painting of his, one of my favourites, you can see that the hills in the distance become slightly bluer:

Notice the blue ting for distant mountains

Notice the blue tinge for distant mountains (Powassan landscape by James M. Lindsay)

I now know that this phenomenon is called Rayleigh scattering (the reason that the sky is blue), and koodos to my Dad for inherently understanding this phenomenon. So since the height of the air column and density of the atmosphere are greater in the valley bottoms, there should be a ‘blue tinge’ in the valleys. We can create this ‘blue tinge’ by transparently overlaying the DEM, displayed in a light-blue-to-white palette (again created using the Palette Manager), and get this:

Final hillshade image

Final hillshade image. (Click to enlarge)

You can play around with the palettes and the levels of transparency (and display minimum and maximum values) until you get things just the way that you want. It’s almost artistic! I think we can all agree however that this final hillshade is a much more visually pleasing and a more effective display of terrain compared to that initial greyscale image. It looks less like a photo taken of the harsh lunar landscape and more like an earthscape. So have some fun with it and don’t always except the default values of things. And terrain displayed in this way lends itself well to overlaying other information:

Columbia Icefield map. (Click to enlarge)

Columbia Icefield map. (Click to enlarge)

Leave your comments below and, as always, best wishes and happy geoprocessing.

John

SRTM DEM data downloader in Whitebox GAT

I just finished developing a tool for Whitebox GAT that will automatically download Shuttle Radar Topography Mission (SRTM) digital elevation models (DEMs) from the USGS SRTM FTP site. SRTM-3 data are among the best global elevation data, with a grid resolution of approximately 90 m. In many areas SRTM data provide the only topographic data set available. Within the United States, the SRTM-1 dataset provides an improved 30 m resolution. Not only does this Whitebox tool retrieve the SRTM tiles contained within the bounding box of a specified area of interest, but it will also import the tiles to Whitebox GAT, fill missing data holes (which are common with SRTM data in rugged terrain) and mosaic the tiles.

Whitebox's new Retrieve SRTM Data tool

Whitebox’s new Retrieve SRTM Data tool

There have been many times in the past when I have needed to download numerous SRTM files, import the files, fill the missing data holes, and finally mosaic the individual tiles. It can be a laborious workflow indeed. This tool will save me a great deal of time, and so I’m rather excited about it. It’s as though the data magically appear in Whitebox!

SRTM data in Whitebox GAT

SRTM data in Whitebox GAT

Of course, with Whitebox GAT’s extensive Terrain Analysis and Hydrological Analysis toolboxes, there’s plenty of interesting things that you can do once those data do magically appear.

90 m DEM of the entire British Isles

90 m DEM of the entire British Isles. (Note that the image is of a coarser resolution than the actual DEM.)

As an example, I created the SRTM-3 90 m DEM of the entire British Isles shown above in 4 minutes, 56.27 seconds, including the time to download 91 individual SRTM tiles, fill their missing data gaps, and mosaic the tiles. My son was even watching Netflix during my downloading, so I can only imagine how much that slowed things down! I only wish that other data providers could follow a similar data sharing model as the USGS and use an anonymous FTP server to distribute their data. If that were the case, we could have many other data sets automatically ingested directly into Whitebox. I’ll release this new SRTM retrieval tool in the next public release of Whitebox GAT (v 3.2.1), which will likely be sometime later this summer. If you are as keen to try it out as I am, email me for a preview copy. If you have any other comments or feedback, please leave them in the comments section below. As always, best wishes and happy geoprocessing.

John Lindsay

*****EDIT*****

Version 3.2.1 of Whitebox has now been released, with the SRTM downloader tool embedded in the Data Layers menu. Please let me know if you have any issues in using the tool.

****EDIT****

I’ve included an image of the DEM of Ireland for John below:

SRTM DEM of Ireland

SRTM DEM of Ireland (Click to enlarge)