Posted on January 10th, 2013

Working In GDAL (WIG) – Creating Contours

In this, the last of our monthly WIG series, I once again 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! For the past months, I focused 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 could emulate my steps in your office and/or home. In this final edition of WIG, I show our readers how to use GDAL to create contour vectors from a raster digital elevation model (DEM). While GDAL is a powerful tool, there are only so many features that warrant an explanation and unfortunately we have reached the end of the list with this edition of WIG.

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.x) 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 are Contours?

In past newsletters, we have covered this topic in some detail (you can read more here), so for this edition of WIG, I will leave it at a high-level introduction only. Contours were made popular by the green and beige United Stated Geological Survey maps many of us saw in high school/college geography classes. Simply put, a contour is a line that connects similar elevations on a map. Every map with contours has a contour interval which defines the elevation between each of the (contour) lines. So taken together, if a contour line has a 100-foot elevation and the contour interval is 25-feet, then the next contour line would be 125 feet and so forth from there. Typically, every fifth contour, or index line, is differentiate (with a wider line, for instance) for ease of reading.

Maps with contours on them create strong visual clues once users understand how to interpret them. When contours are closely spaced, they indicate steeper terrain; while those that are farther apart indicate flatter areas. A series of shrinking, closed circles shows the rise of a hilltop; while V’s indicate the presence of a stream and the direction it is flowing. Given the strong visual clues that are associated with vectors, they are a favorite of many outdoor enthusiasts, engineers and GIS-users alike.

Creating Contour Vectors

In order to create contours from a DEM, I employed the GDAL_CONTOUR function – you can find out more about it here. These are the steps I took to create the contours:

  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: cdtest
    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 color relief, I used the following command: gdal_contour –a elev DTM.tif contours.shp –inodata –snodata -32767 –i 5.0
    1. In this command, ‘gdal_contour’ defines the utility that GDAL uses to process the data
    2. ‘–a elev’ is used to create a new attribute in the output file for the elevation values; this is a crucial step as without it, you will get contours with no elevations tied to them
    3. ‘DTM.tif’ is the full name of the input (original) DEM file
    4. ‘contours.shp’ is the name I chose for the contour vectors
    5. ‘inodata’ tells GDAL to look for no data values (or black fill that is commonly used at the edge of DEM files)
    6. ‘-snodata -32767’ sets the no data value GDAL looks for; I will explain how I found this value below
    7. ‘-i 5.0’ creates a contour interval of 5-meters
  6. There are several other options you can add to the gdal_contours command that are worth looking at, for instance creating contours in a format other than shapefile
    1. A complete list of these commands can be found here

Working With Contours in ArcGIS

After creating contours in GDAL, let’s move over to ArcGIS and look at the data we created. When I was working in GDAL, I also created a 2-meter contour layer by altering the last option above (item g in the fifth step) to ‘-i 2.0’. A quick note before we move on. Since the data I am using is projected in UTM and the units are meters, when I use a value of 2.0 for the interval in GDAL, it creates a 2.0-meter contour. If your data were in geographic projection, the interval value would be equated as degrees and not meters, so 2 degrees which is a huge area!

The first step I actually took in preparing for this WIG was loading the DEM in ArcGIS to figure out what value was used at the edges for black fill or ‘no data’ values. This step is pretty straight forward – I will also show you how to find this no data value in the short movie that is linked below. Simply load your DEM and click on the blue circle in ArcGIS, also called the Identify Tool. Now click around the edges of your DEM and see what value is in each pixel; the edge values should be very distinct from the true elevation values. Most DEMs will have very large negative numbers for the no data values – in my DEM, the no data value is -32,767.

Now, let’s move on to symbology. When you load your contours into ArcGIS, you will start with a single symbol for your polyline file which can be a bit hard to interpret. Further, you will find that creating symbology that makes sense, for instance adding in index lines, is not an easy process in Arc 10.x. In fact, I spent a good time of time testing out scripts that would automatically create a nice symbology scheme for contours, but none of them worked in Arc 10.1. Therefore, here are some some tips on how to get up and running with your contours quickly.

First, adding labels for elevation heights that show up when you are zoomed in is a big help. To do this, you will want to right click on your contours and go to PropertiesLabels. Click the Label Features checkbox at the top of the window; and then change the ‘Label Field’ to your elevation attribute. Be sure to use a small font such as Arial 8, and you might consider adding a Mask if you have lots of other layers loaded at once. Now click on ‘Placement Properties’ and for ‘Placement’ select a Parallel orientation, On the Line position and remove duplicate labels. Under Conflict Detection, you may want to experiment with these settings if you have multiple layers loaded – I simply used the defaults. Finally, set the ‘Scale Range’ to Don’t show labels when zoomed: ‘Out beyond 1:10,000.’ This will prevent your display from being cluttered with impossible to read tiny labels. As a heads up, I will run through these steps in the movie below as well.

Now, let’s add index lines at elevations that make sense for a quick visual reference. I will show you this process in the movie below on a 2-meter contour where I add index lines at 100 meters, 110 meters and so forth. I have seen multiple approaches suggested for the creation of index lines and while the procedure I present in the movie might not be automated, it is very straightforward and quick to implement, hence my reliance upon it.

A Note on Interpolation Versus Over-Interpolation

contours_comparisonClick on the image above to see an animation comparing a 2-meter and 5-meter contour layer created over the same location using the same elevation data. In this animation, the 5-meter contours are shown in red and then the 2-meter contours are in blue. The DTM layer has been added below as a reference where green is the lowest elevations, yellow are mid and purple the highest. You can see the extra information that is added when 5-meter contours are overlain on the DTM; and then even more when the 2-meter contours are overlain. (DTM Source: DigitalGlobe)

When contours are created, you are connecting areas that are equal in elevation by using known or spot elevations. These known elevations will rarely fall exactly on a contour line you are trying to create, for instance at 100 meters, where the closest spot elevations might be 99.2 meters and 101.3 meters. So to create a 100-meter contour, there is a certain amount of interpolation (or estimation) to plot the exact position of this line. I do not make this point to discredit contour creation as they are certainly a valid form of elevation data; I make the point to state that it is possible to over-interpolate or to create information that is an over assumption (perhaps even a fabrication) based on the data provided to you. In the world of elevation data, I suggest you select a contour interval equal to or greater in value than the relative accuracy of your DEM.

Relative accuracy measures how close elevations within the same layer are from each other. For example, take two elevations on a hill, one at 10 meters and the other at 15 meters – the difference in height between these two points is 5 meters. If the relative accuracy of the DEM is 1 meters, that tells you that these points will have between a 4 and 6 meter height difference in the ‘real world;’ even if their absolute heights are much different than 10 and 15 meters.

As a concluding note to this series, if you have followed each of our editions of WIG, then you are armed with the tools you need to understand and utilize the other functionality offered in GDAL. I hope you enjoyed this series and learned a little bit along the way!

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

Share This Article
This entry was posted in The Geospatial Times and tagged , , , , , by Apollo Mapping. Bookmark the permalink.

Leave a Reply

Your email address will not be published. Required fields are marked *

Time limit is exhausted. Please reload CAPTCHA.

This site uses Akismet to reduce spam. Learn how your comment data is processed.

    The Geospatial Times Archive