Working In GDAL (WIG) – Creating a Topographic Position Index Map

In this monthly series, I turn my attention to one of the most powerful software packages for raster processing available – the Geospatial Data Abstraction Library (GDAL) – and as an added bonus, it’s free! Each month, I focus on a different utility offered in GDAL, walking the reader through its applications and offering up a series of tips and illustrations so that you can emulate this process in your office and/or home. In this edition of WIG, I show our readers how to use GDAL to create a topographic position index (TPI) map from a raster digital elevation model (DEM).

About GDAL

GDAL is an open source software application that was launched in 1998 and has been updated multiple times since. While the current version of GDAL (i.e. version 1.9) can be downloaded from multiple locations, I find this site to be the most up-to-date and easy to navigate. I suggest downloading the stable version of the MSVC2010 build for either 32-bit or 64-bit PCs depending on which you have access to. You can also download GDAL as part of an easier to install program, FW Tools, but the version of GDAL included with this installer is usually not the most recent.

GDAL contains a wide array of utilities to help you process raster files. I find that GDAL is more stable and runs quicker than ArcGIS, ENVI and PCI for data production – although the toolkit is more limited than in these applications. GDAL does not rate high for user-friendly functionality as it is command line driven (i.e. words and text, not graphics and mouse clicks); and this is the inspiration for creating the WIG series.

What is the Topographic Position Index (TPI)?

The topographic position index is a relative measure of a location’s elevation with respect to its surroundings. It is calculated in GDAL by subtracting the mean elevation of the 8 cells surrounding a central location from the central location’s elevation. Given that TPI is calculated by comparing a central point to the mean, surrounding elevation, there are negative values unlike the terrain roughness index which I wrote about last month. Here is a quick summary of how you can interpret TPI values:

• Zero or near-zero – these are areas that are flat or have a near continuous slope; for example, the pixel is in the middle of a evenly sloping hill
• Large positive – the central pixel is, on average, much higher than the surrounding areas; for example, on the top of a ridge or hill
• Large negative – the central pixel is, on average, much lower than the surrounding areas; for example, at the bottom of a valley or is a gulley

It is also important to note that TPI is scale-dependent in that you can compare a central point to the mean elevation of just the 8 immediate surrounding pixels (as is done in GDAL); or you can compare this pixel, for example, to the 5 closest neighbors in all directions. By alternating the radius of the mean elevation calculation, you can identify either fine or large scale features. In GDAL, you will be locked into identifying fine scale features such as ditches and mountain peaks given that it will only calculated TPI with a 1-pixel radius (i.e. a central point’s 8 immediate neighbors).

Creating a Topographic Position Index Map

In order to create a TRI map from a DEM, I employed the GDALDEM functions – you can find out more about them here. Here are the steps I took to create this layer:

1. I moved my DEM file, a digital terrain model (DTM) produced by DigitalGlobe with 2-meter resolution over Gadara, India, to a folder called, ‘Test,’ located the root of my C: Drive
1. This is the same data used in last month’s edition of WIG
2. I chose a simple folder name and file location (i.e. C:\test) for the data to make things easier when working with the command line interface of GDAL
   
   
2. I also changed the name of the DTM file to ‘DTM.tif’ to make things easier when working with GDAL
   
   
3. After starting up GDAL, I used the command: cd\test
1. This command points GDAL to the folder I created with the DTM
   
   
4. If you would like to get a list of the files in the folder you are working in, use this command: dir
   
   
5. To calculate TPI, I used the following command: gdaldem TPI DTM.tif TPI.tif
1. In this command, the words ‘gdaldem TPI’ define the utility that GDAL uses to process the data
2. ‘DTM.tif’ is the full name of the input (original) DEM file
3. ‘TPI.tif’ is the name I chose for the output topographic position map as a TIFF file; you can change the file extension ending to any in this list if you would prefer an alternate raster format
   
   
6. There is no extra functionality commands that can be added to TRI in GDAL

Tips for Interpreting and Working with TPI Maps in ArcGIS

As has been the trend with these WIG pieces, creating symbology in ArcGIS required to interpret an output from GDAL is a multi-step process. The first step in the process is determining the maximum and minimum range for TPI. I will follow the same general procedure set out in previous WIGs to determine this. As TPI has both negative and positive values, that will add an extra wrinkle to this process from what I presented last month. As such, I created a short video tutorial using ArcMAP 10.1 to show you the steps I use.

Now that you have seen how to determine your maximum and minimum values for TPI, it is time to move on to creating a color scheme that helps us locate fine geologic structures. To do this, you will want to right click on your TPI layer and go to Properties -> Symbology -> Classified. From here, click on the Set Exclusions… button and enter the correct expression to exclude the values above and below your true TPI range. What you are trying to do here is exclude the ‘fake’ TPI data that is created by GDAL at the borders of your raster. For my exclusions, I used this expression:

-16500 - -0.697, 0.725 - 20600

To create a color scheme that helps find fine features, it makes sense to use a Classified technique whereby ranges of TPI values are represented by easy to spot colors. You can do this from the same menu as the exclusions were entered into, i.e. Properties -> Symbology -> Classified. To do this, I ran multiple classifications offered in the pull down menu with 5 classes; and I found that Natural Breaks (Jenks) made the most intuitive sense with my data.

From there, I hit the Classify... button and changed the Break Value closest to zero to exactly 0.0. What this did is create two categories of TPI with positive values (i.e. hills) and then three with negative values (i.e. gullies). I then selected a simple color ramp of red to green so that the colors red, orange and yellow were negative; and then dark and light green were positive.

As a final step, I added high resolution imagery available through the ArcMap 10.x basemap menu to help bring the TPI map to life so to speak. The animation below shows the result of this color scheme and base layer imagery. By using this combination, I could identify thin ridges and gullies as well as roads that might have been missed by inspecting the high resolution imagery and/or elevation data on its own.

Brock Adam McCarty
Map Wizard
(720) 470-7988
Brock@apollomapping.com